# lmer() parametric bootstrap testing for fixed effects

I am performing a parametric bootstrap to test whether I need a specific fixed effect in my model or not. I have mainly done this for exercise and I am interested if my procedure so far is correct.

First, I fit the two models to be compared. One of them includes the effect to be tested for and the other one does not. As I am testing for fixed effects I set REML=FALSE:

    mod8 <- lmer(log(value)
~ matching
+ (sentence_type | subject)
+ (sentence_type | context),
data = wort12_lmm2_melt,
REML = FALSE)
mod_min <- lmer(log(value)
~ 1
+ (sentence_type | subject)
+ (sentence_type | context),
data = wort12_lmm2_melt,
REML = FALSE)


Both models are fit on a balanced data set which includes few missing values. There are slightly above 11000 observations for 70 subjects. Every subject saw every item only one time. The dependent variable are reading times; sentence_type and matching are two-level factors. Context and subject are treated as random effects. Context has 40 levels.

I call anova():

    anova(mod_min, mod8)


and get the output:

    Data: wort12_lmm2_melt
Models:
mod_min: log(value) ~ 1 + (sentence_type |  subject) + (sentence_type | context)
mod8: log(value) ~ matching + (sentence_type |  subject) + (sentence_type |
mod8:     context)
Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod_min  8 3317.6 3375.8 -1650.8   3301.6
mod8     9 3310.9 3376.4 -1646.4   3292.9 8.6849      1   0.003209 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Mistrusting the almighty p I set up a parametric bootstrap by hand:

    mod <- mod8
modnull <- mod_min
lrt.obs <- anova(mod, modnull)$Chisq[2] # save the observed likelihood ratio test statistic n.sim <- 10000 lrt.sim <- numeric(n.sim) dattemp <- mod@frame # pb <- txtProgressBar(min = 0, max = n.sim, style = 3) # set up progress bar to satisfy need for control for(i in 1:n.sim) { # Sys.sleep(0.1) # progress bar related stuff ysim <- unlist(simulate(modnull) # simulate new observations from the null-model modnullsim <- lmer(ysim ~ 1 + (sentence_type | subject) + (sentence_type | context), data = dattemp, REML = FALSE) # fit the null-model modaltsim <- lmer(ysim ~ matching + (sentence_type | subject) + (sentence_type | context), data = dattemp, REML = FALSE) # fit the alternative model lrt.sim[i] <- anova(modnullsim, modaltsim)$Chisq[2] # save the likelihood ratio test statistic
# setTxtProgressBar(pb, i)
}

# assuming chi-squared distribution for comparison

pchisq((2*(logLik(mod8)-logLik(mod_min))),
df    = 1,
lower = FALSE)


with the output:

    'log Lik.' 0.003208543 (df=9)


compare to parametric bootstrap p-value

    p_mod8_mod_min <- (sum(lrt.sim>=lrt.obs)+1)/(n.sim+1)  # p-value. alternative: mean(lrt.sim>lrt.obs)


with the output:

    [1] 0.00319968


Plot the whole thing:

    xx <- seq(0, 20, 0.1)
hist(lrt.sim,
xlim     = c(0, max(c(lrt.sim, lrt.obs))),
col      = "blue",
xlab     = "likelihood ratio test statistic",
ylab     = "density",
cex.lab  = 1.5,
cex.axis = 1.2,
freq     = FALSE)
abline(v   = lrt.obs,
col = "orange",
lwd = 3)
lines(density(lrt.sim),
col = "blue")
lines(xx,
dchisq(xx, df = 1),
col = "red")
box()


which yields:

I do have some questions though:

(1) Is the procedure correct or did I make a mistake?

(2) How is the histogram to be interpreted?

(3) Is the form of the histogram normal or extreme?

Thanks for any help!

• This all looks fine. The histogram is the null distribution of differences in deviance between the full and reduced model. Because you have a large (40) number of levels in your smallest random effect, the likelihood ratio test is accurate -- the p-values based on parametric bootstrapping and on the LRT match almost exactly. You can also use PBmodcomp from the pbkrtest package to run these sorts of comparisons, or KRmodcomp (same package) to get a better (than the LRT) approximation of the p-value. Apr 14 '14 at 10:55

This all looks fine. The histogram is the null distribution of differences in deviance between the full and reduced model. Because you have a large (40) number of levels in your smallest random effect, the likelihood ratio test is accurate -- the p-values based on parametric bootstrapping and on the LRT match almost exactly. You can also use PBmodcomp from the pbkrtest package to run these sorts of comparisons, or KRmodcomp (same package) to get a better (than the LRT) approximation of the p-value.