In many cases, predictions are made by holding the random effects to
zero, which can be seen as the typical case. When the clusters carry
substantive meaning, for example, school, country, or hospital, it might
be of interest to note which are most ‘typical’, or at any other value
relative to their peers. This will work for any model to which
extract_random_effects applies.
library(lme4)
lmer_model <- lmer(Reaction ~ Days + (1 + Days | Subject), data = sleepstudy)
library(mixedup)
find_typical(lmer_model) # closest RE to zero
[38;5;246m# A tibble: 2 × 7
[39m
group_var effect group value se lower_2.5 upper_97.5
[3m
[38;5;246m<chr>
[39m
[23m
[3m
[38;5;246m<chr>
[39m
[23m
[3m
[38;5;246m<fct>
[39m
[23m
[3m
[38;5;246m<dbl>
[39m
[23m
[3m
[38;5;246m<dbl>
[39m
[23m
[3m
[38;5;246m<dbl>
[39m
[23m
[3m
[38;5;246m<dbl>
[39m
[23m
[38;5;250m1
[39m Subject Days 333 -
[31m0
[39m
[31m.
[39m
[31m224
[39m 2.31 -
[31m4
[39m
[31m.
[39m
[31m74
[39m 4.29
[38;5;250m2
[39m Subject Intercept 335 -
[31m0
[39m
[31m.
[39m
[31m334
[39m 12.1 -
[31m24
[39m
[31m.
[39m
[31m0
[39m 23.3
find_typical(lmer_model, probs = c(.1, .9))
[38;5;246m# A tibble: 4 × 8
[39m
group_var effect group value se lower_2.5 upper_97.5 probs
[3m
[38;5;246m<chr>
[39m
[23m
[3m
[38;5;246m<chr>
[39m
[23m
[3m
[38;5;246m<fct>
[39m
[23m
[3m
[38;5;246m<dbl>
[39m
[23m
[3m
[38;5;246m<dbl>
[39m
[23m
[3m
[38;5;246m<dbl>
[39m
[23m
[3m
[38;5;246m<dbl>
[39m
[23m
[3m
[38;5;246m<chr>
[39m
[23m
[38;5;250m1
[39m Subject Days 310 -
[31m5
[39m
[31m.
[39m
[31m45
[39m 2.31 -
[31m9
[39m
[31m.
[39m
[31m97
[39m -
[31m0
[39m
[31m.
[39m
[31m931
[39m 10%
[38;5;250m2
[39m Subject Days 350 6.61 2.31 2.10 11.1 90%
[38;5;250m3
[39m Subject Intercept 370 -
[31m25
[39m
[31m.
[39m
[31m6
[39m 12.1 -
[31m49
[39m
[31m.
[39m
[31m3
[39m -
[31m1
[39m
[31m.
[39m
[31m96
[39m 10%
[38;5;250m4
[39m Subject Intercept 331 22.3 12.1 -
[31m1
[39m
[31m.
[39m
[31m40
[39m 45.9 90% For nlme, only models with varIdent have
been tested. For glmmTMB, this
applies to models that use the diag
function in the formula.
library(nlme)
lme_1 <- lme(
distance ~ scale(age) + Sex,
data = Orthodont,
random = ~ 1 | Subject,
weights = varIdent(form = ~ 1 | Sex)
)
extract_het_var(lme_1)
[38;5;246m# A tibble: 2 × 2
[39m
group variance
[3m
[38;5;246m<chr>
[39m
[23m
[3m
[38;5;246m<dbl>
[39m
[23m
[38;5;250m1
[39m Male 3.10
[38;5;250m2
[39m Female 0.64May include brms at some point,
but for these, you can extract any parameters associated with sigma
model for distributional models using the component
argument of the other standard functions.
With some packages one can estimate the residual correlation structure. Right now this works for nlme, glmmTMB, and brms objects.
base_model <-
lme(Reaction ~ Days,
random = ~ 1 + Days | Subject,
data = sleepstudy)
lme_corSymm <-
update(
base_model,
corr = corSymm(form = ~ 1 | Subject),
data = dplyr::filter(sleepstudy, Days < 5)
)
lme_corAR <-
update(
base_model,
corr = corAR1(form = ~ Days),
data = dplyr::filter(sleepstudy, Days < 5)
)
extract_cor_structure(lme_corSymm)
V1 V2 V3 V4 V5
V1 1.000 0.263 -0.460 -0.222 -0.092
V2 0.263 1.000 0.272 0.356 0.081
V3 -0.460 0.272 1.000 0.377 -0.417
V4 -0.222 0.356 0.377 1.000 0.572
V5 -0.092 0.081 -0.417 0.572 1.000
extract_cor_structure(lme_corAR)
[38;5;246m# A tibble: 1 × 1
[39m
Phi
[3m
[38;5;246m<dbl>
[39m
[23m
[38;5;250m1
[39m 0.662Just a wrapper for model.frame.
head(extract_model_data(lmer_model))
[38;5;246m# A tibble: 6 × 3
[39m
Reaction Days Subject
[3m
[38;5;246m<dbl>
[39m
[23m
[3m
[38;5;246m<dbl>
[39m
[23m
[3m
[38;5;246m<fct>
[39m
[23m
[38;5;250m1
[39m 250. 0 308
[38;5;250m2
[39m 259. 1 308
[38;5;250m3
[39m 251. 2 308
[38;5;250m4
[39m 321. 3 308
[38;5;250m5
[39m 357. 4 308
[38;5;250m6
[39m 415. 5 308