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(Reaction ~ Days + (1 + Days | Subject), data = sleepstudy)
lmer_model
library(mixedup)
find_typical(lmer_model) # closest RE to zero
38;5;246m# A tibble: 2 × 7
[39m
[.5 upper_97.5
group_var effect group value se lower_238;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;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
[.5 upper_97.5 probs
group_var effect group value se lower_238;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
[3m
[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(
lme_1 ~ scale(age) + Sex,
distance data = Orthodont,
random = ~ 1 | Subject,
weights = varIdent(form = ~ 1 | Sex)
)
extract_het_var(lme_1)
38;5;246m# A tibble: 2 × 2
[39m
[
group variance38;5;246m<chr>
[39m
[23m
[3m
[38;5;246m<dbl>
[39m
[23m
[3m
[38;5;250m1
[39m Male 3.10
[38;5;250m2
[39m Female 0.64
[
May 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 V51.000 0.263 -0.460 -0.222 -0.092
V1 0.263 1.000 0.272 0.356 0.081
V2 -0.460 0.272 1.000 0.377 -0.417
V3 -0.222 0.356 0.377 1.000 0.572
V4 -0.092 0.081 -0.417 0.572 1.000
V5
extract_cor_structure(lme_corAR)
38;5;246m# A tibble: 1 × 1
[39m
[
Phi38;5;246m<dbl>
[39m
[23m
[3m
[38;5;250m1
[39m 0.662
[
Just a wrapper for model.frame
.
head(extract_model_data(lmer_model))
38;5;246m# A tibble: 6 × 3
[39m
[
Reaction Days Subject38;5;246m<dbl>
[39m
[23m
[3m
[38;5;246m<dbl>
[39m
[23m
[3m
[38;5;246m<fct>
[39m
[23m
[3m
[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
[