7
$\begingroup$

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:

enter image description here

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!

$\endgroup$

migrated from stackoverflow.com Apr 18 '14 at 21:07

This question came from our site for professional and enthusiast programmers.

  • 5
    $\begingroup$ 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. $\endgroup$ – Ben Bolker Apr 14 '14 at 10:55
  • $\begingroup$ This question appears to be off-topic because it is about interpretation of statistical analysis: belongs on CrossValidated $\endgroup$ – Ben Bolker Apr 14 '14 at 10:57
  • $\begingroup$ Thank you! I flagged the question and asked for it to be migrated to CrossValidated. $\endgroup$ – lord.garbage Apr 14 '14 at 11:14
  • $\begingroup$ @BenBolker I would suggest to convert your comment to a full reply, so that it could be voted up and eventually accepted by the OP. $\endgroup$ – chl Jul 15 '14 at 15:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.