0
$\begingroup$

I modelled eight lmer models via the lme4 package.

First, I compared the models via a likelihood ratio test yielding this:

npar AIC BIC logLik deviance Chisq Df Pr(>Chisq) baselinemodel 4 63313 63338 -31652 63305 mod1 5 63313 63344 -31651 63303 2.0226 1 0.15497 mod2 6 63232 63270 -31610 63220 82.5963 1 < 2.2e-16 *** mod3 7 63185 63229 -31586 63171 49.0685 1 2.472e-12 *** mod4 8 63186 63237 -31585 63170 0.7713 1 0.37981 mod5 9 63188 63244 -31585 63170 0.3639 1 0.54632 mod6 10 63186 63249 -31583 63166 3.4708 1 0.06246 . mod7 11 63188 63257 -31583 63166 0.3997 1 0.52722

Then, I bootstrapped the lrt values using this kind of code:

set.seed(123)
b1 <- bootMer(baselinemodel,
         FUN = function(x) as.numeric(logLik(x)), 
         nsim = 1000)

set.seed(123)
b2 <- bootMer(mod1, 
          FUN = function(x) as.numeric(logLik(x)), 
          nsim = 1000)


# Observed Likelihood-Ratio
lrt <- as.numeric(-2 * logLik(baselinemodel3) + 2 * logLik(stimM4))

# 1000 bootstrap Likelihood-Ratios
lrt.b <- -2 * b1$t + 2 * b2$t

# 95% CI
qu <- quantile(lrt.b, probs = c(0.025, 0.975))

so far so good. This is what the lrt test gives me, too.

But comparing the other models I get these CIs:

# mod1 vs mod2
#  2.5%     97.5% 
# -169.3204  323.8374 
# -> mod 2 not better than mod1


# mod 2 vs. mod3
#  2.5%    97.5% 
# 46.37453 54.53090 
# -> mod3. better than mod 2


# mod3 vs. mod4
#   2.5%     97.5% 
# 0.7715048 5.9623860 
# -> mod4 better than mod 3

# mod4 vs mod5
#    2.5%     97.5% 
# 0.3611693 5.7582975 
# -> mod5 better than mod4

# mod 5 vs. mod6
# 2.5%    97.5% 
# 3.471502 8.316651 
# -> mod6 better than mod5

# mod 5 vs mod6
#   2.5%     97.5% 
#0.4005908 5.2964914 
# -> mod6 better than mod5

# mod 6 vs mod7
#  2.5%     97.5% 
# -238.5510  255.8117 
# -> mod7 better than mod6

So, according to the lrt test: mod2 > mod1 mod3 > mod2

But, according to the bootstrapped lrt values: mod3> mod2 mod4> mod3 mod5> mod4 mod6> mod5 mod7> mod6

How is this possible? 1. Is it because I got this warning with every bootstrap: "2 warning(s): Model failed to converge with max|grad| = 0.00229957 (tol = 0.002, component 1) (and others)"? 2. Or did I do something wrong?

btw. descriptive data and theory says that mod2> mod1 (which I did not find w/ bootstrapping)

Thanks in advance!!

Edit:

