how to extract the random effect in multilevel modeling using lmer in r?

Question

For example, this is the result of certain multilevel analysis

MLM1<-lmer(y ~ 1 + con + ev1 + ev2 + (1 | pid),data=dat_ind)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: y ~ 1 + con + ev1 + ev2 + (1 | pid)
   Data: dat_ind

REML criterion at convergence: 837

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.57771 -0.52765  0.00076  0.54715  2.27597 

Random effects:
 Groups   Name        Variance Std.Dev.
 pid      (Intercept) 1.4119   1.1882  
 Residual             0.9405   0.9698  
Number of obs: 240, groups:  pid, 120

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)   
(Intercept)   0.1727     0.1385 116.7062   1.247  0.21494   
con           0.3462     0.1044 227.3108   3.317  0.00106 **
ev1          -0.3439     0.2083 116.8432  -1.651  0.10143   
ev2           0.2525     0.1688 117.0168   1.495  0.13753   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
    (Intr) con    ev1   
con  0.031              
ev1  0.171 -0.049       
ev2 -0.423  0.065 -0.407

for example, I can extract fixed effect such as following. summary(MLM1)[['coefficients']]['ev1','Pr(>|t|)']

How can I extract random effect coefficients? for example, I want to extract 1.4119, 1.1882, 0.9405, 0.9698.

Random effects: Groups Name Variance Std.Dev. pid (Intercept) 1.4119 1.1882
Residual 0.9405 0.9698

Answer 1

The random effects results are not coefficients, but to get the variance and standard deviation as reported in the summary output, you can use the VarCorr function.

For example,

library(lme4)
#> Loading required package: Matrix

fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)

summary(fm1)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: Reaction ~ Days + (Days | Subject)
#>    Data: sleepstudy
#> 
#> REML criterion at convergence: 1743.6
#> 
#> Scaled residuals: 
#>     Min      1Q  Median      3Q     Max 
#> -3.9536 -0.4634  0.0231  0.4634  5.1793 
#> 
#> Random effects:
#>  Groups   Name        Variance Std.Dev. Corr
#>  Subject  (Intercept) 612.10   24.741       
#>           Days         35.07    5.922   0.07
#>  Residual             654.94   25.592       
#> Number of obs: 180, groups:  Subject, 18
#> 
#> Fixed effects:
#>             Estimate Std. Error t value
#> (Intercept)  251.405      6.825  36.838
#> Days          10.467      1.546   6.771
#> 
#> Correlation of Fixed Effects:
#>      (Intr)
#> Days -0.138

If you want the results as a table you could do:

cbind(Var    = diag(VarCorr(fm1)$Subject), 
      stddev = attr(VarCorr(fm1)$Subject, "stddev"))
#>                   Var    stddev
#> (Intercept) 612.10016 24.740658
#> Days         35.07171  5.922138

Obviously, you'll need pid instead of Subject in the code above - we don't have your data or model for a demo here.

^{Created on 2022-04-27 by the reprex package (v2.0.1)}

Answer 2

VarCorr(MLM1)$pid is the basic object.

broom.mixed::tidy(MLM1, effects = "ran_pars") may give you a more convenient format.

library(lme4)
fm1 <- lmer(Reaction ~ Days + (1|Subject), sleepstudy)
## RE variance
v1 <- VarCorr(fm1)$Subject
s1 <- attr(VarCorr(fm1)$Subject, "stddev")
## or
s1 <- sqrt(v1)
attr(VarCorr(fm1), "sc") ## residual std dev
## or
sigma(fm1)
## square these values if you want the residual variance

Or:

broom.mixed::tidy(fm1, effects = "ran_pars") ## std devs
broom.mixed::tidy(fm1, effects = "ran_pars", scales = "vcov") ## variances