How to obtain the standard errors of the random effects in lme?

Question

If I fix a linear mixed effects model using R's lme from the nlme package, how do I obtain the standard errors of the random effects estimates?

For example, if lme gives the following results:

null.model <- lme(fixed = fev1 ~ 1, data = Data, random = ~ 1 | conwrd)

null.model

Linear mixed-effects model fit by REML
  Data: Dat 
  Log-restricted-likelihood: -887.7505
  Fixed: fev1 ~ 1 
(Intercept) 
   15.00424 

Random effects:
 Formula: ~1 | conwrd
        (Intercept) Residual
StdDev:    3.010589 4.130609

Number of Observations: 308
Number of Groups: 11 

How do I obtain the standard errors of the level-2 (Intercept) random effect estimate and the Residual effect estimate? For example, Stata's mixed command returns not only these estimates, but standard errors on them, and confidence interval estimates derived from these standard errors as below. NOTE: Stata reports variances, whereas R reports standard deviations, so 3.010589 and 4.130609 from the above R model output equal the square roots of 9.063698 and 17.06193 from the below Stata model output on the same data.

mixed fev1 || conwrd: , reml

[SNIP]


Mixed-effects REML regression                   Number of obs     =        308
Group variable: conwrd                          Number of groups  =         11

                                                Obs per group:
                                                              min =         25
                                                              avg =       28.0
                                                              max =         31

                                                Wald chi2(0)      =          .
Log restricted-likelihood = -887.75054          Prob > chi2       =          .

------------------------------------------------------------------------------
        fev1 |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _cons |   15.00424   .9378441    16.00   0.000      13.1661    16.84238
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
conwrd: Identity             |
                  var(_cons) |   9.063698   4.324303      3.557919    23.08952
-----------------------------+------------------------------------------------
               var(Residual) |   17.06193   1.400088      14.52711    20.03905
------------------------------------------------------------------------------
LR test vs. linear model: chibar2(01) = 94.48         Prob >= chibar2 = 0.0000

Fake data used for both these models:

conwrd <- c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11)
fev1 <- c(23, 18, 22, 19, 16, 13, 15, 17, 16, 23, 25, 18, 21, 20, 17, 21, 17, 19, 20, 21, 21, 20, 17, 15, 15, 17, 13, 14, 19, 9, 9, 11, 10, 19, 13, 16, 12, 10, 9, 11, 11, 9, 11, 10, 12, 9, 7, 8, 11, 14, 16, 13, 9, 10, 9, 8, 16, 14, 13, 9, 11, 9, 12, 12, 13, 11, 17, 16, 17, 19, 23, 24, 28, 26, 22, 25, 19, 24, 22, 23, 20, 27, 12, 12, 10, 9, 10, 9, 4, 4, 4, 8, 9, 6, 8, 6, 9, 11, 9, 8, 9, 11, 14, 17, 11, 12, 13, 10, 10, 9, 14, 13, 15, 15, 20, 12, 13, 6, 15, 16, 12, 7, 10, 7, 15, 17, 15, 18, 20, 18, 16, 21, 22, 16, 12, 15, 11, 13, 8, 17, 19, 20, 16, 20, 18, 12, 11, 8, 12, 11, 11, 16, 17, 16, 17, 17, 14, 20, 24, 24, 24, 23, 20, 21, 25, 13, 14, 14, 15, 21, 16, 17, 15, 14, 11, 8, 11, 13, 14, 13, 15, 13, 12, 15, 17, 19, 16, 14, 16, 16, 14, 14, 11, 17, 7, 10, 16, 12, 18, 18, 15, 11, 13, 9, 12, 11, 13, 9, 11, 16, 15, 15, 18, 24, 28, 24, 24, 27, 23, 23, 21, 23, 23, 22, 15, 10, 11, 13, 17, 15, 13, 10, 15, 13, 11, 13, 18, 18, 15, 22, 18, 19, 18, 20, 17, 19, 18, 14, 13, 10, 7, 11, 14, 19, 18, 15, 14, 9, 14, 15, 14, 19, 18, 14, 10, 17, 23, 25, 26, 24, 24, 26, 25, 25, 20, 20, 20, 20, 17, 15, 14, 12, 11, 11, 11, 11, 9, 10, 11, 13, 13, 17, 16, 11, 11, 11, 12, 19, 15, 13, 15, 15, 12, 9, 12, 10, 8, 8)
Data <- data.frame(conwrd,fev1)

Answer 1

The answer is very well explained here https://stats.oarc.ucla.edu/r/faq/how-can-i-calculate-standard-errors-for-variance-components-from-mixed-models/ based on standard errors of variance of random effects using Fisher information matrix from the package lmeInfo. Basically, the SE is acquired from the Hessian matrix (second-order derivative) of estimated parameters.

However, confidence interval is random effect standard deviation or variance encouraged to report over SE, as the latter usually intrigues readers to conduct hypothesis testing in a t-test style, which is dangerous in this case of a parameter of standard deviation or variance.

Another way, of course, is to obtain standard errors from standard deviation of estimates across bootstrapped samples. This may require a customization of functions from the boot packages, as none of the available packages supposedly designed for longitudinal model bootstrapping that I tested does resampling in the correct way.