1

Is it possible to get prediction intervals from a model average in R?

I've used the MuMIn package to model-average several linear mixed models (that I fit using lme4::lmer()). The MuMIn package supports model predictions & st. errors of estimates (if all of the component models support the estimation of st. errors), which are convenient for getting an [estimated][1] confidence interval on the prediction.

To get a prediction interval from a single linear mixed model fit using lme4::lmer(), I could follow Ben Bolker's instructions: 

library(lme4)
data("Orthodont",package="MEMSS")
fm1 <- lmer(
    formula = distance ~ age*Sex + (age|Subject)
    , data = Orthodont
)
newdat <- expand.grid(
    age=c(8,10,12,14)
    , Sex=c("Female","Male")
    , distance = 0
)
newdat$distance <- predict(fm1,newdat,re.form=NA)
mm <- model.matrix(terms(fm1),newdat)
## or newdat$distance <- mm %*% fixef(fm1)
pvar1 <- diag(mm %*% tcrossprod(vcov(fm1),mm))
tvar1 <- pvar1+VarCorr(fm1)$Subject[1]  ## must be adapted for more complex models
cmult <- 2 ## could use 1.96
newdat <- data.frame(
    newdat
    , plo = newdat$distance-cmult*sqrt(pvar1) # Confidence Interval
    , phi = newdat$distance+cmult*sqrt(pvar1) # Confidence Interval
    , tlo = newdat$distance-cmult*sqrt(tvar1) # Prediction Interval
    , thi = newdat$distance+cmult*sqrt(tvar1) # Prediction Interval
)

But how could I do this for several models that are averaged together? This gives me a [rough][1] confidence interval, but it's unclear to me how to average the prediction interval across models:

library(lme4)
library(MuMIn)
data("Orthodont",package="MEMSS")
fit_full <- lmer(
   formula = distance ~ age*Sex + (age|Subject), 
   data = Orthodont,
   REML = FALSE, 
   na.action = 'na.fail'
)

fit_dredge <- dredge(fit_full)

fit_ma <- model.avg(object = get.models(fit_dredge, subset = delta <= 4))

newdat <- expand.grid(
   age=c(8,10,12,14),
   Sex=c("Female","Male"),
   distance = 0
)

predicted <- predict(fit_ma,newdat,re.form=NA, se.fit = TRUE)
newdat$distance <- predicted$fit
newdat$distance_lower_CI <- predicted$fit - 1.96*predicted$se.fit
newdat$distance_upper_CI <- predicted$fit + 1.96*predicted$se.fit

[1] As Ben Bolker notes here, these confidence intervals only account for uncertainty in the fixed effects, not uncertainty in the random effects. lme4::bootMer() will give a better estimate of the confidence interval, but it only works on a single model, not a model-average.

filups21
  • 1,611
  • 1
  • 19
  • 22
  • 1
    this question is partly programming (how would I do that) and partly statistical/algorithmic (what is even a sensible *theoretical* way to solve this problem ...) – Ben Bolker Feb 06 '18 at 22:00
  • Well, this is beyond my ken (hence the original question), but I thought one possibility might be to use `lme::bootMer()` to get predictions from each candidate model (accounting for uncertainty in fixed and random effects as much as possible), but with `nsim` proportional to each model's AIC weight. Then combine the predictions from all the models and summarize. Is that at all sensible? – filups21 Feb 08 '18 at 19:26

0 Answers0