Standard Error of variance component from the output of lmer

I need to extract the standard error of variance component from the output of lmer .

library(lme4)
model <- lmer(Reaction ~ Days + (1|Subject), sleepstudy)

The following produces estimates of variance component :

s2 <- VarCorr(model)$Subject[1]

It is NOT the standard error of the variance. And I want the standard error . How can I have it ?

EDIT :

Perhaps I am unable to make you understand what I meant by "standard error of the variance component". So I am editing my post .

In Chapter 12 , Experiments with Random Factors , of the book Design and Analysis of Experiments , by Douglas C. Montgomery , at the end of the chapter , Example 12-2 is done by SAS . In Example 12-2 , the model is a two-factor factorial random effect model .The output is given in Table 12-17

I am trying to fit the model in R by lmer .

library(lme4)
fit <- lmer(y~(1|operator)+(1|part),data=dat)

R codes for extracting the Estimate , annotated by 4 in the table 12-17 :

est_ope=VarCorr(fit)$operator[1]
est_part = VarCorr(fit)$part[1]
sig = summary(fit)$sigma
est_res = sig^2

Now I want to extract the results of Std Errors , annotated by 5 in the table 12-17 from lmer output .

Many Thanks !

I think you are looking for the Wald standard error of the variance estimates. Please note that these (as often pointed out by Doug Bates) the Wald standard errors are often very poor estimates of the uncertainty of variances, because the likelihood profiles are often far from quadratic on the variance scale ... I'm assuming you know what you're doing and have some good use for these numbers ...

This can be (now) done using the merDeriv package.

library(lme4)
library(merDeriv)
m1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
sqrt(diag(vcov(m1, full = TRUE)))
vv <- vcov(m1, full = TRUE)
colnames(vv)
## [1] "(Intercept)"                  "Days"                        
## [3] "cov_Subject.(Intercept)"      "cov_Subject.Days.(Intercept)"
## [5] "cov_Subject.Days"             "residual"

Because this represents the full/combined covariance matrix, the first two indices (rows/columns) {(Intercept), Days} represent the fixed-effect intercept and slope. The elements prefaced by cov_ are variance components, in the format cov_<grp_variable>.<term> for variances (e.g. cov_Subject.(Intercept) is the among-subject variance in the intercept) and cov_<grp_variable>.<term1>.<term2> for covariances (cov_Subject.Days.(Intercept) is the among-subject covariance between the slope (Days) and the intercept). residual denotes the residual variance.

If we want the standard errors of the variance components, we take the square root of the diagonal and keep only elements 3 to 5:

sqrt(diag(vv)[3:5])
## [1] 288.78602  46.67876  14.78208

or more generally

sqrt(diag(vv)[grepl("^cov", colnames(vv))])

(weirdly, only colnames() works — the row names of vv are blank for the random effects components).

old answer

library("lme4")
model <- lmer(Reaction ~ Days + (1|Subject), sleepstudy, REML=FALSE)

(At present it's quite a bit harder to do this for the REML estimates ...)

Extract deviance function parameterized in terms of standard deviation and correlation rather than in terms of Cholesky factors (note this an internal function, so there's not a guarantee that it will keep working in the same way in the future ...)

 dd.ML <- lme4:::devfun2(model,useSc=TRUE,signames=FALSE)

Extract parameters as standard deviations on original scale:

 vv <- as.data.frame(VarCorr(model)) ## need ML estimates!
 pars <- vv[,"sdcor"]
 ## will need to be careful about order if using this for
 ## a random-slopes model ...

Now compute the second-derivative (Hessian) matrix:

library("numDeriv")
hh1 <- hessian(dd.ML,pars)
vv2 <- 2*solve(hh1)  ## 2* converts from log-likelihood to deviance scale
sqrt(diag(vv2))  ## get standard errors

These are the standard errors of the standard deviations: double them to get the standard errors of the variances (when you transform a value, its standard errors scale according to the derivative of the transformation).

I think this should do it, but you might want to double-check it ...