CODEWITHSUNDEEP
rlme4
user7898623
user7898623

Reputation: 43

using aictab() for lmer and glmer models

I am trying to do AICc model selection using aictab() output. This posting is similar (Function not defined when calling aictab) 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) says the function can handle lm, lme and glm models: does that include lmer and glmer?

My code actually worked a few days ago, but recently broke and returned this error code

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)

Upvotes: 4

Views: 5508

Answers (1)

Ben Bolker
Ben Bolker

Reputation: 226961

tl;dr you loaded the lmerTest package, so your models have a different class, which is confusing aictab(). You could either make sure you have lme4 and not lmerTest loaded when you fit the models, or use bbmle::AICctab(), which seems slightly more robust. (It might be worth contacting the package maintainer of AICcmodavg about this ...)

Set up with 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 ...

This works fine:

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

Now load lmerTest and refit the models (we could just do e.g. m1 <- as(m1, "lmerModLmerTest") but it's easy enough to refit).

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

The bbmle::AICctab() function is a little bit more robust because it relies only on the logLik method being defined for a class (by default it gives a table with only the delta-AIC and 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

Upvotes: 4

Related Questions