Description Usage Arguments Details Value Note Author(s) References See Also Examples

View source: R/family.binomial.R

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

1 2 3 |

`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 *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.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 | ```
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])
``` |

```
Loading required package: stats4
Loading required package: splines
VGLM linear loop 1 : loglikelihood = -302.14583
VGLM linear loop 2 : loglikelihood = -301.97106
VGLM linear loop 3 : loglikelihood = -301.97102
VGLM linear loop 4 : loglikelihood = -301.97102
(Intercept):1 (Intercept):2 x2:1 x2:2
2.065922 -2.335380 -1.265935 1.295785
logitlink(prob1) logitlink(prob2)
(Intercept) 2.065922 -2.335380
x2 -1.265935 1.295785
prob1 prob2
1 0.8843954 0.09083758
2 0.8351175 0.13223533
3 0.8642996 0.10756377
4 0.8458077 0.12314597
5 0.7802782 0.17976885
6 0.7076150 0.24500856
y1 y2
1 1.0000000 0.1000000
2 1.0000000 0.3333333
3 1.0000000 0.2500000
4 1.0000000 0.3636364
5 0.7142857 0.2000000
6 0.8000000 0.3750000
[,1]
[1,] 10
[2,] 9
[3,] 8
[4,] 11
[5,] 7
[6,] 10
1 2 3 4 5 6
10 9 8 11 5 8
1 2 3 4 5 6
1 3 2 4 1 3
```

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.