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
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 (
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,
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.
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
library(visreg) visreg(fm1, "Days", ylab="Reaction time (ms)")
If we are interested in the fitted value at some specified time, let’s say 5 days, we can use
re.form=NA to yield predictions for the average subject (i.e. not include random effects in the predictions). Note that we actually invoke
fm1 is fitted with
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
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()
The standard error of our predicted value can be estimated simply as the standard deviation of the sampling distribution, that is,
##  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
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
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.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)
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)