использование aictab () для моделей lmer и glmer - PullRequest
Новая
       15

использование aictab () для моделей lmer и glmer

1 голос
/ 22 марта 2019

Я пытаюсь сделать выбор модели AICc, используя вывод aictab(). Эта публикация похожа ( Функция не определена при вызове 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)

1 Ответ

2 голосов
/ 22 марта 2019

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 и восстановите модели (мы могли бы просто сделать, например, 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 (канд.set, modnames = names, second.ord = TRUE,: функция, еще не определенная для этого класса объектов

Функция bbmle::AICctab() немного более устойчива, поскольку она опирается только на logLikметод, определенный для класса (по умолчанию он дает таблицу только с дельта-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
...