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?
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