Beta-binomial Distribution Family Function

Fits a beta-binomial distribution by maximum likelihood estimation. The two parameters here are the mean and correlation coefficient.

Usage

betabinomial(lmu = "logitlink", lrho = "logitlink",
   irho = NULL, imethod = 1,
   ishrinkage = 0.95, nsimEIM = NULL, zero = "rho")

Arguments

lmu, lrho

Link functions applied to the two parameters. See Links for more choices. The defaults ensure the parameters remain in $(0,1)$, however, see the warning below. For lrho, log1plink (with an offset log(size - 1) for $\eta_2$) and cloglink may be very good choices.

irho

Optional initial value for the correlation parameter. If given, it must be in $(0,1)$, and is recyled to the necessary length. Assign this argument a value if a convergence failure occurs. Having irho = NULL means an initial value is obtained internally, though this can give unsatisfactory results.

imethod

An integer with value 1 or 2 or ..., which specifies the initialization method for $\mu$. If failure to converge occurs try the another value and/or else specify a value for irho.

zero

Specifies which linear/additive predictor is to be modelled as an intercept only. If assigned, the single value can be either 1 or 2. The default is to have a single correlation parameter. To model both parameters as functions of the covariates assign zero = NULL. See CommonVGAMffArguments for more information.

ishrinkage, nsimEIM

See CommonVGAMffArguments for more information. The argument ishrinkage is used only if imethod = 2. Using the argument nsimEIM may offer large advantages for large values of $N$ and/or large data sets.

Details

There are several parameterizations of the beta-binomial distribution. This family function directly models the mean and correlation parameter, i.e., the probability of success. The model can be written $T|P=p \sim Binomial(N,p)$ where $P$ has a beta distribution with shape parameters $\alpha$ and $\beta$. Here, $N$ is the number of trials (e.g., litter size), $T=NY$ is the number of successes, and $p$ is the probability of a success (e.g., a malformation). That is, $Y$ is the proportion of successes. Like binomialff, the fitted values are the estimated probability of success (i.e., $E[Y]$ and not $E[T]$) and the prior weights $N$ are attached separately on the object in a slot.

The probability function is $$P(T=t) = {N \choose t} \frac{Be(\alpha+t, \beta+N-t)} {Be(\alpha, \beta)}$$ where $t=0,1,\ldots,N$, and $Be$ is the beta function with shape parameters $\alpha$ and $\beta$. Recall $Y = T/N$ is the real response being modelled.

The default model is $\eta_1 = logit(\mu)$ and $\eta_2 = logit(\rho)$ because both parameters lie between 0 and 1. The mean (of $Y$) is $p=\mu=\alpha/(\alpha+\beta)$ and the variance (of $Y$) is $\mu(1-\mu)(1+(N-1)\rho)/N$. Here, the correlation $\rho$ is given by $1/(1 + \alpha + \beta)$ and is the correlation between the $N$ individuals within a litter. A litter effect is typically reflected by a positive value of $\rho$. It is known as the over-dispersion parameter.

This family function uses Fisher scoring. Elements of the second-order expected derivatives with respect to $\alpha$ and $\beta$ are computed numerically, which may fail for large $\alpha$, $\beta$, $N$ or else take a long time.

Value

An object of class "vglmff" (see vglmff-class). The object is used by modelling functions such as vglm.

Suppose fit is a fitted beta-binomial model. Then depvar(fit) are the sample proportions $y$, fitted(fit) returns estimates of $E(Y)$, and weights(fit, type = "prior") returns the number of trials $N$.

References

Moore, D. F. and Tsiatis, A. (1991). Robust estimation of the variance in moment methods for extra-binomial and extra-Poisson variation. Biometrics, 47, 383–401.

Author

T. W. Yee

Note

This function processes the input in the same way as binomialff. But it does not handle the case $N=1$ very well because there are two parameters to estimate, not one, for each row of the input. Cases where $N=1$ can be omitted via the subset argument of vglm.

The extended beta-binomial distribution of Prentice (1986) implemented by extbetabinomial is the preferred VGAM family function for BBD regression.

