# Zero-Inflated Model

Log likelihood function to estimate parameters for a Zero-inflated Poisson model. With examples and comparison to pscl package output. Also includes approach based on Hilbe GLM text.

Zero-inflated models are applied to situations in which target data has relatively many of one value, usually zero, to go along with the other observed values. They are two-part models, a logistic model for whether an observation is zero or not, and a count model for the other part. The key distinction from hurdle count models is that the count distribution contributes to the excess zeros. While the typical application is count data, the approach can be applied to any distribution in theory.

## Poisson

### Data Setup

Get the data.

library(tidyverse)

fish = haven::read_dta("http://www.stata-press.com/data/r11/fish.dta")

### Function

The log likelihood function.

zip_ll <- function(y, X, par) {
# arguments are response y, predictor matrix X, and parameter named starting points of 'logit' and 'pois'

# Extract parameters
logitpars = par[grep('logit', names(par))]
poispars  = par[grep('pois', names(par))]

# Logit part; in this function Xlogit = Xpois but one could split X argument into Xlogi and Xpois for example
Xlogit  = X
LPlogit = Xlogit %*% logitpars
logi0   = plogis(LPlogit)  # alternative 1/(1+exp(-LPlogit))

# Poisson part
Xpois  = X
mupois = exp(Xpois %*% poispars)

# LLs
logliklogit = log( logi0 + exp(log(1 - logi0) - mupois) )
loglikpois  = log(1 - logi0) + dpois(y, lambda = mupois, log = TRUE)

# Hilbe formulation
# logliklogit = log(logi0 + (1 - logi0)*exp(- mupois) )
# loglikpois = log(1-logi0) -mupois + log(mupois)*y     #not necessary: - log(gamma(y+1))

y0 = y == 0  # 0 values
yc = y > 0   # Count part

loglik = sum(logliklogit[y0]) + sum(loglikpois[yc])
-loglik
}

### Estimation

Get starting values or simply do zeros.

# for zip: need 'logit', 'pois'
initial_model = glm(
count ~ persons + livebait,
data = fish,
x = TRUE,
y = TRUE,
"poisson"
)

# starts = c(logit = coef(initial_model), pois = coef(initial_model))
starts = c(rep(0, 3), rep(0, 3))

names(starts) = c(paste0('pois.', names(coef(initial_model))),
paste0('logit.', names(coef(initial_model))))

Estimate with optim.

fit = optim(
par = starts ,
fn  = zip_ll,
X   = initial_model$x, y = initial_model$y,
method  = "BFGS",
control = list(maxit = 5000, reltol = 1e-12),
hessian = TRUE
)

# fit

### Comparison

Extract for clean display.

B  = fit$par se = sqrt(diag(solve((fit$hessian))))
Z  = B/se
p  = pnorm(abs(Z), lower = FALSE)*2

Results from pscl.

library(pscl)

fit_pscl = zeroinfl(count ~ persons + livebait, data = fish, dist = "poisson") 

Compare.

coef B se Z p
pscl X.Intercept. -2.006 0.324 -6.196 0.000
pscl persons 0.747 0.043 17.516 0.000
pscl livebait 1.809 0.292 6.195 0.000
pscl X.Intercept..1 0.303 0.674 0.449 0.654
pscl persons.1 -0.069 0.129 -0.537 0.591
pscl livebait.1 -0.031 0.558 -0.056 0.956
zip_poisson_ll pois.(Intercept) -2.006 0.324 -6.196 0.000
zip_poisson_ll pois.persons 0.747 0.043 17.516 0.000
zip_poisson_ll pois.livebait 1.809 0.292 6.195 0.000
zip_poisson_ll logit.(Intercept) 0.303 0.674 0.449 0.654
zip_poisson_ll logit.persons -0.069 0.129 -0.537 0.591
zip_poisson_ll logit.livebait -0.031 0.558 -0.056 0.956

## Negative Binomial

### Function

zinb_ll <- function(y, X, par) {
# arguments are response y, predictor matrix X, and parameter named starting points of 'logit', 'negbin', and 'theta'

# Extract parameters
logitpars  = par[grep('logit', names(par))]
negbinpars = par[grep('negbin', names(par))]
theta = exp(par[grep('theta', names(par))])

# Logit part; in this function Xlogit = Xnegbin but one could split X argument into Xlogit and Xnegbin for example
Xlogit  = X
LPlogit = Xlogit %*% logitpars
logi0   = plogis(LPlogit)

# Negbin part
Xnegbin = X
munb = exp(Xnegbin %*% negbinpars)

# LLs
logliklogit  =
log(
logi0 +
exp(
log(1 - logi0) +
suppressWarnings(dnbinom(0, size = theta, mu   = munb, log  = TRUE))
)
)

logliknegbin = log(1 - logi0) +
suppressWarnings(dnbinom(y,
size = theta,
mu   = munb,
log  = TRUE))

# Hilbe formulation
# theta part
# alpha = 1/theta
# m = 1/alpha
# p = 1/(1 + alpha*munb)

# logliklogit = log( logi0 + (1 - logi0)*(p^m) )
# logliknegbin = log(1-logi0) + log(gamma(m+y)) - log(gamma(m)) + m*log(p) + y*log(1-p)   # gamma(y+1) not needed

y0 = y == 0   # 0 values
yc = y > 0    # Count part

loglik = sum(logliklogit[y0]) + sum(logliknegbin[yc])
-loglik
}

### Estimation

Get starting values or simply do zeros.

# for zinb: 'logit', 'negbin', 'theta'
initial_model = model.matrix(count ~ persons + livebait, data = fish) # to get X matrix

