bootstrapping for lmer with interaction term

Question

I am running a mixed model using lme4 in R:

full_mod3=lmer(logcptplus1 ~ logdepth*logcobb + (1|fyear) + (1 |flocation),
data=cpt, REML=TRUE)

summary:

Formula: logcptplus1 ~ logdepth * logcobb + (1 | fyear) + (1 | flocation)
   Data: cpt

REML criterion at convergence: 577.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.7797 -0.5431  0.0248  0.6562  2.1733 

 Random effects:
 Groups    Name        Variance Std.Dev.
 fyear     (Intercept) 0.2254   0.4748  
 flocation (Intercept) 0.1557   0.3946  
 Residual              0.9663   0.9830  
Number of obs: 193, groups:  fyear, 16; flocation, 16

Fixed effects:
                  Estimate Std. Error t value
(Intercept)        4.3949     1.2319   3.568
logdepth           0.2681     0.4293   0.625
logcobb           -0.7189     0.5955  -1.207
logdepth:logcobb   0.3791     0.2071   1.831

I have used the effects package and function in R to calculate the 95% confidence intervals for the model output. I have calculated and extracted the 95% CI and standard error using the effects package so that I can examine the relationship between the predictor variable of importance and the response variable by holding the secondary predictor variable (logdepth) constant at the median (2.5) in the data set:

gm=4.3949 + 0.2681*depth_median + -0.7189*logcobb_range + 0.3791*
(depth_median*logcobb_range)

ef2=effect("logdepth*logcobb",full_mod3,
     xlevels=list(logcobb=seq(log(0.03268),log(0.37980),,200)))

I have attempted to bootstrap the 95% CIs using code from here. However, I need to calculate the 95% CIs for only the median depth (2.5). Is there a way to specify in the confint() code so that I can calculate the CIs needed to visualize the bootstrapped results as in the plot above?

confint(full_mod3,method="boot",nsim=200,boot.type="perc")

Answer 1

You can do this by specifying a custom function:

library(lme4)
?confint.merMod

FUN: bootstrap function; if ‘NULL’, an internal function that returns the fixed-effect parameters as well as the random-effect parameters on the standard deviation/correlation scale will be used. See ‘bootMer’ for details.

So FUN can be a prediction function (?predict.merMod) that uses a newdata argument that varies and fixes appropriate predictor variables.

An example with built-in data (not quite as interesting as yours since there's a single continuous predictor variable, but I think it should illustrate the approach clearly enough):

fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
pframe <- data.frame(Days=seq(0,20,by=0.5))
## predicted values at population level (re.form=NA)
pfun <- function(fit) {
    predict(fit,newdata=pframe,re.form=NA)
}
set.seed(101)
cc <- confint(fm1,method="boot",FUN=pfun)

Picture:

par(las=1,bty="l")
matplot(pframe$Days,cc,lty=2,col=1,type="l",
     xlab="Days",ylab="Reaction")