Description Usage Arguments Details Value Warning Note Author(s) References See Also Examples

View source: R/family.binomial.R

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

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

`lmu, lrho` |
Link functions applied to the two parameters.
See |

`irho` |
Optional initial value for the correlation parameter.
If given, it must be in |

`imethod` |
An integer with value |

`zero` |
Specifyies which
linear/additive predictor is to be modelled as an intercept only.
If assigned, the single value can be either |

`ishrinkage, nsimEIM` |
See |

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 ~ 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) = choose(N,t) Be(alpha+t, beta+N-t) / Be(alpha, beta)*

where *t=0,1,…,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 *eta1 =logit(mu)*
and *eta2 = 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
`fit@y`

contains the sample proportions *y*,
`fitted(fit)`

returns estimates of *E(Y)*, and
`weights(fit, type="prior")`

returns the number
of trials *N*.

If the estimated rho parameter is close to 0 then it pays to
try `lrho = "rhobitlink"`

. One day this may become the default
link function.

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.

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)
is currently not implemented in the VGAM package as it has
range-restrictions for the correlation parameter that are currently
too difficult to handle in this package.
However, try `lrho = "rhobitlink"`

.

T. W. Yee

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.

Prentice, R. L. (1986).
Binary regression using an extended beta-binomial distribution,
with discussion of correlation induced by
covariate measurement errors.
*Journal of the American Statistical Association*,
**81**, 321–327.

`betabinomialff`

,
`Betabinom`

,
`binomialff`

,
`betaff`

,
`dirmultinomial`

,
`lirat`

,
`simulate.vlm`

.

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 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 | ```
# 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, data = bdata, trace = TRUE)
coef(fit, matrix = TRUE)
Coef(fit)
head(cbind(depvar(fit), weights(fit, type = "prior")))
# Example 2
fit <- vglm(cbind(R, N-R) ~ 1, betabinomial, lirat,
trace = TRUE, subset = N > 1)
coef(fit, matrix = TRUE)
Coef(fit)
t(fitted(fit))
t(depvar(fit))
t(weights(fit, type = "prior"))
# Example 3, which is more complicated
lirat <- transform(lirat, fgrp = factor(grp))
summary(lirat) # Only 5 litters in group 3
fit2 <- vglm(cbind(R, N-R) ~ fgrp + hb, betabinomial(zero = 2),
data = lirat, trace = TRUE, subset = N > 1)
coef(fit2, matrix = TRUE)
## Not run: with(lirat, plot(hb[N > 1], fit2@misc$rho,
xlab = "Hemoglobin", ylab = "Estimated rho",
pch = as.character(grp[N > 1]), col = grp[N > 1]))
## End(Not run)
## Not run: # 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)
}
## End(Not run)
``` |

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.