Coding a Zero Inflated Poisson with regression on p and lambda in r

Question

I am trying to fit a zero-inflated Poisson with regression parameters for lambda as well as p. I am following the framework of:

"Zero-Inflated Poisson Regression with an Application to Defects in Manufacturing, Diane Lambert, Technometrics, Vol 34, No 1 Feb 1992 pp 1-14"

However I cannot seem to get my code to fit right. I thought I had it earlier in the day, but now it is not working. Note that I used this post for some help as well. Link

Here is the code that I am using to simulate ZIP data in r

n<-50 
covL<-seq(0,1,length.out=n)
covp<-seq(0,1,length.out=n)
trueMeans<-exp(1.5-2*covL)
probability<-exp(-1.5+covp*2)/(1+exp(-1.5+covp*2))
U<-runif(n,0,1)
y <- rpois(n,trueMeans)
y[U<probability] <- 0

y is my simulated data. I was working on making sure my coding was correct before putting the code into a function. The following is the code I created:

#initialize values
initial.l<-glm(y[y>0] ~ covL[y>0], family = "poisson")$coefficients
Beta0.l<-as.numeric(initial.l[1])
Beta1.l<-as.numeric(initial.l[2])

Beta0.p<-(sum(y==0)-sum(exp(Beta0.l+Beta1.l*covL)))/length(y)
Beta1.p<-0

#going to use the EM algorithm, so zhat is going to be my latent variable, initialize it here.  

zhat <- rep(0,length(y))
for (i in 1:100) {
#rep through the EM algorithm 100 times, I know this is excessive


#E step
zhat[y==0]<-((1+exp(-(Beta0.p+Beta1.p*covp)-exp(Beta0.l+Beta1.l*covL)))^(-1))[y==0]


  #M step for the coefficients that estimate the mean (lambda)
 lValues<-glm(y ~ covL, family = "poisson",weights=1-zhat)$coefficients
 Beta0.l<-as.numeric(lValues[1])
 Beta1.l<-as.numeric(lValues[2])


#M step for the coefficients that estimate the p value
y_star<-as.numeric(y!=0)
y_star<-c(y_star,y[y==0])
weightsV<-c(1-zhat,zhat[y==0])
covariates<-c(covp,covp[y==0])

 pValues<- glm(y_star~covariates,family=binomial ,weights=weightsV)$coefficients

 Beta0.p<-as.numeric(pValues[1])
 Beta1.p<-as.numeric(pValues[2])

 cat("Iteration: ",i, "  Beta0.l: ",Beta0.l, "  Beta1.l: ",Beta1.l,"  Beta0.p: ", Beta0.p,"  Beta1.p: ", Beta1.p,"\n")
}

The values I should be getting are (values that I am getting are in parentheses):

Beta0.l: 1.5 (1.73)

Beta1.l: -2 (-2.81)

Beta0.p: -1.5 (1.13)

Beta1.p: 2 (-2.62)

What concerns me are the Beta.p values, the signs should be reversed. Any help would be appreciated. Sorry if this is in the wrong place, this is my first time posting. I will try and answer anything that I left out.

Answer 1

While it is possible to use the EM algorithm to estimate ZIP models, my experience is that direct optimization works as reliably in most situations but is somewhat faster (which is not unusual in comparison to EM). The zeroinfl() function in the pscl package implements both approaches (defaulting to EM = FALSE) but it was important to employ analytical gradients and start from reasonable initial values (estimated with glm()).

For your simple data it is sufficient though to use a gradient-free optimizer (below I use Nelder-Mead) and employ dumb starting values (all zero). To make everything exactly reproducible I fix the random seed:

set.seed(1071)
n<-50 
covL<-seq(0,1,length.out=n)
covp<-seq(0,1,length.out=n)
trueMeans<-exp(1.5-2*covL)
probability<-exp(-1.5+covp*2)/(1+exp(-1.5+covp*2))
U<-runif(n,0,1)
y <- rpois(n,trueMeans)
y[U<probability] <- 0

Then, the negative log-likelihood can be computed as (see also Equation 7 in http://dx.doi.org/10.18637/jss.v027.i08):

nll_zip <- function(par) {
  lambda <- exp(par[1] + par[2] * covL)
  p <- exp(par[3] + par[4] * covp) / (1 + exp(par[3] + par[4] * covp))
  lik <- p * (y == 0) + (1 - p) * dpois(y, lambda)
  -sum(log(lik))
}

And Nelder-Mead optimization can be invoked as follows, also computing a numerical Hessian for the standard errors:

opt <- optim(rep(0, 4), nll_zip, hessian = TRUE)

To compare the results to zeroinfl() we can do:

library("pscl")
m <- zeroinfl(y ~ covL | covp)

Then the estimated parameters from both approaches are virtually identical and relatively close to the true parameters we simulated from:

cbind("True" = c(1.5, -2, -1.5, 2), "zeroinfl" = coef(m), "optim" = opt$par)
##                   True  zeroinfl     optim
## count_(Intercept)  1.5  1.350510  1.350757
## count_covL        -2.0 -1.877199 -1.877871
## zero_(Intercept)  -1.5 -1.452235 -1.453011
## zero_covp          2.0  1.739965  1.739817

Also, the asymptotic standard errors essentially agree:

cbind("zeroinfl" = sqrt(diag(vcov(m))), "optim" = sqrt(diag(solve(opt$hessian))))
##                    zeroinfl     optim
## count_(Intercept) 0.2566790 0.2566963
## count_covL        0.7814104 0.7816395
## zero_(Intercept)  0.8410982 0.8414687
## zero_covp         1.8843705 1.8856983