Confidence intervals on predictions from mixed-effects models

PUBLISHED ON JUN 1, 2017 — R

With simple linear regression, standard errors and confidence intervals for fitted (and predicted) values are easily computed. In R, we can use the se.fit argument in predict.lm, which returns the standard error for the fitted values, and interval = "confidence" to return confidence intervals. With linear mixed-effects models, however, it is not so easy. Neither predict.lme (from nlme) nor predict.merMod (from lme4) provide these methods, as confidence intervals on mixed-effects model predictions are harder to produce.

The solution is to use the parametric bootstrap, which is conveniently implemented in bootMer to be applied to models fit with the lme4 package (lmer, not glmer). Here I describe a simple wrapper around bootMer, providing an alternative for predict.merMod that calculates standard errors (and confidence intervals) for predictions.

A side-effect of this implementation is that confidence intervals now appear in visreg plots of mixed-effects models.

The code shown below is implemented as a very simple package, bootpredictlme4, which can be installed with,

devtools::install_github("remkoduursma/bootpredictlme4")

Note that the package cannot be used in conjunction with the lmerTest package, since both packages replace the predict.merMod function from lme4. For this reason I will probably not attempt to submit this package to CRAN.

Bootstrapped confidence intervals

Suppose we fit a linear mixed-effects model, using the built-in sleepstudy dataset (in the lme4 package), which fits the reaction time of subjects as a function of the number of days since a sleep deprivation (3 hours per night) experiment.

library(lme4)
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)

The structure allows for a random intercept and slope of reaction vs. days by subject. A visreg plot (Fig. 1) shows the fitted model (for the average subject). Note that I have not yet loaded the bootpredictlme4 package.

library(visreg)
visreg(fm1, "Days", ylab="Reaction time (ms)") Figure 1: A simple visreg plot of the sleepstudy mixed-effects model.

If we are interested in the fitted value at some specified time, let’s say 5 days, we can use predict, setting re.form=NA to yield predictions for the average subject (i.e. not include random effects in the predictions). Note that we actually invoke predict.merMod, because fm1 is fitted with lmer.

predict(fm1, newdata = data.frame(Days=5), re.form=NA)
##        1
## 303.7415

Now we would like to calculate a confidence interval on this prediction, using the bootstrap. We can generate 100 bootstrap resamples with the bootMer function, as follows. In this case, we apply the predict.merMod to each of the bootstrapped fits. I am just using 100 fits to speed up this code, but you probably want to use 1000 samples or more.

b <- bootMer(fm1, nsim=100,
FUN=function(x)predict(x, newdata=data.frame(Days=5), re.form=NA))

The resulting object b contains lots of information, but most useful is the element b$t, which lists the output of the function we applied to each of the bootstrap resamples, in this case the fitted value at Days=5. The following histogram shows the distribution of the fitted values, that is, the sampling distribution. Note that if you use more samples (which you should), the distribution will look more like a normal distribution. hist(b$t, breaks=seq(250,350,by=5),
ylim=c(0,25),
main="", xlab="Reaction time at 5 Days (ms)",
col="cornflowerblue")
box() Figure 2: The sampling distribution of the fitted Reaction time at 5 days, a histogram of bootstrap estimates.

The standard error of our predicted value can be estimated simply as the standard deviation of the sampling distribution, that is,

sd(b$t) ##  9.184015 and the confidence interval as the quantiles of the distribution (in this case, a 95% confidence interval). quantile(b$t, probs=c(0.025, 0.975))
##     2.5%    97.5%
## 286.5467 321.5505

A simple R package

The example above is implemented in the bootpredictlme4 package (as mentioned, install with devtools::install_github("remkoduursma/bootpredictlme4")). When loading the package, the predict.merMod function is replaced with a function that takes an se.fit argument.

library(bootpredictlme4)
predict(fm1, newdata=data.frame(Days=5), re.form=NA, se.fit=TRUE, nsim=100)
## $fit ## 1 ## 303.7415 ## ##$se.fit
##        1
## 10.50841
##
## $ci.fit ## 1 ## 2.5% 283.6878 ## 97.5% 323.8666 Two standard errors are computed, the se.boot is the one described above (the standard deviation of the sampling distribution), and se.fit is an effective standard error that reconstructs the bootstrapped confidence interval if we were to assume asymptotic normality. In the case of normality, the half-width of the confidence interval (for 95% coverage) is of course 1.96 times the standard error. Thus se.fit is calculated as the half-width of the bootstrapped confidence interval, divided by 1.96. It is this standard error that is used by visreg (see below), making sure that the resulting confidence interval is actually equal to the bootstrapped confidence interval. In most cases se.fit and se.boot will be very similar, unless you used very few bootstrap samples, or the sampling distribution is very non-normal (which is certainly possible). A major advantage is that visreg automatically recognizes that we can compute standard errors, and now a confidence interval appears on the standard effects plot. Note that bootpredictlme4 warns about using few bootstrap resamples. The default has a low number to allow rapid testing. visreg(fm1, "Days", ylab="Reaction time (ms)") ## Number of bootstrap replicates very low. ## Set to higher value with e.g. options(bootnsim = 500) Figure 3: A simple visreg plot of the sleepstudy mixed-effects model - with the bootpredictlme4 package loaded. The new predict method also allows more complicated visreg plots, like the overlay plot (my personal favorite). # Add a fake grouping variable to sleepstudy high <- with(sleepstudy, levels(reorder(Subject,Reaction,mean)))[1:9] sleepstudy$Group <- factor(ifelse(sleepstudy\$Subject %in% high, "A", "B"))

fm2 <- lmer(Reaction ~ Days*Group + (Days | Subject), sleepstudy)
visreg(fm2, "Days", by="Group", overlay=TRUE, ylab="Reaction time (ms)")
## Number of bootstrap replicates very low.
##  Set to higher value with e.g. options(bootnsim = 500)
## Number of bootstrap replicates very low.
##  Set to higher value with e.g. options(bootnsim = 500) Figure 4: A visreg plot with a factor variable included in the mixed-effects model. Confidence intervals are calculated with the bootstrap via bootpredictlme4.