The standard errors of variance components in a mixed-effects model can provide valuable information about the contribution of the random effects to the model. Typically, the reported parameter of a random effect is the standard deviation of the random intercepts or random slopes. R presents these standard deviations, but does not report their standard errors. The standard errors of a random effects parameter, if very large, can be a red flag suggesting a problem with the model specification or data. Otherwise, these values indicate how certain you are of your parameter values indicating how groups or subjects differ in their intercepts or slopes.
When fitting a mixed-effects model in R using the nlme package, the information provided in the
summary command includes a section for random effects. Below, we
use an example dataset from Singer and Willet’s Applied Longitudinal Data Analysis.
The random effects output is italicized.
alcohol1 <- read.table("https://stats.idre.ucla.edu/stat/r/examples/alda/data/alcohol1_pp.txt", header=T, sep=",")
attach(alcohol1)
library(nlme)
library(msm)
library(lmeInfo)
model.c <- lme(alcuse ~ coa*age_14 , data=alcohol1, random= ~ age_14 | id, method="ML")
summary(model.c)
Linear mixed-effects model fit by maximum likelihood
Data: alcohol1
AIC BIC logLik
637.2026 665.2453 -310.6013
Random effects:
Formula: ~age_14 | id
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 0.6982690 (Intr)
age_14 0.3880684 -0.219
Residual 0.5807688
Fixed effects: alcuse ~ coa * age_14
Value Std.Error DF t-value p-value
(Intercept) 0.3159517 0.13177099 162 2.397734 0.0176
coa 0.7432120 0.19616697 80 3.788671 0.0003
age_14 0.2929552 0.08492089 162 3.449742 0.0007
coa:age_14 -0.0494299 0.12642140 162 -0.390993 0.6963
Correlation:
(Intr) coa age_14
coa -0.672
age_14 -0.460 0.309
coa:age_14 0.309 -0.460 -0.672
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.5480533 -0.3879877 -0.1057547 0.3602391 2.3960739
Number of Observations: 246
Number of Groups: 82
While the standard errors of these estimated standard deviations are not
reported, they can be generated using the intervals command.
intervals(model.c)
Approximate 95% confidence intervals
Fixed effects:
lower est. upper
(Intercept) 0.05786567 0.31595172 0.5740378
coa 0.35601420 0.74321203 1.1304099
age_14 0.12662951 0.29295517 0.4592808
coa:age_14 -0.29703830 -0.04942994 0.1981784
attr(,"label")
[1] "Fixed effects:"
Random Effects:
Level: id
lower est. upper
sd((Intercept)) 0.5397860 0.6982690 0.9032833
sd(age_14) 0.2686708 0.3880684 0.5605266
cor((Intercept),age_14) -0.5622649 -0.2189939 0.1886537
Within-group standard error:
lower est. upper
0.4982931 0.5807688 0.6768956
Note that the intervals for the random effects standard deviations are NOT symmetric about the estimate. Because standard deviations must be non-negative, the actual model-estimated value is the log of the standard deviation. This intervals presented are based on the log(sd) scale. We can see this by looking one random effect, sd((Intercept)), and noting the symmetry of the logged interval and estimate values.
re1 <- c(0.5397860, 0.6982690, 0.9032833) log(re1)[2] - log(re1)[1] [1] 0.2574317 log(re1)[3] - log(re1)[2] [1] 0.2574318
These differences can be divided by 1.96 to find the standard error in the log(sd) scale. If we wish to calculate standard errors in the standard deviation scale, we can use the delta method and the variance-covariance matrix of these random effects parameters. To see the variance-covariance matrix of these parameters, we can look at the apVar object of our model and then the “Pars” attribute within that.
var <-model.c$apVar
var
reStruct.id1 reStruct.id2 reStruct.id3 lSigma
reStruct.id1 0.017251545 0.008342683 -0.03183090 -0.003538851
reStruct.id2 0.008342683 0.035194863 -0.04416056 -0.006864807
reStruct.id3 -0.031830896 -0.044160558 0.17807321 0.011314163
lSigma -0.003538851 -0.006864807 0.01131416 0.006106870
attr(,"Pars")
reStruct.id1 reStruct.id2 reStruct.id3 lSigma
-0.3591508 -0.9465736 -0.4451982 -0.5434025
attr(,"natural")
[1] TRUE
par<-attr(var, "Pars")
par
reStruct.id1 reStruct.id2 reStruct.id3 lSigma
-0.3591508 -0.9465736 -0.4451982 -0.5434025
These are logged standard deviations, so we will transform them to variances:
vc<-exp(par)^2 vc reStruct.id1 reStruct.id2 reStruct.id3 lSigma 0.4875796 0.1505971 0.4104930 0.3372924
We can square the standard deviations in our random effects output to match the first, second, and fourth values in this vector. The third value relates to the correlation of the random intercepts and random slopes.
0.6982690*0.6982690 [1] 0.4875796 0.3880684*0.3880684 [1] 0.1505971 0.5807688*0.5807688 [1] 0.3372924
Thus, to estimate the standard errors of these variances, we can use the delta method with the variance/covariance matrix entries saved as var, the list of untransformed random effects parameters saved as par, and the indicated transform, exp(x)^2. For more on the delta method in R, see FAQ: How can I estimate the standard error of a transformed parameter in R using the delta method?.
deltamethod (~ exp(x1)^2, par, var) [1] 0.1280824 deltamethod (~ exp(x2)^2, par, var) [1] 0.05650492 deltamethod (~ exp(x4)^2, par, var) [1] 0.05271642
These values are the standard errors of the variances of the random intercept, random slope, and model residuals from our model.
Note that the variance covariance matrix of the log transformed of the standard deviations of random effects, var, are already approximated using delta method and we are using delta method one more time to approximate the standard errors of the variances of random components. This might not be the most accurate and effective way. We can extract the standard errors of variance of random effects directly using fisher information matrix from the package lmeInfo.
I < Fisher_info(model.c, type = "expected") sqrt(diag(solve(I))) Tau.id.var((Intercept)) Tau.id.cov(age_14,(Intercept)) Tau.id.var(age_14) sigma_sq 0.12781822 0.06572766 0.05638669 0.05267622
More directly we can use the variance covariance of variance components.
sqrt(diag(varcomp_vcov(model.c))) Tau.id.var((Intercept)) Tau.id.cov(age_14,(Intercept)) Tau.id.var(age_14) sigma_sq 0.12781822 0.06572766 0.05638669 0.05267622