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,
  y   = fish$count,
  method  = "BFGS",
  control = list(maxit = 5000, reltol = 1e-12),
  hessian = TRUE
)
# fit_nb

Comparison

Extract for clean display.

B  = fit_nb$par
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