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

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. • Look carefully at the formula for your E-step, especially the parentheses. Are you sure you got them right? You have, simplified, a term that looks like$1 + \exp(... -\exp(...))$. That caught my eye, but I haven't checked the math to make sure I'm right! Commented Jun 14, 2017 at 19:28 • It does look peculiar, but in the Lambert paper it does have it in this form, I have compared this to another method of calculating zhat, and it seems to be working all right. Thanks for pointing this out though! Commented Jun 14, 2017 at 21:39 ## 1 Answer 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

• Thanks for the info, I am going to play around with this and see how this works. I think the zeroinfl() function is exactly what I need though. I wish I could have found that sooner. Commented Jun 14, 2017 at 21:40