betabinomial.RdFits a beta-binomial distribution by maximum likelihood estimation. The two parameters here are the mean and correlation coefficient.
betabinomial(lmu = "logitlink", lrho = "logitlink",
irho = NULL, imethod = 1,
ishrinkage = 0.95, nsimEIM = NULL, zero = "rho")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.
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.
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.
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.
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.
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.
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\).
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.
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.
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.
# 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