Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

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

share|improve this question
3  
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

1 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 :)

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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