Extract random effect variances from lme4 mer model object

Question

I have a mer object that has fixed and random effects. How do I extract the variance estimates for the random effects? Here is a simplified version of my question.

study <- lmer(Reaction ~ Days + (1|Subject), data = sleepstudy)
study

This gives a long output - not too long in this case. Anyway, how do I explicitly select the

Random effects:
Groups   Name        Variance Std.Dev.
Subject  (Intercept) 1378.18  37.124  
Residual              960.46  30.991  

part of the output? I want the values themselves.

I have taken long looks at

str(study)

and there's nothing there! Also checked any extractor functions in the lme4 package to no avail. Please help!

Prison

Answer 1

lmer returns an S4 object, so this should work:

remat <- summary(study)@REmat
print(remat, quote=FALSE)

Which prints:

 Groups   Name        Variance Std.Dev.
 Subject  (Intercept) 1378.18  37.124  
 Residual              960.46  30.991  

...In general, you can look at the source of the print and summary methods for "mer" objects:

class(study) # mer
selectMethod("print", "mer")
selectMethod("summary", "mer")

Horn

If you want the values then VarCorr() is much more efficient. Have a look at the post of Ben Bolker – Footcloth

this is somewhat out of date now (although the original question does refer to "mer objects", which are by definition associated with pre-1.0 lme4 -- the class is now called merMod. – Swords

Answer 2

Some of the other answers are workable, but I claim that the best answer is to use the accessor method that is designed for this -- VarCorr (this is the same as in lme4's predecessor, the nlme package).

UPDATE in recent versions of lme4 (version 1.1-7, but everything below is probably applicable to versions >= 1.0), VarCorr is more flexible than before, and should do everything you want without ever resorting to fishing around inside the fitted model object.

library(lme4)
study <- lmer(Reaction ~ Days + (1|Subject), data = sleepstudy)
VarCorr(study)
##  Groups   Name        Std.Dev.
##  Subject  (Intercept) 37.124  
##  Residual             30.991

By default VarCorr() prints standard deviations, but you can get variances instead if you prefer:

print(VarCorr(study),comp="Variance")
##  Groups   Name        Variance
##  Subject  (Intercept) 1378.18 
##  Residual              960.46 

(comp=c("Variance","Std.Dev.") will print both).

For more flexibility, you can use the as.data.frame method to convert the VarCorr object, which gives the grouping variable, effect variable(s), and the variance/covariance or standard deviation/correlations:

as.data.frame(VarCorr(study))
##        grp        var1 var2      vcov    sdcor
## 1  Subject (Intercept) <NA> 1378.1785 37.12383
## 2 Residual        <NA> <NA>  960.4566 30.99123

Finally, the raw form of the VarCorr object (which you probably shouldn't mess with you if you don't have to) is a list of variance-covariance matrices with additional (redundant) information encoding the standard deviations and correlations, as well as attributes ("sc") giving the residual standard deviation and specifying whether the model has an estimated scale parameter ("useSc").

unclass(VarCorr(fm1))
## $Subject
##             (Intercept)      Days
## (Intercept)  612.089748  9.604335
## Days           9.604335 35.071662
## attr(,"stddev")
## (Intercept)        Days 
##   24.740448    5.922133 
## attr(,"correlation")
##             (Intercept)       Days
## (Intercept)  1.00000000 0.06555134
## Days         0.06555134 1.00000000
## 
## attr(,"sc")
## [1] 25.59182
## attr(,"useSc")
## [1] TRUE
## 

Swords

VarCorr seems only to provide the std deviations not the variance estimates which is what in general people would like to report right? – Unification

(1) it's pretty easy to square the standard deviations; (2) print(VarCorr(fitted_model),comp="Variance") or as.data.frame(VarCorr(fitted_model)) will easily retrieve the variances; (3) reporting variances vs standard deviations is context dependent -- I generally prefer variances if trying to think about var decomposition/proportion explained and std devs if trying to compare to the magnitude of fixed effects – Swords

Answer 3

> attributes(summary(study))$REmat
 Groups     Name          Variance  Std.Dev.
 "Subject"  "(Intercept)" "1378.18" "37.124"
 "Residual" ""            " 960.46" "30.991"

Whicker

I might be wrong, by there looks to be no REmat in attributes(summary(study)) anymore – Cyanite

Answer 4

This package is useful for things like this (https://easystats.github.io/insight/reference/index.html)

library("insight")

get_variance_random(study) #Where study is your fit mixed model

Congenital

Answer 5

This answer is heavily based on that on @Ben Bolker's, but if people are new to this and want the values themselves, instead of just a printout of the values (as OP seems to have wanted), then you can extract the values as follows:

Convert the VarCorr object to a data frame.

re_dat = as.data.frame(VarCorr(study))

Then access each individual value:

int_vcov = re_dat[1,'vcov']
resid_vcov = re_dat[2,'vcov']

With this method (specifying rows and columns in the date frame you created) you can access whichever values you'd like.

Avis

Answer 6

Another possibility is

sum <- summary (study)
var <- data.frame (sum$varcor)

Tica

Answer 7

Try

attributes(study)

As an example:

> women
   height weight
1      58    115
2      59    117
3      60    120
4      61    123
5      62    126
6      63    129
7      64    132
8      65    135
9      66    139
10     67    142
11     68    146
12     69    150
13     70    154
14     71    159
15     72    164

> lm1 <- lm(height ~ weight, data=women)
> attributes(lm1)
$names
 [1] "coefficients"  "residuals"     "effects"       "rank"         
 [5] "fitted.values" "assign"        "qr"            "df.residual"  
 [9] "xlevels"       "call"          "terms"         "model"        

$class
[1] "lm"

> lm1$coefficients
(Intercept)      weight 
 25.7234557   0.2872492 

> lm1$coefficients[[1]]

[1] 25.72346


> lm1$coefficients[[2]]

[1] 0.2872492

The end.

Gurango

Err, your code uses lm(), the questions was about lmer() which is not the same thing. – Glutton