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)
```

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.