对 lmer 和 glmer 模型使用 aictab()答案

【问题标题】：using aictab() for lmer and glmer models对 lmer 和 glmer 模型使用 aictab()
【发布时间】：2019-03-21 23:57:55
【问题描述】：

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

我的代码几天前确实有效，但最近发生故障并返回此错误代码

aictab.default() 中的错误：尚未为此对象类定义函数

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)

【问题讨论】：

标签： r lme4

【解决方案1】：

tl;dr你加载了lmerTest 包，所以你的模型有一个不同的类，这很容易混淆aictab()。您可以确保在适合模型时加载了 lme4 和 not 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)

aictab.default 中的错误（cand.set，modnames = 名称，second.ord = TRUE，：尚未为此对象类定义函数

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

【讨论】：

非常感谢您的回复。分离 lmerTest 包效果很好。