seq2binomial: The Two-stage Sequential Binomial Distribution Family...

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

View source: R/family.binomial.R

Description

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

Usage

1
2
3
seq2binomial(lprob1 = "logitlink", lprob2 = "logitlink",
             iprob1 = NULL,    iprob2 = NULL,
             parallel = FALSE, zero = NULL)

Arguments

lprob1, lprob2

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

iprob1, iprob2

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

parallel, zero

Details at Links. If parallel = TRUE then the constraint also applies to the intercept.

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

Value

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

Note

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.

Author(s)

Thomas W. Yee

References

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

See Also

binomialff, cfibrosis.

Examples

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

Example output

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 

VGAM documentation built on Jan. 16, 2021, 5:21 p.m.