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)")
```

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()
```

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

`sd(b$t)`

`## [1] 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.24986
##
## $se.boot
## 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)
```

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