startlogi  = glm(count == 0 ~ persons + livebait, data = fish, family = "binomial")
startcount = glm(count ~ persons + livebait, data = fish, family = "poisson")

starts = c(
negbin = coef(startcount),
logit = coef(startlogi),
theta = 1
)
# starts = c(negbin = rep(0, 3),
#            logit = rep(0, 3),
#            theta = log(1))

Estimate with optim.

fit_nb = optim(
par = starts ,
fn  = zinb_ll,
X   = initial_model,
se = sqrt(diag(solve((fit_nb$hessian)))) Z = B/se p = pnorm(abs(Z), lower = FALSE)*2 Results from pscl. # pscl results library(pscl) fit_pscl = zeroinfl(count ~ persons + livebait, data = fish, dist = "negbin") Compare. coef B se Z p pscl X.Intercept. -2.803 0.558 -5.026 0.000 pscl persons 0.849 0.124 6.833 0.000 pscl livebait 1.791 0.511 3.504 0.000 pscl Log.theta. -0.969 0.302 -3.206 0.001 pscl X.Intercept..1 -4.276 4.278 -1.000 0.318 pscl persons.1 0.560 0.517 1.084 0.279 pscl livebait.1 1.168 3.661 0.319 0.750 zinb_ll negbin.(Intercept) -2.803 0.558 -5.026 0.000 zinb_ll negbin.persons 0.849 0.124 6.834 0.000 zinb_ll negbin.livebait 1.791 0.511 3.504 0.000 zinb_ll logit.(Intercept) -4.276 4.278 -1.000 0.318 zinb_ll logit.persons 0.560 0.517 1.084 0.279 zinb_ll logit.livebait 1.168 3.661 0.319 0.750 zinb_ll theta -0.969 0.302 -3.206 0.001 ## Supplemental Example This supplemental example uses the bioChemists data. It contains a sample of 915 biochemistry graduate students with the following information: • art: count of articles produced during last 3 years of Ph.D. • fem: factor indicating gender of student, with levels Men and Women • mar: factor indicating marital status of student, with levels Single and Married • kid5: number of children aged 5 or younger • phd: prestige of Ph.D. department • ment: count of articles produced by Ph.D. mentor during last 3 years data("bioChemists", package = "pscl") initial_model = model.matrix(art ~ fem + mar + kid5 + phd + ment, data = bioChemists) # to get X matrix startlogi = glm(art==0 ~ fem + mar + kid5 + phd + ment, data = bioChemists, family = "binomial") startcount = glm(art ~ fem + mar + kid5 + phd + ment, data = bioChemists, family = "quasipoisson") starts = c( negbin = coef(startcount), logit = coef(startlogi), theta = summary(startcount)$dispersion
)

# starts = c(negbin = rep(0, 6),
#            logit = rep(0, 6),
#            theta = 1)

fit_nb_pub = optim(
par = starts ,
fn  = zinb_ll,
X   = initial_model,
y   = bioChemists$art, method = "BFGS", control = list(maxit = 5000, reltol = 1e-12), hessian = TRUE ) # fit_nb_pub B = fit_nb_pub$par
se = sqrt(diag(solve((fit_nb_pub$hessian)))) Z = B/se p = pnorm(abs(Z), lower = FALSE)*2 library(pscl) fit_pscl = zeroinfl(art ~ . | ., data = bioChemists, dist = "negbin") summary(fit_pscl)$coefficients
$count Estimate Std. Error z value Pr(>|z|) (Intercept) 0.4167465901 0.143596450 2.90220678 3.705439e-03 femWomen -0.1955076374 0.075592558 -2.58633447 9.700275e-03 marMarried 0.0975826042 0.084451953 1.15548073 2.478936e-01 kid5 -0.1517320709 0.054206071 -2.79917119 5.123397e-03 phd -0.0006997593 0.036269674 -0.01929323 9.846072e-01 ment 0.0247861500 0.003492672 7.09661548 1.278491e-12 Log(theta) 0.9763577454 0.135469554 7.20721163 5.710921e-13$zero
Estimate Std. Error    z value    Pr(>|z|)
(Intercept) -0.19160645  1.3227962 -0.1448496 0.884829645
femWomen     0.63587048  0.8488959  0.7490559 0.453823498
marMarried  -1.49943716  0.9386562 -1.5974296 0.110169987
kid5         0.62840922  0.4427746  1.4192531 0.155825245
phd         -0.03773288  0.3080059 -0.1225070 0.902497523
ment        -0.88227364  0.3162186 -2.7900755 0.005269575
round(data.frame(B, se, Z, p), 4)
                         B     se       Z      p
negbin.(Intercept)  0.4167 0.1436  2.9021 0.0037
negbin.femWomen    -0.1955 0.0756 -2.5863 0.0097
negbin.marMarried   0.0976 0.0845  1.1555 0.2479
negbin.kid5        -0.1517 0.0542 -2.7992 0.0051
negbin.phd         -0.0007 0.0363 -0.0195 0.9845
negbin.ment         0.0248 0.0035  7.0967 0.0000
logit.(Intercept)  -0.1916 1.3229 -0.1449 0.8848
logit.femWomen      0.6359 0.8490  0.7491 0.4538
logit.marMarried   -1.4995 0.9387 -1.5974 0.1102
logit.kid5          0.6284 0.4428  1.4192 0.1558
logit.phd          -0.0377 0.3080 -0.1225 0.9025
logit.ment         -0.8823 0.3162 -2.7901 0.0053
theta               0.9763 0.1355  7.2071 0.0000