lmer() parametric bootstrap testing for fixed effects

Question

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!

Answer 1

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.

Answer 2

While assuming the model is correct regarding all assumptions, then anova function in lmerTest package will give us exact p-values for testing fixed effects. The parametric booktstrap method is mainly used for testing the random effects if you want to obtain a less conservative p-values of exact likelihood ratio tests (LRTs) than those of the asymptotic LRTs given by ranova function in lmerTest.