The Two-stage Sequential Binomial Distribution Family Function
Estimation of the probabilities of a two-stage binomial distribution.
seq2binomial(lprob1 = "logitlink", lprob2 = "logitlink",
iprob1 = NULL, iprob2 = NULL,
parallel = FALSE, zero = NULL)lprob1, lprob2 |
Parameter link functions applied to the two probabilities,
called p and q below.
See |
iprob1, iprob2 |
Optional initial value for the first and second probabilities respectively.
A |
parallel, zero |
Details at |
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 y1 spores that germinate, each has a probability q of bending in a particular direction. Let y2 be the number that bend in the specified direction. The probability model for this data is P(y1,y2) =
{choose(m,y1)} p^{y1} (1-p)^{m-y1} {choose(y1,y2)} q^{y2} (1-q)^{y1-y2}
for 0 < p < 1, 0 < q < 1,
y1=1,…,m
and
y2=1,…,y1.
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.
An object of class "vglmff" (see vglmff-class).
The object is used by modelling functions such as vglm
and vgam.
The response must be a two-column matrix of sample proportions
corresponding to y1 and y2.
The m values should be inputted with the weights
argument of vglm
and vgam.
The fitted value is a two-column matrix of estimated probabilities
p and q.
A common form of error is when there are no trials
for y1, e.g., if mvector below has some values
which are zero.
Thomas W. Yee
Crowder, M. and Sweeting, T. (1989). Bayesian inference for a bivariate binomial distribution. Biometrika, 76, 599–603.
sdata <- data.frame(mvector = round(rnorm(nn <- 100, m = 10, sd = 2)),
x2 = runif(nn))
sdata <- transform(sdata, prob1 = logitlink(+2 - x2, inverse = TRUE),
prob2 = logitlink(-2 + x2, 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) ~ x2, seq2binomial, weight = mvector,
data = sdata, trace = TRUE)
coef(fit)
coef(fit, matrix = TRUE)
head(fitted(fit))
head(depvar(fit))
head(weights(fit, type = "prior")) # Same as with(sdata, mvector)
# Number of first successes:
head(depvar(fit)[, 1] * c(weights(fit, type = "prior")))
# Number of second successes:
head(depvar(fit)[, 2] * c(weights(fit, type = "prior")) *
depvar(fit)[, 1])Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.