Warning

If the estimated rho parameter is close to 0 then a good solution is to use extbetabinomial. Or you could try lrho = "rhobitlink".

This family function is prone to numerical difficulties due to the expected information matrices not being positive-definite or ill-conditioned over some regions of the parameter space. If problems occur try setting irho to some numerical value, nsimEIM = 100, say, or else use etastart argument of vglm, etc.

Examples

# Example 1
bdata <- data.frame(N = 10, mu = 0.5, rho = 0.8)
bdata <- transform(bdata,
            y = rbetabinom(100, size = N, prob = mu, rho = rho))
fit <- vglm(cbind(y, N-y) ~ 1, betabinomial, bdata, trace = TRUE)
#> Iteration 1: loglikelihood = -181.28351
#> Iteration 2: loglikelihood = -181.28342
#> Iteration 3: loglikelihood = -181.28342
coef(fit, matrix = TRUE)
#>             logitlink(mu) logitlink(rho)
#> (Intercept)      0.123523       1.127352
Coef(fit)
#>        mu       rho 
#> 0.5308416 0.7553500 
head(cbind(depvar(fit), weights(fit, type = "prior")))
#>   [,1] [,2]
#> 1  0.4   10
#> 2  0.0   10
#> 3  0.0   10
#> 4  0.2   10
#> 5  0.5   10
#> 6  1.0   10


# Example 2
fit <- vglm(cbind(R, N-R) ~ 1, betabinomial, lirat,
            trace = TRUE, subset = N > 1)
#> Iteration 1: loglikelihood = -122.69536
#> Iteration 2: loglikelihood = -122.69531
#> Iteration 3: loglikelihood = -122.69531
coef(fit, matrix = TRUE)
#>             logitlink(mu) logitlink(rho)
#> (Intercept)    -0.1190049      0.4069599
Coef(fit)
#>        mu       rho 
#> 0.4702838 0.6003587 
t(fitted(fit))
#>           [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
#> [1,] 0.4702838 0.4702838 0.4702838 0.4702838 0.4702838 0.4702838 0.4702838
#>           [,8]      [,9]     [,10]     [,11]     [,12]     [,13]     [,14]
#> [1,] 0.4702838 0.4702838 0.4702838 0.4702838 0.4702838 0.4702838 0.4702838
#>          [,15]     [,16]     [,17]     [,18]     [,19]     [,20]     [,21]
#> [1,] 0.4702838 0.4702838 0.4702838 0.4702838 0.4702838 0.4702838 0.4702838
#>          [,22]     [,23]     [,24]     [,25]     [,26]     [,27]     [,28]
#> [1,] 0.4702838 0.4702838 0.4702838 0.4702838 0.4702838 0.4702838 0.4702838
#>          [,29]     [,30]     [,31]     [,32]     [,33]     [,34]     [,35]
#> [1,] 0.4702838 0.4702838 0.4702838 0.4702838 0.4702838 0.4702838 0.4702838
#>          [,36]     [,37]     [,38]     [,39]     [,40]     [,41]     [,42]
#> [1,] 0.4702838 0.4702838 0.4702838 0.4702838 0.4702838 0.4702838 0.4702838
#>          [,43]     [,44]     [,45]     [,46]     [,47]     [,48]     [,49]
#> [1,] 0.4702838 0.4702838 0.4702838 0.4702838 0.4702838 0.4702838 0.4702838
#>          [,50]     [,51]     [,52]     [,53]     [,54]     [,55]     [,56]
#> [1,] 0.4702838 0.4702838 0.4702838 0.4702838 0.4702838 0.4702838 0.4702838
#>          [,57]
#> [1,] 0.4702838
t(depvar(fit))
#>        1         2    3 4 5         6 7 8 9  10 11  12 13        14        15
#> [1,] 0.1 0.3636364 0.75 1 1 0.8181818 1 1 1 0.7  1 0.9  1 0.8181818 0.6666667
#>             16 17        18        19        20        21 22        23
#> [1,] 0.7777778  1 0.5833333 0.8181818 0.6153846 0.3571429  1 0.8333333
#>             24 25        26 27   28 29        30 31  32        33         34 35
#> [1,] 0.6153846  1 0.2142857  1 0.75  1 0.3846154  1 0.1 0.3333333 0.07692308  0
#>             36        37        38     39 40 41 43 44         45 46         47
#> [1,] 0.2857143 0.2222222 0.1538462 0.0625  0  0  0  0 0.09090909  0 0.07142857
#>      48 49 50        51        52 53 54         55 56 57 58
#> [1,]  0  0  0 0.2222222 0.1176471  0  0 0.07142857  0  0  0
t(weights(fit, type = "prior"))
#>       1  2  3 4  5  6 7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
#> [1,] 10 11 12 4 10 11 9 11 10 10 12 10  8 11  6  9 14 12 11 13 14 10 12 13 10
#>      26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 43 44 45 46 47 48 49 50 51
#> [1,] 14 13  4  8 13 12 10  3 13 12 14  9 13 16 11  4 12  8 11 14 14 11  3 13  9
#>      52 53 54 55 56 57 58
#> [1,] 17 15  2 14  8  6 17


