How to loop through a list of lmerTest outputs to extract deviance component

Question

I am estimating multilevel models using 80 multiply imputed data sets with the package mitml. I can use the testModels() command to compare nested models, but I want to view model fit components (specifically, deviance) for each of the 80 individual imputation models and calculate an overall averaged deviance value.

My model estimates are saved in a mitml.result list called modelt1.

I can extract the deviance value for the first model (of 80) using indexing:

> modelt1[[1]]@devcomp[["cmp"]][["dev"]]
[1] 22637.1

However, I am not sure how to efficiently extract and average all 80 of these values. I know I need to use a loop, but I'm not sure how to combine a loop with indexing like this.

My attempt was something like:

> for(i in modelt1){print(modelt1[[1]]@devcomp[["cmp"]][["dev"]])}
[1] 22637.1

This, unsurprisingly, returns only the deviance for the first model within modelt1.

I tried replacing [[1]] with [[i]], and got an error.

I also attempted to loop through all models like this:

> for(i in modelt1){print(modelt1)}

But this of course provides the full summary output for all 80 models, when I just want the deviance values.

How can I write a loop that will print all 80 deviance values?

Answer 1

You were close. The trick is to use a sequence i in 1:length(fit). Just i in fit yields only a single value and that's the reason why you only get one coefficient.

for (i in 1:length(fit)) print(fit[[i]]@devcomp[["cmp"]][["dev"]])
# [1] 8874.517
# [1] 8874.517
# [1] 8874.517
# [1] 8874.517
# [1] 8874.517

However, since R is a vectorized language, I suggest (in most cases) not using for loops and getting used to sapply & Co. for speed and convenience reasons.

Example:

library(mitml)
fml <- ReadDis + SES ~ ReadAchiev + (1|ID)
imp <- panImpute(studentratings, formula=fml, n.burn=1000, n.iter=100, m=5)
implist <- mitmlComplete(imp, print=1:5)

library(lme4)
fit <- with(implist, lmer(ReadAchiev ~ (1|ID), REML=FALSE))

sapply(seq(fit), function(i) fit[[i]]@devcomp[["cmp"]][["dev"]])
# [1] 8874.517 8874.517 8874.517 8874.517 8874.517