seq2binomial: The Two-stage Sequential Binomial Distribution Family Function

Description

Estimation of the probabilities of a two-stage binomial distribution.

Usage

seq2binomial(lprob1 = "logit", lprob2 = "logit", eprob1 = list(),
             eprob2 = list(), iprob1 = NULL, iprob2 = NULL, zero = NULL)

Arguments

lprob1, lprob2

Parameter link functions applied to the two probabilities, called $p$ and $q$ below. See Links for more choices.

eprob1, eprob2

Lists. Extra arguments for the links. See earg in Links for general information.

iprob1, iprob2

Optional initial value for the first and second probabilities respectively. A NULL means a value is obtained in the initialize slot.

zero

An integer-valued vector specifying which linear/additive predictors are modelled as intercepts only. If used, the value must be from the set {1,2} which correspond to the first and second probabilities respectively. A NULL value mean

Value

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

Details

This VGAM family function fits the model described by Crowder and Sweeting (1989) which is described as follows. Each of $m$ spores has a probability $p$ of germinating. Of the $y_1$ spores that germinate, each has a probability $q$ of bending in a particular direction. Let $y_2$ be the number that bend in the specified direction. The probability model for this data is $P(y_1,y_2) =$

(\binom{m}{y_{1}}) p^{y_{1}} (1 - p)^{m - y_{1}} (\binom{y_{1}}{y_{2}}) q^{y_{2}} (1 - q)^{y_{1} - y_{2}}

for $0 < p < 1$, $0 < q < 1$, $y_1=1,\ldots,m$ and $y_2=1,\ldots,y_1$. Here, $p$ is prob1, $q$ is prob2.

Although the Authors refer to this as the bivariate binomial model, I have named it the (two-stage) sequential binomial model. Fisher scoring is used.

References

Crowder, M. and Sweeting, T. (1989). Bayesian inference for a bivariate binomial distribution. Biometrika, 76, 599--603.

Examples

Run this code

sdata = data.frame(mvector = round(rnorm(nn <- 100, m = 10, sd = 2)),
                   x = runif(nn))
sdata = transform(sdata, prob1 = logit(+2 - x, inverse = TRUE),
                         prob2 = logit(-2 + x, inverse = TRUE))
sdata = transform(sdata, successes1 = rbinom(nn, size=mvector, prob=prob1))
sdata = transform(sdata, successes2 = rbinom(nn, size=successes1, prob=prob2))
sdata = transform(sdata, y1 = successes1 / mvector)
sdata = transform(sdata, y2 = successes2 / successes1)
fit = vglm(cbind(y1,y2) ~ x, seq2binomial,  weight=mvector,
           data = sdata, trace=TRUE)
coef(fit)
coef(fit, matrix = TRUE)
head(fitted(fit))

Run the code above in your browser using DataLab

Last chance! 50% off unlimited learning