zibinomial: Zero-Inflated Binomial Distribution Family Function

Description

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

Usage

zibinomial(lpstr0 = "logit", lprob = "logit",
           ipstr0 = NULL, zero = 1, mv = 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.

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.

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

zero, imethod

See CommonVGAMffArguments for information.

Value

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

Warning

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

Details

This function uses Fisher scoring and is based on

P (Y = 0) = ϕ + (1 - ϕ) (1 - μ)^{N},

for $y=0$, and

P (Y = y) = (1 - ϕ) (\binom{N}{N y}) μ^{N y} (1 - μ)^{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, the two linear/additive predictors are $(logit(\phi), logit(\mu))^T$.

Examples

Run this code

size <- 10  # Number of trials; N in the notation above
nn <- 200
zibdata <- data.frame(pstr0 = logit( 0, inverse = TRUE), # 0.50
                      mubin = logit(-1, inverse = TRUE), # Mean of usual binomial
                      sv    = rep(size, length = nn))
zibdata <- transform(zibdata, 
                     y = rzibinom(nn, size = sv, prob = mubin, pstr0 = pstr0))
with(zibdata, table(y))
fit <- vglm(cbind(y, sv - y) ~ 1, zibinomial, zibdata, trace = TRUE)
fit <- vglm(cbind(y, sv - y) ~ 1, zibinomial, zibdata, trace = TRUE, stepsize = 0.5)

coef(fit, matrix = TRUE)
Coef(fit)  # Useful for intercept-only models
fit@misc$pobs0  # Estimate of P(Y = 0)
head(fitted(fit))
with(zibdata, mean(y))  # Compare this with fitted(fit)
summary(fit)

Run the code above in your browser using DataLab

Get 50% off unlimited learning