I've used MuMIn::model.avg() to average a subset of models derived from MuMIn::dredge(). It seems only model averaged fixed effect estimates are held in the model averaged object. But you can get subject-level predictions, so it's using random effect estimates somehow. Does anyone know how to inspect random effect estimates from a MuMIn::model.avg() object, or how to calculate them from the model list? A working example...
library(glmmTMB)
library(MuMIn)
data(iris)
#make more covariates, factor species, center/scale
data <- iris |> dplyr::mutate(Species=as.factor(Species),
Sepal.Length.2=Sepal.Length^2,
Sepal.Width.2=Sepal.Width^2,
Petal.Width.2=Petal.Width^2,
Sepal.Length.Log=log(Sepal.Length),
Sepal.Width.Log=log(Sepal.Width),
Petal.Width.Log=log(Petal.Width)) |>
dplyr::mutate(dplyr::across(c(-Species,-Petal.Length), ~c(scale(.x))))
#full glmm
mod <- glmmTMB::glmmTMB(Petal.Length~Sepal.Length+Sepal.Length.2 + Sepal.Length.Log +
Sepal.Width+Sepal.Width.2+Sepal.Width.Log +
Petal.Width+Petal.Width.2+Petal.Width.Log +
(1|Species), data, na.action='na.fail')
#dredge model subsets with random effect fixed
mmodels <- MuMIn::dredge(mod, m.lim=c(2,3), fixed = ~(1|Species))
#Get subset of top models
confset.95p <- MuMIn::get.models(mmodels, subset = cumsum(weight) <= .95)
#get model average from the subset
avgm <- MuMIn::model.avg(confset.95p)
#coefficients (fixed effects only)
avgm$coefficients
# cond((Int)) cond(Petal.Width) cond(Sepal.Length) cond(Sepal.Length.2) cond(Sepal.Length.Log)
# full 3.758 0.435179 0.07070698 0.7536782 -0.3883600
# subset 3.758 0.435179 0.13420163 0.9513564 -0.7888739
# cond(Sepal.Width.2) cond(Sepal.Width) cond(Sepal.Width.Log)
# full -0.008693878 -0.00390457 -0.001576163
# subset -0.080420866 -0.07202914 -0.059919304
#could maybe use weights and random effect estimates below to get weighted mean?
avgm$msTable$weight
# [1] 0.31908561 0.28451108 0.20778561 0.10810475 0.05420819 0.02630476
#random intercept estimates for each model
lapply(confset.95p,ranef)
# $`26`
# $Species
# (Intercept)
# setosa -1.3814974
# versicolor 0.4405681
# virginica 0.9409315
# ...
#demonstrate predictions are using random effect estimates..but what are they?
test <- data.frame(avgm$x,Species=data$Species)[-1]
predict(avgm, test, type='response', re.form=~0)[1]#population-level predictions not using random effect estimates (i.e., set to zero)
predict(avgm, test, type='response', re.form=NULL)[1]#subject-level prediction using random effect estimates