Zero Inflated Poisson model, estimation of vectors beta and gamma 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!
 A: 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.
