3

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)
Alexis
  • 784
  • 8
  • 32
  • Can you supply a [reproducible example](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) with sample input data and the code you ran? You should also supply the value that STATA returns for those values so possible solutions can be tested. – MrFlick Feb 17 '16 at 20:57
  • @MrFlick TY. Code and Stata output added. Data not provided (large and not public), but this is less a question about specific output values than about the method of obtaining a specific *kind* of output. :) – Alexis Feb 17 '16 at 21:13
  • What you have added doesn't help much. You can (and should) include sample data rather than your real data. If it's not reproducible, it doesn't make it any easier to help you. – MrFlick Feb 17 '16 at 22:07
  • @MrFlick OK! Done. (Caveat: If what I did doesn't help much, why did you request it?) – Alexis Feb 18 '16 at 00:56
  • I requested a reproducible example. Your first edit did not help with that because it did not have data. Today i don't have as much time to look into this. There is a discussion of a similar topic on [stats.se]: http://stats.stackexchange.com/questions/48254/standard-error-of-random-effects-in-r-lme4-vs-stata-xtmixed It depends on how you want to estimate the standard errors of the parameters. – MrFlick Feb 18 '16 at 18:59
  • @MrFlick TY! I missed that other question, although they are asking about SEs on the *fixed* effects, not the *random* effects. Stata reports them, `lme` seemingly does not. :( – Alexis Feb 18 '16 at 23:00
  • Have you tried `intervals` (for CI, not SE)? – aosmith Feb 24 '16 at 21:11
  • @aosmith intervals does not report the SEs. I want the SEs. – Alexis Feb 25 '16 at 00:15
  • I believe the covariance matrix of the log variances can be found in `model$apVar`. You could probably use the delta method on these to get the variances of the variances. – aosmith Apr 07 '16 at 19:06
  • @aosmith That's pretty weak sauce (not your comment, but the seeming lack of functionality of the programs). – Alexis Apr 07 '16 at 22:45
  • 1
    [This answer](http://stackoverflow.com/a/31704646/2461552) seems related to what you are after, and has a nice caveat about the weak sauce that is the Wald standard errors of variance components. The reason they aren't easy to get in `lme` is likely because of how the original package authors viewed them. – aosmith Apr 08 '16 at 14:37
  • 2
    See also [this question](https://stats.stackexchange.com/questions/161225/estimates-of-the-variance-of-the-variance-component-of-a-mixed-effects-model/189208) on CrossValidated. – Livius Jan 15 '19 at 21:43

1 Answers1

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. However, CI is encouraged to report over SE for random effect standard deviation or variance.

DrJerryTAO
  • 21
  • 3