对 lmer 和 glmer 模型使用 aictab()

using aictab() for lmer and glmer models

我正在尝试使用 aictab() 输出进行 AICc 模型选择。这个帖子是相似的() but does not apply to my problem because it used glm.nb models and I am using lmer models and also glmer(family=binomial). The aictab documentation (https://www.rdocumentation.org/packages/AICcmodavg/versions/2.2-1/topics/aictab)说函数可以处理lmlmeglm模型:那包括lmerglmer吗?

我的代码在几天前确实有效,但最近出现故障并返回此错误代码

Error in aictab.default() : Function not yet defined for this object class

m1 <- lmer(y ~ x + (1|z), data=df)
m2 <- lmer(y ~ w + (1|z), data=df)
m3 <- lmer(y ~ v + (1|z), data=df)
cand.set <- list(m1, m2, m3)
names <- list("x", "w", "v")
aictab(cand.set, modnames=names, second.ord=TRUE, nobs=NULL, sort=TRUE)

tl;dr 您加载了 lmerTest 包,因此您的模型具有不同的 class,这令人困惑 aictab()。您可以确保在拟合模型时加载了 lme4 而不是 lmerTest,或者使用 bbmle::AICctab(),这似乎更稳健。 (可能值得就此联系 AICcmodavg 的包维护者......)

设置 lme4:

library(lme4)
library(AICcmodavg)

ss <- transform(sleepstudy, junk1=rnorm(nrow(sleepstudy)),
                junk2=rnorm(nrow(sleepstudy)))
m1 <- lmer(Reaction ~ Days + (1|Subject), data=ss, REML=FALSE)
m2 <- update(m1, . ~ . - Days + junk1)
m3 <- update(m1, . ~ . - Days + junk2)
cand.set <- list(m1, m2, m3)
names <- c("Days", "junk1", "junk2")  ## this should be a vector, not a list ...

这很好用:

aictab(cand.set, modnames=names, second.ord=TRUE, nobs=NULL, sort=TRUE)
##       K    AICc Delta_AICc AICcWt Cum.Wt      LL
## Days  4 1802.31       0.00      1      1 -897.04
## junk2 4 1918.47     116.16      0      1 -955.12
## junk1 4 1918.71     116.40      0      1 -955.24

现在加载 lmerTest 并改装模型(我们可以这样做,例如 m1 <- as(m1, "lmerModLmerTest") 但改装很容易)。

library(lmerTest)
m1 <- lmer(Reaction ~ Days + (1|Subject), data=ss, REML=FALSE)
m2 <- update(m1, . ~ . - Days + junk1)
m3 <- update(m1, . ~ . - Days + junk2)
cand.set <- list(m1, m2, m3)
aictab(cand.set, modnames=names, second.ord=TRUE, nobs=NULL, sort=TRUE)

Error in aictab.default(cand.set, modnames = names, second.ord = TRUE, : Function not yet defined for this object class

bbmle::AICctab() 函数更健壮一点,因为它仅依赖于为 class 定义的 logLik 方法(默认情况下它给出 table只有 delta-AIC 和 df)

library(bbmle)
AICctab(m1, m2, m3, mnames=names, base=TRUE, weights=TRUE, logLik=TRUE)
##       logLik AICc   dLogLik dAICc  df weight
## Days  -897.0 1802.3   58.2     0.0 4  1     
## junk2 -955.1 1918.5    0.1   116.2 4  <0.001
## junk1 -955.2 1918.7    0.0   116.4 4  <0.001