View source: R/family.binomial.R
betabinomial | R Documentation |
Fits 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")
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 \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
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
.
# 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)
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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.