How to inspect random effect estimates from MuMIn::model.avg() object?

14 views Asked by At

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
0

There are 0 answers