3
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ 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! $\endgroup$
    – jbowman
    Commented Jun 14, 2017 at 19:28
  • $\begingroup$ 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! $\endgroup$
    – mprice
    Commented Jun 14, 2017 at 21:39

1 Answer 1

5
$\begingroup$

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
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – mprice
    Commented Jun 14, 2017 at 21:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.