View source: R/family.binomial.R
seq2binomial | R Documentation |
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 |
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 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) =
{m \choose y_1} p^{y_1} (1-p)^{m-y_1}
{y_1 \choose 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.
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 y_1
and y_2
.
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 y_1
,
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.
binomialff
,
cfibrosis
.
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])
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.