R: comparing jglmm models

Question

When using lmer, one can compare two models using anova:

mod.ld.frq = lmer(ln_rt_offset ~ 1 + lfrq + 
                  (1 + lfrq | PID) + (1| Item), 
                  data=ldfw.std)
mod.ld.frq.dur = lmer(ln_rt_offset ~ 1 + lfrq + duration + 
                  (1 + lfrq | PID) + (1| Item), 
                  data=ldfw.std)
anova(mod.ld.frq, mod.ld.frq.dur)

What is the equivalent for jglmm, e.g.:

jmod.ld.frq = jglmm(ln_rt_offset ~ 1 + lfrq + 
                  (1 + lfrq | PID) + (1| Item), 
                  data=ldfw.std)
jmod.ld.frq.dur = jglmm(ln_rt_offset ~ 1 + lfrq + duration + 
                  (1 + lfrq | PID) + (1| Item), 
                  data=ldfw.std)

If I try using anova, I get an error like this:

anova(jmod.ld.frq, jmod.ld.frq.dur.phon)
Error in UseMethod("anova") : 
  no applicable method for 'anova' applied to an object of class "jglmm"

Thank you for any advice!

Answer 1

This doesn't seem to exist, but I wrote a crude anova function based on the existing jglmm::extractAIC.jglmm method, which has to extract the same information (df, log-likelihood) from the model.

library(JuliaCall)
my_anova <- function(m1, m2) {
   julia_assign("model1", m1$model)
   julia_assign("model2", m2$model)
   df1 <- julia_eval("dof(model1);")
   df2 <- julia_eval("dof(model2);")
   loglik1 <- julia_eval("loglikelihood(model1);")
   loglik2 <- julia_eval("loglikelihood(model2);")
   ddf <- df1 - df2
   ddev <- 2*(loglik1 - loglik2)
   c(ddev, ddf, pchisq(ddev, ddf, lower.tail = FALSE))
}

jglmm results:

m1 <- jglmm(Reaction ~ Days + (Days|Subject),lme4::sleepstudy)
m2 <- jglmm(Reaction ~ 1 + (Days|Subject),lme4::sleepstudy)
my_anova(m1, m2)
## [1] 2.353654e+01 1.000000e+00 1.225640e-06

Compare with lme4:

library(lme4)
m1B <- lmer(Reaction ~ Days + (Days|Subject),lme4::sleepstudy)
m2B <- lmer(Reaction ~ 1 + (Days|Subject),lme4::sleepstudy)
anova(m1B, m2B)

Results:

refitting model(s) with ML (instead of REML)
Data: lme4::sleepstudy
Models:
m2B: Reaction ~ 1 + (Days | Subject)
m1B: Reaction ~ Days + (Days | Subject)
    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
m2B    5 1785.5 1801.4 -887.74   1775.5                         
m1B    6 1763.9 1783.1 -875.97   1751.9 23.537  1  1.226e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The lme4 results are much prettier, but the key components (change in deviance, change in df, p-value) are the same.

On second thought, we could have done this more easily by calling the logLik accessor method for the fit (the returned information includes both the log-likelihood and the df) and putting that into the my_anova() function ...