How to get standard errors from R zero-inflated count data regression?

The following code

PredictNew <- predict (glm.fit, newdata = Predict, X1 =X1, Y1= Y1,
type = "response", se.fit = TRUE)


produces a 3-column data.frame--PredictNew, the fitted values, the standard errors and a residual scale term.

Perfect... However using a model fitted with zeroinfl {pscl}:

PredictNew <- predict (zeroinfl.fit, newdata = Predict, X1 =X1, Y1= Y1,
type = "response", se.fit = TRUE)


or

PredictNew <- predict (zeroinfl.fit, newdata = Predict, X1 =X1, Y1= Y1,
type = "response", se.fit = TRUE, MC = 2500, conf = .95))


produce a single column vector of fitted values only. I would, however, be very keen to have standard errors. Everything I have read says that they should be produced..

(Code has been simplified somewhat, I actually have four variables and an offset - no probs with the predict.glm and se.fit = TRUE producing SEs.)

-
Have a look at this thread on R-Help: stat.ethz.ch/pipermail/r-help/2008-December/thread.html#182806 (particularly the message from Achim Zeileis who provides code to do what I think you're trying to do). It doesn't look like standard errors have been implemented into the predict() function for zeroinfl() at this time. –  smillig Jun 26 '12 at 10:15
Thanks , that code seemed to produce pretty reasonable results. Others should note that the predict() parameter in the new zeroinfl.predict function for se.fit= TRUE was changed to se= TRUE, in order to extract the predicted intervals and s.e. –  KalahariKev Jun 28 '12 at 2:11

To my knowledge, the predict method for results from zeroinfl does not include standard errors. If your goal is to construct confidence intervals, one attractive alternative is to use bootstrapping. I say attractive because bootstrapping has the potential to be more robust (at a loss of efficiency if all the assumptions for the SEs are met).

Here is some rough code to do what you want. It will not work exactly, but hopefully you can make the necessary corrections.

## load boot package
require(boot)
## output coefficients from your original model
## these can be used as starting values for your bootstrap model
## to help speed up convergence and the bootstrap
dput(round(coef(zeroinfl.fit, "count"), 3))
dput(round(coef(zeroinfl.fit, "zero"), 3))

## function to pass to the boot function to fit your model
## needs to take data, an index (as the second argument!) and your new data
f <- function(data, i, newdata) {
require(pscl)
m <- zeroinfl(count ~ child + camper | persons, data = data[i, ], start = list(count = c(1.598, -1.0428, 0.834), zero = c(1.297, -0.564)))
mparams <- as.vector(t(do.call(rbind, coef(summary(m)))[, 1:2]))
yhat <- predict(m, newdata, type = "response")
return(c(mparams, yhat))
}

## set the seed and do the bootstrap, make sure to set your number of cpus
set.seed(10)
res <- boot(dat, f, R = 1200, newdata = Predict, parallel = "snow", ncpus = 4)

## get the bootstrapped percentile CIs
## the 10 here is because in my initial example, there were 10 parameters before predicted values
yhat <- t(sapply(10 + (1:nrow(Predict)), function(i) {
out <- boot.ci(res, index = i, type = c("perc"))
with(out, c(Est = t0, pLL = percent[4], pUL = percent[5]))
}))

## merge CIs with predicted values
Predict<- cbind(Predict, yhat)


I drew this code from two pages I wrote, one bootstrapping parameters from a zero-inflated poisson regression with zeroinfl Zero-inflated poisson and one demonstrating how to get bootstrapped confidence intervals for predicted values from a zero-truncated negative binomial model Zero-truncated negative binomial. Combined, hopefully that provides you sufficient examples to get it working with predicted values from a zero-inflated poisson. You may also get some graphing ideas :)

-