使用aictab()函数对lmer和glmer模型进行分析。

Question

我正在尝试使用 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是否包括在内呢？

我的代码几天前实际上是可行的，但最近出现了以下错误代码：

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)

- user7898623

Answer 1

您加载了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并重新拟合模型（我们可以只是执行例如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 = names, second.ord = TRUE)函数发生错误：此对象类别尚未定义该函数

由于只依赖于类别的logLik方法被定义（默认情况下仅显示delta-AIC和df），因此bbmle::AICctab()函数更加健壮。

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

- Ben Bolker

---

非常感谢您的回复。分离lmerTest包完美解决了问题。 - user7898623