I adapted ping`s code

boot01 <- numeric(100)

for(i in 1:100){
  rmath <- unlist(simulate(baselinemodel))
  bmod <- refit(mod1, rmath)
  smod <- refit(baelinemodel, rmath)
  boot01[i] <- 2*(logLik(bmod)-logLik(smod))
}

pvalue01 <- mean(boot01 > lr01)

# 0.011

boot12 <- numeric(1000)
for(i in 1:1000){
  rmath <- unlist(simulate(mod1))
  bmod <- refit(mod2, rmath)
  smod <- refit(mod1, rmath)
  boot12[i] <- 2*(logLik(bmod)-logLik(smod))
}

pvalue12 <- mean(boot12 > lr12)

# 0

boot23 <- numeric(100)
for(i in 1:100){
  rmath <- unlist(simulate(mod2))
  bmod <- refit(mod3, rmath)
  smod <- refit(mod2, rmath)
  boot23[i] <- 2*(logLik(bmod)-logLik(smod))
}

pvalue23 <- mean(boot23 > lr23)

# 0

According to this code none of my models is significant. Why?

Edit 2:

I fitted the models without REML= FALSE (with ML) and the results corresponded to my lrt-test. But this is still wrong because I am not resampling from the same distribution... or is set.seed(123) accounting for that?

$\endgroup$
  • $\begingroup$ Did you set REML=FALSE? For model comparison you should use ML. $\endgroup$ – user289381 Jul 5 at 9:45
  • $\begingroup$ I used REML= False because this is needed for an lrt-Test. Is it not needed with an lrt test via bootstrapping? $\endgroup$ – a.henrietty Jul 5 at 9:46
  • $\begingroup$ Yes it should be REML=FALSE for all models. $\endgroup$ – user289381 Jul 5 at 9:48
  • $\begingroup$ Your bootstrap procedure is not correct. You are not estimating the logLik with the same resampled data for the various models. $\endgroup$ – user289381 Jul 5 at 10:45
  • $\begingroup$ What is the correct way? I am running the models with ML right now $\endgroup$ – a.henrietty Jul 5 at 10:50
1
$\begingroup$

This is a similar solution to that used here using the pbkrtest R package.

A mistake is that you are combining bootstrapped statistics with different sampled datasets.

If you want to aggregate the bootstrapped log-likelihoods of two models, you need to estimate them with the same resampled datasets.

This is a snippet to perform a parametric bootstrap:

library(lme4)

data(cake)

fm1 <- lmer(angle ~ recipe + (1|recipe:replicate), cake, REML=FALSE)
fm2 <- lmer(angle ~ recipe + temperature + (1|recipe:replicate), cake, REML=FALSE)
fm3 <- lmer(angle ~ recipe * temperature + (1|recipe:replicate), cake, REML=FALSE)

lr12 <- anova(fm1, fm2)$Chisq[2]
lr23 <- anova(fm2, fm3)$Chisq[2]

boot12 <- numeric(1000)
for(i in 1:1000){
  rmath <- unlist(simulate(fm1))
  bmod <- refit(fm2, rmath)
  smod <- refit(fm1, rmath)
  boot12[i] <- 2*(logLik(bmod)-logLik(smod))
}

pvalue12 <- mean(boot12 > lr12)

boot23 <- numeric(1000)
for(i in 1:1000){
  rmath <- unlist(simulate(fm2))
  bmod <- refit(fm3, rmath)
  smod <- refit(fm2, rmath)
  boot23[i] <- 2*(logLik(bmod)-logLik(smod))
}

pvalue23 <- mean(boot23 > lr23)

Comparison of the results:

fm1, fm2

Data: cake
Models:
fm1: angle ~ recipe + (1 | recipe:replicate)
fm2: angle ~ recipe + temperature + (1 | recipe:replicate)
    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
fm1    5 1785.7 1803.7 -887.84   1775.7                         
fm2   10 1709.6 1745.6 -844.79   1689.6 86.106  5  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

fm2, fm3

Data: cake
Models:
fm2: angle ~ recipe + temperature + (1 | recipe:replicate)
fm3: angle ~ recipe * temperature + (1 | recipe:replicate)
    npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
fm2   10 1709.6 1745.6 -844.79   1689.6                    
fm3   20 1719.0 1791.0 -839.53   1679.0 10.53 10     0.3953

and the corresponding results obtained with boostrapping:

mean(boot12 > lr12)
[1] 0

mean(boot23 > lr23)
[1] 0.421
| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you! I adapted my code (see edit above). Also, I fitted without REML=FALSE for my initial attemp and that yielded my results from the lrt test. $\endgroup$ – a.henrietty Jul 5 at 13:51

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.