Zero-Inflated Binomial Distribution Family Function

Fits a zero-inflated binomial distribution by maximum likelihood estimation.

zibinomial(lpstr0 = "logitlink", lprob = "logitlink",
           type.fitted = c("mean", "prob", "pobs0", "pstr0", "onempstr0"),
           ipstr0 = NULL, zero = NULL, multiple.responses = FALSE,
           imethod = 1)
zibinomialff(lprob = "logitlink", lonempstr0 = "logitlink",
             type.fitted = c("mean", "prob", "pobs0", "pstr0", "onempstr0"),
             ionempstr0 = NULL, zero = "onempstr0",
             multiple.responses = FALSE, imethod = 1)

Arguments

lpstr0, lprob

Link functions for the parameter $\phi$ and the usual binomial probability $\mu$ parameter. See Links for more choices. For the zero-deflated model see below.

type.fitted

See CommonVGAMffArguments and fittedvlm.

ipstr0

Optional initial values for $\phi$, whose values must lie between 0 and 1. The default is to compute an initial value internally. If a vector then recyling is used.

lonempstr0, ionempstr0

Corresponding arguments for the other parameterization. See details below.

multiple.responses

Logical. Currently it must be FALSE to mean the function does not handle multiple responses. This is to remain compatible with the same argument in binomialff.

zero, imethod

See CommonVGAMffArguments for information. Argument zero changed its default value for version 0.9-2.

Details

These functions are based on $$P(Y=0) = \phi + (1-\phi) (1-\mu)^N,$$ for $y=0$, and $$P(Y=y) = (1-\phi) {N \choose Ny} \mu^{Ny} (1-\mu)^{N(1-y)}.$$ for $y=1/N,2/N,\ldots,1$. That is, the response is a sample proportion out of $N$ trials, and the argument size in rzibinom is $N$ here. The parameter $\phi$ is the probability of a structural zero, and it satisfies $0 < \phi < 1$. The mean of $Y$ is $E(Y)=(1-\phi) \mu$ and these are returned as the fitted values by default. By default, the two linear/additive predictors for zibinomial() are $(logit(\phi), logit(\mu))^T$.

The VGAM family function zibinomialff() has a few changes compared to zibinomial(). These are: (i) the order of the linear/additive predictors is switched so the binomial probability comes first; (ii) argument onempstr0 is now 1 minus the probability of a structural zero, i.e., the probability of the parent (binomial) component, i.e., onempstr0 is 1-pstr0; (iii) argument zero has a new default so that the onempstr0 is intercept-only by default. Now zibinomialff() is generally recommended over zibinomial(). Both functions implement Fisher scoring.

Value

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

References

Welsh, A. H., Lindenmayer, D. B. and Donnelly, C. F. (2013). Fitting and interpreting occupancy models. PLOS One, 8, 1–21.

Author

T. W. Yee

Note

The response variable must have one of the formats described by binomialff, e.g., a factor or two column matrix or a vector of sample proportions with the weights argument specifying the values of $N$.

To work well, one needs large values of $N$ and $\mu>0$, i.e., the larger $N$ and $\mu$ are, the better. If $N = 1$ then the model is unidentifiable since the number of parameters is excessive.

Setting stepsize = 0.5, say, may aid convergence.

Estimated probabilities of a structural zero and an observed zero are returned, as in zipoisson.

The zero-deflated binomial distribution might be fitted by setting lpstr0 = identitylink, albeit, not entirely reliably. See zipoisson for information that can be applied here. Else try the zero-altered binomial distribution (see zabinomial).

Warning

Numerical problems can occur. Half-stepping is not uncommon. If failure to converge occurs, make use of the argument ipstr0 or ionempstr0, or imethod.

Examples

size <- 10  # Number of trials; N in the notation above
nn <- 200
zdata <- data.frame(pstr0 = logitlink( 0, inverse = TRUE),  # 0.50
                    mubin = logitlink(-1, inverse = TRUE),  # Mean of usual binomial
                    sv    = rep(size, length = nn))
zdata <- transform(zdata,
                   y = rzibinom(nn, size = sv, prob = mubin, pstr0 = pstr0))
with(zdata, table(y))
#> y
#>   0   1   2   3   4   5   6 
#> 106  14  25  27  20   5   3 
fit <- vglm(cbind(y, sv - y) ~ 1, zibinomialff, data = zdata, trace = TRUE)
#> Iteration 1: loglikelihood = -347.32035
#> Iteration 2: loglikelihood = -383.10832
#> Taking a modified step.
#> Iteration  2 :  loglikelihood = -302.03626
#> Iteration 3: loglikelihood = -291.05676
#> Iteration 4: loglikelihood = -289.20817
#> Iteration 5: loglikelihood = -289.20714
#> Iteration 6: loglikelihood = -289.20714
fit <- vglm(cbind(y, sv - y) ~ 1, zibinomialff, data = zdata, trace = TRUE,
            stepsize = 0.5)
#> Taking a modified step.
#> Iteration  2 :  loglikelihood = -302.03626
#> Taking a modified step.
#> Iteration  3 :  loglikelihood = -290.91281
#> Taking a modified step.
#> Iteration  4 :  loglikelihood = -289.56623
#> Taking a modified step.
#> Iteration  5 :  loglikelihood = -289.29091
#> Taking a modified step.
#> Iteration  6 :  loglikelihood = -289.22744
#> Taking a modified step.
#> Iteration  7 :  loglikelihood = -289.21214
#> Taking a modified step.
#> Iteration  8 :  loglikelihood = -289.20838
#> Taking a modified step.
#> Iteration  9 :  loglikelihood = -289.20745
#> Taking a modified step.
#> Iteration  10 :  loglikelihood = -289.20722
#> Taking a modified step.
#> Iteration  11 :  loglikelihood = -289.20716
#> Taking a modified step.
#> Iteration  12 :  loglikelihood = -289.20715
#> Taking a modified step.
#> Iteration  13 :  loglikelihood = -289.20714
#> Taking a modified step.
#> Iteration  14 :  loglikelihood = -289.20714

coef(fit, matrix = TRUE)
#>             logitlink(prob) logitlink(onempstr0)
#> (Intercept)      -0.9774094          -0.03968286
Coef(fit)  # Useful for intercept-only models
#>      prob onempstr0 
#> 0.2734061 0.4900806 
head(fitted(fit, type = "pobs0"))  # Estimate of P(Y = 0)
#>        [,1]
#> 1 0.5300189
#> 2 0.5300189
#> 3 0.5300189
#> 4 0.5300189
#> 5 0.5300189
#> 6 0.5300189
head(fitted(fit))
#>       [,1]
#> 1 0.133991
#> 2 0.133991
#> 3 0.133991
#> 4 0.133991
#> 5 0.133991
#> 6 0.133991
with(zdata, mean(y))  # Compare this with fitted(fit)
#> [1] 1.34
summary(fit)
#> 
#> Call:
#> vglm(formula = cbind(y, sv - y) ~ 1, family = zibinomialff, data = zdata, 
#>     trace = TRUE, stepsize = 0.5)
#> 
#> Coefficients: 
#>               Estimate Std. Error z value Pr(>|z|)    
#> (Intercept):1 -0.97741    0.07824 -12.493   <2e-16 ***
#> (Intercept):2 -0.03968    0.14835  -0.267    0.789    
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Names of linear predictors: logitlink(prob), logitlink(onempstr0)
#> 
#> Log-likelihood: -289.2071 on 398 degrees of freedom
#> 
#> Number of Fisher scoring iterations: 14 
#> 
#> No Hauck-Donner effect found in any of the estimates
#>