r - 对 lmer 和 glmer 模型使用 aictab()

Question

我正在尝试使用 aictab() 输出进行 AICc 模型选择。该帖子与 (Function not defined when calling aictab) 类似，但不适用于我的问题，因为它使用了 glm.nb 模型，而我使用的是 lmer 模型以及 glmer(家庭=二项式)。 aictab 文档 ( https://www.rdocumentation.org/packages/AICcmodavg/versions/2.2-1/topics/aictab ) 说该函数可以处理 lm、lme 和 glm 模型:是吗包括 lmer 和 glmer？

我的代码在几天前实际上可以工作，但最近坏了并返回了这个错误代码

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)

Answer 1

tl;dr 您加载了 lmerTest包，因此您的模型具有不同的类，这令人困惑 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并 retrofit 模型(我们可以做例如 m1 <- as(m1, "lmerModLmerTest")，但 retrofit 起来很容易)。

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()函数更健壮一点，因为它只依赖于 logLik为一个类定义的方法(默认情况下它给出一个只有 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