Extract random effect variances from lme4 mer model object
Solution 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")
Solution 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
##
Solution 3
> attributes(summary(study))$REmat
Groups Name Variance Std.Dev.
"Subject" "(Intercept)" "1378.18" "37.124"
"Residual" "" " 960.46" "30.991"
Solution 4
Another possibility is
sum <- summary (study)
var <- data.frame (sum$varcor)
Solution 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.
Related videos on Youtube
dynamo
Updated on June 27, 2021Comments
-
dynamo about 3 years
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!
-
Dirk Eddelbuettel over 12 yearsErr, your code uses
lm()
, the questions was aboutlmer()
which is not the same thing. -
Thierry over 12 yearsIf you want the values then VarCorr() is much more efficient. Have a look at the post of Ben Bolker
-
user1320502 almost 10 yearsVarCorr seems only to provide the std deviations not the variance estimates which is what in general people would like to report right?
-
Ben Bolker almost 10 years(1) it's pretty easy to square the standard deviations; (2)
print(VarCorr(fitted_model),comp="Variance")
oras.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 -
user1320502 almost 10 yearsThanks for your comment Ben, very helpful!
-
Ben Bolker almost 10 yearsthis 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 calledmerMod
. -
blazej over 6 yearsI might be wrong, by there looks to be no
REmat
inattributes(summary(study))
anymore