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

Question

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.)

Answer 1

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
## note this requires a fairly recent version of R
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 :)