Linear mixed model covariance decomposition with random intercept- lme4
Daniel
2022-06-17
Linear mixed model covariance decomposition with random intercept- lme4
Load data
data(Orthodont, package="nlme")
Data <- Orthodontplot means and variances by higher level variable (grouping)
barplot(with(Data, tapply(distance, list(subject = Subject), mean)))
barplot(with(Data, tapply(distance, list(subject = Subject), var)))
Using glm
fit.lm <- lm(distance ~ age+ Sex + Subject, data = Data)
anova(fit.lm)## Analysis of Variance Table
##
## Response: distance
## Df Sum Sq Mean Sq F value Pr(>F)
## age 1 235.36 235.356 114.8383 < 2.2e-16 ***
## Sex 1 140.46 140.465 68.5376 2.297e-12 ***
## Subject 25 377.91 15.117 7.3759 3.081e-12 ***
## Residuals 80 163.96 2.049
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Using glm with random slopes
fit.aov <- aov(distance ~ age+ Sex + Error(Subject), data = Data)
summary(fit.aov)##
## Error: Subject
## Df Sum Sq Mean Sq F value Pr(>F)
## Sex 1 140.5 140.46 9.292 0.00538 **
## Residuals 25 377.9 15.12
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## age 1 235.4 235.36 114.8 <2e-16 ***
## Residuals 80 164.0 2.05
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1how can we compute \(σ^2\) a from this output? note that \(σ^2\)= MSE , and E(MSA) = \((σ_a)^2\) +n·\(σ^2\) ,therefore\((σ_a)^2\) = (MSA−MSE)/n
Using lmm
library(lme4) ## Loading required package: Matrixfit.lmer <- lmer(distance ~ age + Sex +(1 | Subject), data = Data)
summary(fit.lmer)## Linear mixed model fit by REML ['lmerMod']
## Formula: distance ~ age + Sex + (1 | Subject)
## Data: Data
##
## REML criterion at convergence: 437.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7489 -0.5503 -0.0252 0.4534 3.6575
##
## Random effects:
## Groups Name Variance Std.Dev.
## Subject (Intercept) 3.267 1.807
## Residual 2.049 1.432
## Number of obs: 108, groups: Subject, 27
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 17.70671 0.83392 21.233
## age 0.66019 0.06161 10.716
## SexFemale -2.32102 0.76142 -3.048
##
## Correlation of Fixed Effects:
## (Intr) age
## age -0.813
## SexFemale -0.372 0.000default is compound symmetry
Get x,y,z
X=(getME(fit.lmer, "X"))
head(X)## (Intercept) age SexFemale
## 1 1 8 0
## 2 1 10 0
## 3 1 12 0
## 4 1 14 0
## 5 1 8 0
## 6 1 10 0y=getME(fit.lmer, "y")
head(y)## [1] 26.0 25.0 29.0 31.0 21.5 22.5Z <- getME(fit.lmer, "Z")
head(Z)## 6 x 27 sparse Matrix of class "dgCMatrix"## [[ suppressing 27 column names 'M16', 'M05', 'M02' ... ]]##
## 1 . . . . . . . . . . . . . . 1 . . . . . . . . . . . .
## 2 . . . . . . . . . . . . . . 1 . . . . . . . . . . . .
## 3 . . . . . . . . . . . . . . 1 . . . . . . . . . . . .
## 4 . . . . . . . . . . . . . . 1 . . . . . . . . . . . .
## 5 . . 1 . . . . . . . . . . . . . . . . . . . . . . . .
## 6 . . 1 . . . . . . . . . . . . . . . . . . . . . . . .Get fixed effect and random effect coefficients
bhat <- getME(fit.lmer, "fixef") #getME(fit.lmer, "beta") # fixef(fit.lmer)
bhat ## (Intercept) age SexFemale
## 17.7067130 0.6601852 -2.3210227uhat <- ranef(fit.lmer)
uhat## $Subject
## (Intercept)
## M16 -1.70183357
## M05 -1.70183357
## M02 -1.37767479
## M11 -1.16156894
## M07 -1.05351602
## M08 -0.94546309
## M03 -0.62130432
## M12 -0.62130432
## M13 -0.62130432
## M14 -0.08103969
## M09 0.13506616
## M15 0.78338371
## M06 1.21559540
## M04 1.43170125
## M01 2.40417758
## M10 3.91691853
## F10 -3.58539251
## F09 -1.31628108
## F06 -1.31628108
## F01 -1.10017524
## F05 -0.01964599
## F07 0.30451279
## F02 0.30451279
## F08 0.62867156
## F03 0.95283034
## F04 1.92530666
## F11 3.22194176
##
## with conditional variances for "Subject"coef(fit.lmer)## $Subject
## (Intercept) age SexFemale
## M16 16.00488 0.6601852 -2.321023
## M05 16.00488 0.6601852 -2.321023
## M02 16.32904 0.6601852 -2.321023
## M11 16.54514 0.6601852 -2.321023
## M07 16.65320 0.6601852 -2.321023
## M08 16.76125 0.6601852 -2.321023
## M03 17.08541 0.6601852 -2.321023
## M12 17.08541 0.6601852 -2.321023
## M13 17.08541 0.6601852 -2.321023
## M14 17.62567 0.6601852 -2.321023
## M09 17.84178 0.6601852 -2.321023
## M15 18.49010 0.6601852 -2.321023
## M06 18.92231 0.6601852 -2.321023
## M04 19.13841 0.6601852 -2.321023
## M01 20.11089 0.6601852 -2.321023
## M10 21.62363 0.6601852 -2.321023
## F10 14.12132 0.6601852 -2.321023
## F09 16.39043 0.6601852 -2.321023
## F06 16.39043 0.6601852 -2.321023
## F01 16.60654 0.6601852 -2.321023
## F05 17.68707 0.6601852 -2.321023
## F07 18.01123 0.6601852 -2.321023
## F02 18.01123 0.6601852 -2.321023
## F08 18.33538 0.6601852 -2.321023
## F03 18.65954 0.6601852 -2.321023
## F04 19.63202 0.6601852 -2.321023
## F11 20.92865 0.6601852 -2.321023
##
## attr(,"class")
## [1] "coef.mer"Get random effect covariance structure (covariance matrix expandation)
# theta <- getME(fit.lmer, "theta") # $r$
# theta**2+ getME(fit.lmer, "sigma")**2
#
# theta**2
# lwr <- getME(fit.lmer, "lower")
# lwr
# Lambda <- getME(fit.lmer, "Lambda") # matrix
# # head(Lambda)
# # dim(Lambda)
# Lambda%*%t(Lambda)the relationship $\Sigma = TSST'$
the relationship of theta and covariance $$ = (
\[\begin{array}{cc}
\theta_1 & 0 \\
\theta_2 & \theta_3
\end{array}\]
) (
\[\begin{array}{cc}
\theta_1 & \theta_2 \\
0 & \theta_3
\end{array}\]
) = (
\[\begin{array}{cc}
\theta_1^2 & \theta_1 \theta_2 \\
\theta_1 \theta_2 & \theta_2^2 + \theta_3^2
\end{array}\]
) = (
\[\begin{array}{cc}
\sigma_1^2 & \sigma_{12} \\
\sigma_{12} & \sigma_2^2
\end{array}\]
)
$$
we can equate θ1 with σ1 (the standard deviation of the intercept), θ2 with σ12/σ1 (the covariance scaled by the SD) ... θ3 is a little more complicated. in particular, in a model with only random-intercept terms, the θ vector is just the vector of RE standard deviations.
The way lme4 works is that, given a value of θ, we can determine Σ
vc <- VarCorr(fit.lmer)
Lambda_new <-vc[["Subject"]][1]*diag(length(levels(Data$Subject)))
head(Lambda_new)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 3.266784 0.000000 0.000000 0.000000 0.000000 0.000000 0 0 0 0
## [2,] 0.000000 3.266784 0.000000 0.000000 0.000000 0.000000 0 0 0 0
## [3,] 0.000000 0.000000 3.266784 0.000000 0.000000 0.000000 0 0 0 0
## [4,] 0.000000 0.000000 0.000000 3.266784 0.000000 0.000000 0 0 0 0
## [5,] 0.000000 0.000000 0.000000 0.000000 3.266784 0.000000 0 0 0 0
## [6,] 0.000000 0.000000 0.000000 0.000000 0.000000 3.266784 0 0 0 0
## [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22]
## [1,] 0 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0 0 0 0 0 0
## [,23] [,24] [,25] [,26] [,27]
## [1,] 0 0 0 0 0
## [2,] 0 0 0 0 0
## [3,] 0 0 0 0 0
## [4,] 0 0 0 0 0
## [5,] 0 0 0 0 0
## [6,] 0 0 0 0 0dim(Lambda_new)## [1] 27 27“theta”: random-effects parameter estimates: these are parameterized as the relative Cholesky factors of each random effect term
use other cov structure by nlme
lme4 does not allow to specify it, and it can only fit LMMs with independent residual errors. However, the range of available variance-covariance matrices for the random effects are restricted to diagonal or general matrices.
Get residual variance and its identity matrix
sigma <- getME(fit.lmer, "sigma")**2
sigma## [1] 2.049456head(sigma*diag(nrow(Data)))## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 2.049456 0.000000 0.000000 0.000000 0.000000 0.000000 0 0 0 0
## [2,] 0.000000 2.049456 0.000000 0.000000 0.000000 0.000000 0 0 0 0
## [3,] 0.000000 0.000000 2.049456 0.000000 0.000000 0.000000 0 0 0 0
## [4,] 0.000000 0.000000 0.000000 2.049456 0.000000 0.000000 0 0 0 0
## [5,] 0.000000 0.000000 0.000000 0.000000 2.049456 0.000000 0 0 0 0
## [6,] 0.000000 0.000000 0.000000 0.000000 0.000000 2.049456 0 0 0 0
## [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22]
## [1,] 0 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0 0 0 0 0 0
## [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34]
## [1,] 0 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0 0 0 0 0 0
## [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46]
## [1,] 0 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0 0 0 0 0 0
## [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55] [,56] [,57] [,58]
## [1,] 0 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0 0 0 0 0 0
## [,59] [,60] [,61] [,62] [,63] [,64] [,65] [,66] [,67] [,68] [,69] [,70]
## [1,] 0 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0 0 0 0 0 0
## [,71] [,72] [,73] [,74] [,75] [,76] [,77] [,78] [,79] [,80] [,81] [,82]
## [1,] 0 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0 0 0 0 0 0
## [,83] [,84] [,85] [,86] [,87] [,88] [,89] [,90] [,91] [,92] [,93] [,94]
## [1,] 0 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0 0 0 0 0 0
## [,95] [,96] [,97] [,98] [,99] [,100] [,101] [,102] [,103] [,104] [,105]
## [1,] 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0 0 0 0 0
## [,106] [,107] [,108]
## [1,] 0 0 0
## [2,] 0 0 0
## [3,] 0 0 0
## [4,] 0 0 0
## [5,] 0 0 0
## [6,] 0 0 0Get y covariance matrix
VM <- Z%*%Lambda_new%*%t(Z)+sigma*diag(nrow(Data))
head(VM)## 6 x 108 sparse Matrix of class "dgCMatrix"## [[ suppressing 108 column names '1', '2', '3' ... ]]##
## 1 5.316240 3.266784 3.266784 3.266784 . . . . . . .
## 2 3.266784 5.316240 3.266784 3.266784 . . . . . . .
## 3 3.266784 3.266784 5.316240 3.266784 . . . . . . .
## 4 3.266784 3.266784 3.266784 5.316240 . . . . . . .
## 5 . . . . 5.316240 3.266784 3.266784 3.266784 . . .
## 6 . . . . 3.266784 5.316240 3.266784 3.266784 . . .
##
## 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
##
## 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
##
## 1 . . . . . . . . . . . . . . . . . . .
## 2 . . . . . . . . . . . . . . . . . . .
## 3 . . . . . . . . . . . . . . . . . . .
## 4 . . . . . . . . . . . . . . . . . . .
## 5 . . . . . . . . . . . . . . . . . . .
## 6 . . . . . . . . . . . . . . . . . . .dim(VM)## [1] 108 108Get fixed effect coefficients covariance matrix
vcov <- vcov(fit.lmer) #fixed cov
vcov## 3 x 3 Matrix of class "dpoMatrix"
## (Intercept) age SexFemale
## (Intercept) 0.69542669 -4.174818e-02 -2.361967e-01
## age -0.04174818 3.795289e-03 -3.051416e-17
## SexFemale -0.23619673 -3.051416e-17 5.797556e-01# computer correlation coefficients
vcov@x[2]/prod(sqrt( diag(vcov(fit.lmer))[-3] ))## [1] -0.8126236Get random effect coefficients covariance matrix (standard deviation)
vc <- VarCorr(fit.lmer)
vc ## default print method: standard dev and corr## Groups Name Std.Dev.
## Subject (Intercept) 1.8074
## Residual 1.4316as.matrix(Matrix::bdiag(vc)) #random effect covariance matrix## (Intercept)
## (Intercept) 3.266784# https://stackoverflow.com/questions/47307340/extracting-the-i-estimated-variance-covariance-matrix-of-random-effects-and-or
# https://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_mixed_sect022.htmapproximate above sd
uintercept <- (uhat[["Subject"]])
sd(uintercept[,1])## [1] 1.647809Computer fixed effect coefficients
library(matlib)
inv(t(as.matrix(X))%*%inv(as.matrix(VM))%*%as.matrix(X))%*%t(as.matrix(X))%*%inv(as.matrix(VM))%*%y## [,1]
## [1,] 17.7067076
## [2,] 0.6601871
## [3,] -2.3210230bhat## (Intercept) age SexFemale
## 17.7067130 0.6601852 -2.3210227Computer covariance of fixed efffect coefficients
inv(t(as.matrix(X))%*%inv(as.matrix(VM))%*%as.matrix(X))##
## [1,] 0.69542669 -0.04174818 -0.2361967
## [2,] -0.04174818 0.00379529 0.0000000
## [3,] -0.23619674 0.00000000 0.5797556# standard error
sqrt(inv(t(as.matrix(X))%*%inv(as.matrix(VM))%*%as.matrix(X)))## Warning in sqrt(inv(t(as.matrix(X)) %*% inv(as.matrix(VM)) %*% as.matrix(X))):
## NaNs produced##
## [1,] 0.8339225 NaN NaN
## [2,] NaN 0.06160593 0.0000000
## [3,] NaN 0.00000000 0.7614169# the following equals lmm summary
vcov(fit.lmer) ## 3 x 3 Matrix of class "dpoMatrix"
## (Intercept) age SexFemale
## (Intercept) 0.69542669 -4.174818e-02 -2.361967e-01
## age -0.04174818 3.795289e-03 -3.051416e-17
## SexFemale -0.23619673 -3.051416e-17 5.797556e-01sqrt(vcov(fit.lmer))## Warning in sqrt(x@x): NaNs produced## 3 x 3 Matrix of class "dsyMatrix"
## (Intercept) age SexFemale
## (Intercept) 0.8339225 NaN NaN
## age NaN 0.06160592 NaN
## SexFemale NaN NaN 0.7614168default is compound symmetry
Computer random effect coefficients
comput_uhat <- (as.matrix(Lambda_new))%*%t(Z)%*%inv(as.matrix(VM))%*%(y-as.matrix(X)%*%(bhat))
cbind((comput_uhat@x),(uhat[["Subject"]]))## (comput_uhat@x) (Intercept)
## M16 -1.70183353 -1.70183357
## M05 -1.70183353 -1.70183357
## M02 -1.37767476 -1.37767479
## M11 -1.16156892 -1.16156894
## M07 -1.05351600 -1.05351602
## M08 -0.94546307 -0.94546309
## M03 -0.62130431 -0.62130432
## M12 -0.62130431 -0.62130432
## M13 -0.62130431 -0.62130432
## M14 -0.08103969 -0.08103969
## M09 0.13506615 0.13506616
## M15 0.78338369 0.78338371
## M06 1.21559538 1.21559540
## M04 1.43170123 1.43170125
## M01 2.40417753 2.40417758
## M10 3.91691845 3.91691853
## F10 -3.58539243 -3.58539251
## F09 -1.31628106 -1.31628108
## F06 -1.31628106 -1.31628108
## F01 -1.10017521 -1.10017524
## F05 -0.01964599 -0.01964599
## F07 0.30451278 0.30451279
## F02 0.30451278 0.30451279
## F08 0.62867155 0.62867156
## F03 0.95283032 0.95283034
## F04 1.92530662 1.92530666
## F11 3.22194169 3.22194176- covariance matrix of random effect coefficients
head(Lambda_new)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 3.266784 0.000000 0.000000 0.000000 0.000000 0.000000 0 0 0 0
## [2,] 0.000000 3.266784 0.000000 0.000000 0.000000 0.000000 0 0 0 0
## [3,] 0.000000 0.000000 3.266784 0.000000 0.000000 0.000000 0 0 0 0
## [4,] 0.000000 0.000000 0.000000 3.266784 0.000000 0.000000 0 0 0 0
## [5,] 0.000000 0.000000 0.000000 0.000000 3.266784 0.000000 0 0 0 0
## [6,] 0.000000 0.000000 0.000000 0.000000 0.000000 3.266784 0 0 0 0
## [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22]
## [1,] 0 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0 0 0 0 0 0
## [,23] [,24] [,25] [,26] [,27]
## [1,] 0 0 0 0 0
## [2,] 0 0 0 0 0
## [3,] 0 0 0 0 0
## [4,] 0 0 0 0 0
## [5,] 0 0 0 0 0
## [6,] 0 0 0 0 0Computer predicted values
yhat <- X%*%(bhat)+Z%*%(uhat[["Subject"]][["(Intercept)"]])
head(yhat)## 6 x 1 Matrix of class "dgeMatrix"
## [,1]
## 1 25.39237
## 2 26.71274
## 3 28.03311
## 4 29.35348
## 5 21.61052
## 6 22.93089head(fitted(fit.lmer))## 1 2 3 4 5 6
## 25.39237 26.71274 28.03311 29.35348 21.61052 22.93089