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

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

CodePudding user response：

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

CodePudding user response：

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)}