R:比较jglmm模型

Question

在使用lmer时，可以比较两种使用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)

什么是jglmm的等价物，例如：

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)

如果我尝试使用anova，就会得到这样的错误：

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"

谢谢你的建议！

Answer 1

这似乎不存在，但我根据现有的anova方法编写了一个粗略的jglmm::extractAIC.jglmm函数，该方法必须从模型中提取相同的信息(df，log-似然)。

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结果：

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

与lme4的比较

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

结果：

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

lme4结果要漂亮得多，但关键部件(偏差变化、df变化、p值变化)是相同的。

再想一想，我们可以通过调用logLik访问器方法(返回的信息包括日志似然和df)并将其放入my_anova()函数中，可以更容易地做到这一点.