1
$\begingroup$

I'm working on zero inflated poisson models but I have a doubt on the estimation of the coefficients. Suppose that I have a small sample of data (just for an example, 10 policies with just 4 covariates (maybe sex, age, family composition and HP_car)). I have understood how the zero inflated poisson works and how to interpret the results in R, but I would like to understend all the single passage of the estimation.

It means that: logit(p)=G * gamma and log(lambda)=B*beta

where G are the covariates sex and age for example and B family composition and HP_car.

How can I estimate beta and gamma coefficients BY HAND? The maximum likelihood I think is too complex and for the EM algorithm how can I procede? Have I to set an initial value for each beta and gamma and compute the Z (defined as zero if the dependent variable is >0 and p/(p+e^(-λ)∙(1-p)) if zero)?

Thank you!

$\endgroup$
1
$\begingroup$

The estimation via full maximum likelihood is easiest, I would say, at least if you are willing to use an all-purpose optimizer like optim() in R. Then you just need to set up the (negative) log-likelihood and some starting values. For full EM you additionally need to iterate between weighted count GLM and modified binary GLM.

For illustration I replicate the following (over-)simplified regression result using zeroinfl() from pscl:

library("pscl")
data("bioChemists", package = "pscl")
m <- zeroinfl(art ~ fem | ment, data = bioChemists)
coef(m)
## count_(Intercept)    count_femWomen  zero_(Intercept)         zero_ment 
##         0.8484406        -0.2197297        -0.3632207        -0.1658713 
logLik(m)
## 'log Lik.' -1640.912 (df=4)

To initialize our replication we first fit a Poisson GLM for the count regression model (ignoring that there are excess zeros) and a binary GLM for the zeros (rather than the excess zeros). We also extract the corresponding regressor matrices x and z, respectively, and the response y.

## count initialization
mp <- glm(art ~ fem, data = bioChemists, family = poisson)
x <- model.matrix(mp)
## binary initialization
mb <- glm(factor(art == 0) ~ ment, data = bioChemists, family = binomial)
z <- model.matrix(mb)
## response
y <- bioChemists$art

Then we implement the (negative) log-likelihood of the zero-inflated Poisson model, see Section 2.3 of vignette("countreg", package = "pscl") for the equations and explanations. In the code below nll() assumes a parameter vector of length 4, uses exp() and plogis() for the inverse link functions of log and logit, respectively, and dpois() for the density of the Poisson distribution:

nll <- function(par) {
  lambda <- exp(x %*% par[1:2])
  ziprob <- plogis(z %*% par[3:4])
  zidens <- ziprob * (y == 0) + (1 - ziprob) * dpois(y, lambda = lambda)
  -sum(log(zidens))
}

With the log-likelihood nll() and the starting values from the models mp and mb you can directly call an all-purpose optimizer like optim(). For more challenging data sets, this might have convergence problems and better starting values or analytical gradients would be useful...but here direct optimization already works well enough and replicates the results from above:

optim(c(coef(mp), coef(mb)), nll, method = "BFGS")
## $par
## (Intercept)    femWomen (Intercept)        ment 
##   0.8484319  -0.2197344  -0.3631972  -0.1658844 
## 
## $value
## [1] 1640.912
## 
## $counts
## function gradient 
##       31       10 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

For the EM algorithm you need to take the fitted values from the initial mp and mb models and compute the updated zero-inflation probabilities by the equation you indicated in your question.

lambda <- fitted(mp)
ziprob <- fitted(mb)
ziprob <- ziprob/(ziprob + (1 - ziprob) * dpois(0, lambda))
ziprob[y > 0] <- 0

Then you need the (negative) log-likelihood at the starting values and a suboptimal value so that the iteration below is started:

nll_new <- nll(c(coef(mp), coef(mb)))
nll_old <- 2 * nll_new

Finally, you run a loop of weighted Poisson GLMs and binary GLMs for the zero-inflation probability (throwing a warning in R because the response is not a binary variable).

while(abs((nll_old - nll_new)/nll_old) > 1e-7) {
  nll_old <- nll_new
  mp <- glm(art ~ fem, data = bioChemists, family = poisson, weights = 1 - ziprob)
  mb <- glm(ziprob ~ ment, data = bioChemists, family = binomial)
  lambda <- fitted(mp)
  ziprob <- fitted(mb)
  ziprob <- ziprob/(ziprob + (1 - ziprob) * dpois(0, lambda))
  ziprob[y > 0] <- 0
  nll_new <- nll(c(coef(mp), coef(mb)))
}

This also replicates the results:

c(coef(mp), coef(mb))
## (Intercept)    femWomen (Intercept)        ment 
##   0.8486745  -0.2197011  -0.3641460  -0.1653609 
nll_new
## [1] 1640.912

The zeroinfl() function in pscl (or in the successor package countreg on R-Forge) essentially implements these two strategies for estimating the model. It just offers a lot more options in terms of distributions, link functions, etc.; uses analytical gradients; and returns a nicely classed object with many methods etc.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you very much, your help has been really useful! $\endgroup$ – Castels Jan 10 at 12:55
  • $\begingroup$ You're welcome. Please also accept the answer so that it is marked as resolved here on StackExchange. $\endgroup$ – Achim Zeileis Jan 10 at 16:11
  • $\begingroup$ How can I do that? $\endgroup$ – Castels Jan 11 at 8:53
  • $\begingroup$ See stackoverflow.com/help/someone-answers $\endgroup$ – Achim Zeileis Jan 11 at 9:06
  • $\begingroup$ I have another question about the topic. I have perfectly understood how to estimate the parameters for the Poisson part while is not really clear for me the estimation for the logistic part. Using the EM algorithm, once you have the initial paramaters, for the Poisson part you run a new glm but weighted (1-z).. but for the Logistic part? Lambert makes a long reasoning but at the end she says that gamma can be found again with a logistic regression weighted. But again with (1-z)? Because in the code you simple write (mb <- glm(ziprob ~ ment, data = bioChemists, family = binomial)) $\endgroup$ – Castels Jan 11 at 11:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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