# Example 3, which is more complicated
lirat <- transform(lirat, fgrp = factor(grp))
summary(lirat)  # Only 5 litters in group 3
#>        N               R                hb              grp        fgrp  
#>  Min.   : 1.00   Min.   : 0.000   Min.   : 3.100   Min.   :1.000   1:31  
#>  1st Qu.: 9.00   1st Qu.: 0.250   1st Qu.: 4.325   1st Qu.:1.000   2:12  
#>  Median :11.00   Median : 3.500   Median : 6.500   Median :1.000   3: 5  
#>  Mean   :10.47   Mean   : 4.603   Mean   : 7.798   Mean   :1.897   4:10  
#>  3rd Qu.:13.00   3rd Qu.: 9.000   3rd Qu.:10.775   3rd Qu.:2.750         
#>  Max.   :17.00   Max.   :14.000   Max.   :16.600   Max.   :4.000         
fit2 <- vglm(cbind(R, N-R) ~ fgrp + hb, betabinomial(zero = 2),
             data = lirat, trace = TRUE, subset = N > 1)
#> Iteration 1: loglikelihood = -104.03084
#> Iteration 2: loglikelihood = -94.204685
#> Iteration 3: loglikelihood = -93.00179
#> Iteration 4: loglikelihood = -92.911711
#> Iteration 5: loglikelihood = -92.908006
#> Iteration 6: loglikelihood = -92.907857
#> Iteration 7: loglikelihood = -92.907851
#> Iteration 8: loglikelihood = -92.907851
coef(fit2, matrix = TRUE)
#>             logitlink(mu) logitlink(rho)
#> (Intercept)     2.0895865      -1.171752
#> fgrp2          -2.4529617       0.000000
#> fgrp3          -2.8840242       0.000000
#> fgrp4          -2.3637974       0.000000
#> hb             -0.1609718       0.000000
if (FALSE)  with(lirat, plot(hb[N > 1], fit2@misc$rho,
         xlab = "Hemoglobin", ylab = "Estimated rho",
         pch = as.character(grp[N > 1]), col = grp[N > 1]))  # \dontrun{}
if (FALSE)   # cf. Figure 3 of Moore and Tsiatis (1991)
with(lirat, plot(hb, R / N, pch = as.character(grp), col = grp,
         xlab = "Hemoglobin level", ylab = "Proportion Dead",
         main = "Fitted values (lines)", las = 1))
smalldf <- with(lirat, lirat[N > 1, ])
for (gp in 1:4) {
  xx <- with(smalldf, hb[grp == gp])
  yy <- with(smalldf, fitted(fit2)[grp == gp])
  ooo <- order(xx)
  lines(xx[ooo], yy[ooo], col = gp, lwd = 2)
}  # \dontrun{}
#> Error in plot.xy(xy.coords(x, y), type = type, ...): plot.new has not been called yet