prediction intervals for GLMMs

[bbolker]

I may be completely missing something, but I'm puzzled at the way that predictInterval handles prediction intervals for GLMMs (let's assume logit/probit binomial models, include.resid.var is TRUE [which is what I would normally mean by "prediction intervals" in contrast to "confidence intervals"] and which = "fixed" for now). As far as I can tell, predictInterval sets sigmahat to 1 (this makes some sense since the dispersion parameter of a binomial is fixed to 1), then adds a Gaussian deviate with that standard deviation to the simulated predicted values (which already include the uncertainty due to sampling variance of the fixed-effect parameters).

I can't understand how that makes sense.

• If we wanted to generate prediction intervals based on a Gaussian approximation, we should be setting the standard deviation to something like sqrt(N*p*(1-p)), not 1
• If we wanted to generate prediction intervals based on the assumed conditional distribution, we should be somehow using quantiles of the conditional distribution ...

I apologize if I'm completely misunderstanding/missing something. If I'm not, is there any chance you could point me to a reference?

[jknowles]

Edited to be more thankful!

Thanks @bbolker for this. I really appreciate you taking the time to a) review predictInterval() and b) write up such a detailed issue. I'm really grateful.

I don't think you are misunderstanding or missing anything. The addition of the GLMM specific code was not something I remember looking very deeply into - at the time I built predictInterval for the LMM case for some large models that I couldn't feasibly fit with Stan or approximate with bootMer().

This is not my area of expertise and I didn't document the decisions - my guess is the choice was done to approximate bootMer() output on a very small test set of models, and without understanding of how to assign uncertainty to each portion of the prediction.

I don't use predictInteval() myself very often anymore so this probably explains in part how this has remained in the code this way for a long time. I'm open to suggestions on how to remedy this and improve the interval. I think I could probably implement your first suggestion on my own at some point (time permitting), but I'm not sure I know how I would implement the second at the moment.

I'm open to any guidance you have or your opinion on which approach makes the most sense to focus on implementing.

[bbolker]

In the short term I think I'd disable this feature (stop("prediction intervals are only available for LMMs")), since I don't see how it can make any sense. Since you're already using a simulation-based approach, something like the following would make sense. Generalizing this is a little tricky since you have to have a lookup for which r* function corresponds to each family (the $simulate() method of the family isn't flexible enough). simulate.merMod has a newparams argument that will let you short-circuit this: simulate.glmmTMB could, without a huge amount of extra trouble, but we haven't implemented that yet ...

library(lme4)
library(glmmTMB) ## for simulate_new
library(MASS)    ## for mvrnorm

## simulate example (like ?glmer.nb example, but better)
set.seed(101)
dd <- expand.grid(f1 = factor(1:3),
                  f2 = LETTERS[1:2], g=1:9, rep=1:15,
                  KEEP.OUT.ATTRS=FALSE)
dd$y <-simulate_new( ~ f1 * f2 + (1|g),
                    family = nbinom2,
                    newdata = dd,
                    newparams = list(theta = 1,
                                     beta = rep(1, 6),
                                     betadisp = 1))[[1]]

## fit model
m <- glmmTMB(y ~ f1 * f2 + (1|g), data = dd, family = nbinom2)

## sim fixed effect parameters conditional on theta etc.
sim_parms <- mvrnorm(1000, fixef(m)$cond, vcov(m)$cond)

## data frame for predictions
pdata <- expand.grid(f1 = levels(dd$f1),
                     f2 = levels(dd$f2))

## construct simulated linear predictor
X <- model.matrix(~ f1 * f2, data = pdata)
eta <- t(X %*% t(sim_parms))
mu <- exp(eta)  ## inverse-link function
## simulate observed values and reshape to matrix
R <- rnbinom(length(eta), mu = exp(eta), size = sigma(m))
dim(R) <- dim(eta)
## prediction intervals via quantile
intervals <- t(apply(R, MARGIN = 2, quantile, c(0.025, 0.975)))
print(intervals)
##      2.5%   97.5%
## [1,]    0  26.025
## [2,]    0  93.300
## [3,]    0  70.000
## [4,]    0  76.025
## [5,]    2 581.875
## [6,]    2 653.700

## non-integer values due to quantile interpolation