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.
