Bootstrapping for lmer with interaction time
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 used a package and a function effects
in R to compute 95% confidence intervals for model inference. I computed and extracted the 95% CI and standard error using the package effects
so that I can examine the relationship between the severity predictor variable and the response variable by keeping the secondary predictor variable ( logdepth
) constant at the median (2.5) in the dataset:
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 tried to load 95% CI using the code here . However, I need to calculate the 95% CI for only the average depth (2.5). Is there a way to specify in the code confint()
so that I can calculate the CI needed to render the loaded results like in the above graph?
confint(full_mod3,method="boot",nsim=200,boot.type="perc")
source to share
You can do this by specifying a custom function:
library(lme4)
?confint.merMod
FUN: boot function; if "NULL" an internal function will be used that returns the fixed effect parameters as well as the random effect parameters in the standard deviation / correlation scale. See "BootMer" for details.
So there FUN
could be a prediction function ( ?predict.merMod
) that takes an argument newdata
that changes and fixes the appropriate predictor variables.
An example with inline data (not as interesting as yours since there is one continuous predictor variable, but I think this should illustrate the approach well):
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)
Photo:
par(las=1,bty="l")
matplot(pframe$Days,cc,lty=2,col=1,type="l",
xlab="Days",ylab="Reaction")
source to share