Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Pass a named list of models to anova.merMod

I want to be able to pass a named list of models (merMod objects) to anova() and preserve the model names in the output. This is particularly useful in the context of using mclapply() to run a batch of slow models like glmers more efficiently in parallel. The best I've come up with is to use do.call on a de-named version of the model list, but that's not ideal because I might have models named (say) "mod12", "mod17", and "mod16" and these model names get translated to "MODEL1", "MODEL2", and "MODEL3" in the output. (That might seem trivial when looking at a single batch, but over the course of a long modeling session with dozens of models it's a surefire recipe for confusion.)

Note that this isn't the same issue as Create and Call Linear Models from List because I'm not trying to compare pairs of models across lists. It's also more complex than Using lapply on a list of models because I'm using anova() in a non-unary way.

Here's a minimal reprex:

library(lme4)

formList <- c(mod12 = angle ~ recipe + temp + (1|recipe:replicate),
              mod17 = angle ~ recipe + temperature + (1|recipe:replicate),
              mod16 = angle ~ recipe * temperature + (1|recipe:replicate))
modList <- lapply(formList, FUN=lmer, data=cake)

# Fails because modList is named so it's interpreted as arg-name:arg pairs
do.call(anova, modList)

# Suboptimal because model names aren't preserved
do.call(anova, unname(modList))

# Fails because object isn't merMod (or some other class covered by methods("anova"))
do.call(anova, list(object=modList[1], ...=modList[-1], model.names=names(modList)))

The second do.call returns this:

Data: ..1
Models:
MODEL1: angle ~ recipe + temp + (1 | recipe:replicate)
MODEL2: angle ~ recipe + temperature + (1 | recipe:replicate)
MODEL3: angle ~ recipe * temperature + (1 | recipe:replicate)
       Df    AIC    BIC  logLik deviance   Chisq Chi Df Pr(>Chisq)
MODEL1  6 1708.2 1729.8 -848.08   1696.2                          
MODEL2 10 1709.6 1745.6 -844.79   1689.6  6.5755      4     0.1601
MODEL3 20 1719.0 1791.0 -839.53   1679.0 10.5304     10     0.3953

Ideally, the output would look like this:

Data: ..1
Models:
mod12: angle ~ recipe + temp + (1 | recipe:replicate)
mod17: angle ~ recipe + temperature + (1 | recipe:replicate)
mod16: angle ~ recipe * temperature + (1 | recipe:replicate)
      Df    AIC    BIC  logLik deviance   Chisq Chi Df Pr(>Chisq)
mod12  6 1708.2 1729.8 -848.08   1696.2                          
mod17 10 1709.6 1745.6 -844.79   1689.6  6.5755      4     0.1601
mod16 20 1719.0 1791.0 -839.53   1679.0 10.5304     10     0.3953

How do I do this? I'm more than happy with an ugly wrapper around anova() if it means I get a more intelligible output.

like image 952
Dan Villarreal Avatar asked Oct 18 '18 07:10

Dan Villarreal


1 Answers

You can pass a list of symbols like this:

with(modList,
     do.call(anova, 
             lapply(names(modList), as.name)))
#refitting model(s) with ML (instead of REML)
#Data: ..1
#Models:
#mod12: angle ~ recipe + temp + (1 | recipe:replicate)
#mod17: angle ~ recipe + temperature + (1 | recipe:replicate)
#mod16: angle ~ recipe * temperature + (1 | recipe:replicate)
#      Df    AIC    BIC  logLik deviance   Chisq Chi Df Pr(>Chisq)
#mod12  6 1708.2 1729.8 -848.08   1696.2                          
#mod17 10 1709.6 1745.6 -844.79   1689.6  6.5755      4     0.1601
#mod16 20 1719.0 1791.0 -839.53   1679.0 10.5304     10     0.3953
like image 128
Roland Avatar answered Oct 16 '22 00:10

Roland