Decomposition of Linear Mixed Models with Sliding and Sum Contrasts
Reinhold Kliegl
2018-06-09 (revised: 2018-09-27 21:12:40)
- 1 Overview
- 2 General forms for sliding contrast
- 3 General forms for sum contrast
- 4 A few co(variance) rules
- 5 Illustrations for two levels (k = 2)
- 6 Illustration for three levels (k = 3)
- 7 Model comparisons
- 8 Images of random-effect structures
- 9 Conclusion
- 10 References
- 11 Acknowledgment
- 12 Appendix
1 Overview
In this script we decompose LMMs with sliding and sum contrasts for a fixed factor with 2 or 3 levels. The decomposition for the treatment contrast is already available. As an illustration we use the Machines data from package MEMMS. This is a follow-up to an RPub about a decomposition of linear mixed models with treatmen contrasts (Kliegl, 2018f).
2 General forms for sliding contrast
Here are the \(\Sigma\)s for sliding contrasts for factors with two and three levels.
LMM m1, k=2
\[\mathbf\Sigma = \begin{bmatrix} \sigma_1^2 + \sigma_2^2/4 - \sigma_{12} & \sigma_1^2 - \sigma_2^2/4 \\ \sigma_1^2 - \sigma_2^2/4 & \sigma_1^2 + \sigma_2^2/4 + \sigma_{12} \\ \end{bmatrix} = \begin{bmatrix} \sigma_1^2 & \sigma_1^2 \\ \sigma_1^2 & \sigma_1^2 \\ \end{bmatrix} + 1/4 \begin{bmatrix} \sigma_2^2 & - \sigma_2^2 \\ - \sigma_2^2 & \sigma_2^2 \\ \end{bmatrix} + \begin{bmatrix} - \sigma_{12} & 0 \\ 0 & \sigma_{12} \\ \end{bmatrix} \]
LMM m3, k=2
\[\mathbf\Sigma = \begin{bmatrix} \sigma_1^2 & \sigma_1^2 \\ \sigma_1^2 & \sigma_1^2 \\ \end{bmatrix} + 1/4 \begin{bmatrix} \sigma_2^2 & - \sigma_2^2 \\ - \sigma_2^2 & \sigma_2^2 \\ \end{bmatrix}\]
LMM m1, k=3
\[\mathbf\Sigma = 1/9 \begin{bmatrix} 9 \sigma_1^2 + 4 \sigma_2^2 + \sigma_3^2 - 12 \sigma_{12} - 6 \sigma_{13} + 4 \sigma_{23} & 9 \sigma_1^2 - 2 \sigma_2^2 + \sigma_3^2 - 3 \sigma_{12} - 6 \sigma_{13} + \sigma_{23} & 9 \sigma_1^2 - 2 \sigma_2^2 - 2 \sigma_3^2 - 3 \sigma_{12} + 3 \sigma_{13} - 5 \sigma_{23}\\ 9 \sigma_1^2 - 2 \sigma_2^2 + \sigma_3^2 - 3 \sigma_{12} - 6 \sigma_{13} + \sigma_{23} & 9 \sigma_1^2 + \sigma_2^2 + \sigma_3^2 + 6 \sigma_{12} - 6 \sigma_{13} - 2 \sigma_{23} & 9 \sigma_1^2 + \sigma_2^2 - 2 \sigma_3^2 + 6 \sigma_{12} + 3 \sigma_{13} + \sigma_{23}\\ 9 \sigma_1^2 - 2 \sigma_2^2 - 2 \sigma_3^2 - 3 \sigma_{12} + 3 \sigma_{13} - 5 \sigma_{23} & 9 \sigma_1^2 + \sigma_2^2 - 2 \sigma_3^2 + 6 \sigma_{12} + 3 \sigma_{13} + \sigma_{23} & 9 \sigma_1^2 + \sigma_2^2 + 4 \sigma_3^2 + 6 \sigma_{12} + 12 \sigma_{13} +4 \sigma_{23}\\ \end{bmatrix}\]
This matrix corresponds to the following decomposition:
\[\mathbf \Sigma = \begin{bmatrix} \sigma_1^2 & \sigma_1^2 & \sigma_1^2 \\ \sigma_1^2 & \sigma_1^2 & \sigma_1^2 \\ \sigma_1^2 & \sigma_1^2 & \sigma_1^2 \\ \end{bmatrix} + 1/9 \begin{bmatrix} 4 \sigma_2^2 + \sigma_3^2 & - 2 \sigma_2^2 + \sigma_3^2 & - 2 \sigma_2^2 - 2 \sigma_3^2 \\ - 2 \sigma_2^2 + \sigma_3^2 & \sigma_2^2 + \sigma_3^2 & \sigma_2^2 - 2 \sigma_3^2\\ - 2 \sigma_2^2 - 2 \sigma_3^2 & \sigma_2^2 - 2 \sigma_3^2 & \sigma_2^2 + 4 \sigma_3^2\\ \end{bmatrix} + 1/9 \begin{bmatrix} - 12 \sigma_{12} - 6 \sigma_{13} + 4 \sigma_{23} & - 3 \sigma_{12} - 6 \sigma_{13} + \sigma_{23} & - 3 \sigma_{12} + 3 \sigma_{13} - 5 \sigma_{23}\\ - 3 \sigma_{12} - 6 \sigma_{13} + \sigma_{23} & 6 \sigma_{12} - 6 \sigma_{13} - 2 \sigma_{23} & 6 \sigma_{12} + 3 \sigma_{13} + \sigma_{23}\\ - 3 \sigma_{12} + 3 \sigma_{13} - 5 \sigma_{23} & 6 \sigma_{12} + 3 \sigma_{13} + \sigma_{23} & 6 \sigma_{12} + 12 \sigma_{13} +4 \sigma_{23}\\ \end{bmatrix}\]
LMM m3, k=3
\[\mathbf \Sigma = \begin{bmatrix} \sigma_1^2 & \sigma_1^2 & \sigma_1^2 \\ \sigma_1^2 & \sigma_1^2 & \sigma_1^2 \\ \sigma_1^2 & \sigma_1^2 & \sigma_1^2 \\ \end{bmatrix} + 1/9 \begin{bmatrix} 4 \sigma_2^2 + \sigma_3^2 & - 2 \sigma_2^2 + \sigma_3^2 & - 2 \sigma_2^2 - 2 \sigma_3^2 \\ - 2 \sigma_2^2 + \sigma_3^2 & \sigma_2^2 + \sigma_3^2 & \sigma_2^2 - 2 \sigma_3^2\\ - 2 \sigma_2^2 - 2 \sigma_3^2 & \sigma_2^2 - 2 \sigma_3^2 & \sigma_2^2 + 4 \sigma_3^2\\ \end{bmatrix}\].
The decomposition shows that for both the two-level and the three-level case the maxLMM m1 and the zcpLMM m3 differ by a component that assembles all the correlation parameters. Thus, the zcpLMM is nested under the maxLMM in both cases.
3 General forms for sum contrast
Here are the \(\Sigma\)s for sum contrasts for factors with two and three levels.
\(\Sigma\) for LMM m1, k=2
\[\mathbf\Sigma = \begin{bmatrix} \sigma_1^2 + \sigma_2^2 + 2\sigma_{12} & \sigma_1^2 - \sigma_2^2 \\ \sigma_1^2 - \sigma_2^2 & \sigma_1^2 + \sigma_2^2 - 2\sigma_{12} \\ \end{bmatrix} = \begin{bmatrix} \sigma_1^2 & \sigma_1^2 \\ \sigma_1^2 & \sigma_1^2 \\ \end{bmatrix} + \begin{bmatrix} \sigma_2^2 & -\sigma_2^2 \\ - \sigma_2^2 & \sigma_2^2 \\ \end{bmatrix} + \begin{bmatrix} 2\sigma_{12} & 0 \\ 0 & - 2\sigma_{12} \\ \end{bmatrix}\]
\(\Sigma\) for LMM m3, k=2
\[\mathbf\Sigma = \begin{bmatrix} \sigma_1^2 & \sigma_1^2 \\ \sigma_1^2 & \sigma_1^2 \\ \end{bmatrix} + \begin{bmatrix} \sigma_2^2 & -\sigma_2^2 \\ -\sigma_2^2 & \sigma_2^2 \\ \end{bmatrix}\]
\(\Sigma\) for LMM m1, k=3
\[\mathbf\Sigma = \begin{bmatrix} \sigma_1^2 + \sigma_2^2 + 2\sigma_{12} & \sigma_1^2 + \sigma_{12}+ \sigma_{13} + \sigma_{23} & \sigma_1^2 - \sigma_2^2 - \sigma_{13} - \sigma_{23} \\ \sigma_1^2 + \sigma_{12}+ \sigma_{13} + \sigma_{23} & \sigma_1^2 + \sigma_3^2 + 2\sigma_{13} & \sigma_1^2 - \sigma_3^2 - \sigma_{12} - \sigma_{23} \\ \sigma_1^2 - \sigma_2^2 - \sigma_{13} - \sigma_{23} & \sigma_1^2 - \sigma_3^2 - \sigma_{12} - \sigma_{23} & \sigma_1^2 + \sigma_2^2 + \sigma_3^2 - 2\sigma_{12} - 2\sigma_{13} + 2\sigma_{23} \\ \end{bmatrix}\]
This matrix corresponds to the following decomposition:
\[\mathbf \Sigma = \begin{bmatrix} \sigma_1^2 & \sigma_1^2 & \sigma_1^2 \\ \sigma_1^2 & \sigma_1^2 & \sigma_1^2 \\ \sigma_1^2 & \sigma_1^2 & \sigma_1^2 \\ \end{bmatrix} + \begin{bmatrix} \sigma_2^2 & 0 & - \sigma_2^2\\ 0 & \sigma_3^2 & - \sigma_3^2\\ - \sigma_2^2 & - \sigma_3^2 & \sigma_2^2 + \sigma_3^2\\ \end{bmatrix} + \begin{bmatrix} 2\sigma_{12} & \sigma_{12}+ \sigma_{13} + \sigma_{23} & - \sigma_{13} - \sigma_{23}\\ \sigma_{12}+ \sigma_{13} + \sigma_{23} & 2\sigma_{13} & - \sigma_{12} - \sigma_{23}\\ - \sigma_{13} - \sigma_{23} & - \sigma_{12} - \sigma_{23} & - 2\sigma_{12} - 2\sigma_{13} + 2\sigma_{23}\\ \end{bmatrix}\]
\(\Sigma\) for LMM m3, k=3
\[\mathbf \Sigma = \begin{bmatrix} \sigma_1^2 & \sigma_1^2 & \sigma_1^2 \\ \sigma_1^2 & \sigma_1^2 & \sigma_1^2 \\ \sigma_1^2 & \sigma_1^2 & \sigma_1^2 \\ \end{bmatrix} + \begin{bmatrix} \sigma_2^2 & 0 & - \sigma_2^2\\ 0 & \sigma_3^2 & - \sigma_3^2\\ - \sigma_2^2 & - \sigma_3^2 & \sigma_2^2 + \sigma_3^2\\ \end{bmatrix}\]
The decomposition shows that for both the two-level and the three-level case the maxLMM m1 and the zcpLMM m3 differ by a component that assembles all the correlation parameters. Thus, the zcpLMM is nested under the maxLMM in both cases.
4 A few co(variance) rules
V, W, X, Y, Z are random variables; a, b, c, d are constants.
var(aX) = cov(aX, aX) = a^2*cov(X,X) = a^2*var(X),
var (X + Y) = var(X) + var(Y) + 2*cov(X,Y)
var (X - Y) = var(X) + var(Y) - 2*cov(X,Y)
var (X + Y + Z) = var(X) + var(Y) + var(Z) + 2*cov(X,Y) + 2*cov(X,Z) + 2*cov(Y,Z)
cov(aX + bY, cV + dW) = ac*cov(X,V) + ad*cov(X,W) + bc*cov(Y,V) + bd*cov(Y,W)5 Illustrations for two levels (k = 2)
5.1 LMM m1a
# LMM x0
x0 <- lmer(score ~ 0 + Machine + (0 + Machine | Worker), data=Machines2, REML=FALSE)
(vc0 <- as.data.frame(VarCorr(x0)))## grp var1 var2 vcov sdcor
## 1 Worker MachineA <NA> 13.737191 3.7063717
## 2 Worker MachineB <NA> 61.866449 7.8655228
## 3 Worker MachineA MachineB 23.537283 0.8073833
## 4 Residual <NA> <NA> 1.160278 1.0771619A <- vc0[1, "vcov"]
B <- vc0[2, "vcov"]
AB <- vc0[3, "vcov"]5.2 Sliding contrast
var(GM) = var(1/2*A + 1/2*B) = 1/4*var(A+B) = 1/4*(A^2 + B^2 + 2*AB)
var(dBA) = var(B - A) = A^2 + B^2 - 2*AB
cov(GM, dBA) = 1/2*cov(A + B, B - A) = 1/2*(B^2 -A^2 + AB - AB) = 1/2*( B^2 - A^2) contrasts(Machines2$Machine) <- MASS::contr.sdif(2)
contrasts(Machines2$Machine)## 2-1
## A -0.5
## B 0.5mm2 <- model.matrix(~ 1 + Machine, Machines2)
dBA <- mm2[, 2]
var_GM_1 <- 1/4 *(A + B + 2*AB)
var_dBA_1 <- A + B - 2*AB
cov_GM_dBA_1 <- 1/2*(B - A)
cat("var_GM: ", var_GM_1, "var_dBA: ", var_dBA_1, "cov_GM_dBA: ", cov_GM_dBA_1)## var_GM: 30.66955 var_dBA: 28.52907 cov_GM_dBA: 24.06463x1 <- lmer(score ~ 1 + dBA + (1 + dBA | Worker), data=Machines2, REML=FALSE)
(vc1 <- data.frame(VarCorr(x1)))## grp var1 var2 vcov sdcor
## 1 Worker (Intercept) <NA> 30.669553 5.5380099
## 2 Worker dBA <NA> 28.529073 5.3412614
## 3 Worker (Intercept) dBA 24.064629 0.8135451
## 4 Residual <NA> <NA> 1.160278 1.0771619(sigmas_x1 <- setNames(c(vc1[1, "vcov"] + 1/4*vc1[2, "vcov"] - vc1[3, "vcov"],
vc1[1, "vcov"] + 1/4*vc1[2, "vcov"] + vc1[3, "vcov"],
vc1[1, "vcov"] - 1/4*vc1[2, "vcov"]),
c("sigma_1^2", "sigma_2^2", "sigma_12")))## sigma_1^2 sigma_2^2 sigma_12
## 13.73719 61.86645 23.53728V1 <- getME_V(x1)
round(V1[seq(1,6,3), seq(1,6,3)], 2)## 2 x 2 sparse Matrix of class "dgCMatrix"
## 1 19
## 1 13.74 23.54
## 19 23.54 61.87# save
x1_sld_2 <- x1
V1_sld_2 <- V1
# LMM x3
x3 <- lmer(score ~ 1 + dBA + (1 + dBA || Worker), data=Machines2, REML=FALSE)
(vc3 <- data.frame(VarCorr(x3)))## grp var1 var2 vcov sdcor
## 1 Worker (Intercept) <NA> 30.669553 5.538010
## 2 Worker.1 dBA <NA> 28.529074 5.341261
## 3 Residual <NA> <NA> 1.160278 1.077162(sigmas_x3 <- setNames(c(vc3[1, "vcov"] + 1/4*vc3[2, "vcov"],
vc3[1, "vcov"] + 1/4*vc3[2, "vcov"],
vc3[1, "vcov"] - 1/4*vc3[2, "vcov"]),
c("sigma_1^2", "sigma_2^2", "sigma_12")))## sigma_1^2 sigma_2^2 sigma_12
## 37.80182 37.80182 23.53728V3 <- getME_V(x3)
round(V3[seq(1,6,3), seq(1,6,3)], 4)## 2 x 2 sparse Matrix of class "dgCMatrix"
## 1 19
## 1 37.8018 23.5373
## 19 23.5373 37.8018# save
x3_sld_2 <- x3
V3_sld_2 <- V35.3 Sum contrast
var(GM) = var(1/2*A + 1/2*B) = 1/4*var(A + B) = 1/4*(A^2 + B^2 + 2*AB)
var(dA) = var(1/2*A - 1/2*B) = 1/4*var(A - B) = 1/4*(A^2 + B^2 - 2*AB)
cov(GM, dA) = 1/4* cov(A + B, A - B) = 1/4*(A^2 - AB + AB - B^2) = 1/4*(A^2 - B^2)contrasts(Machines2$Machine) <- contr.sum(2)
contrasts(Machines2$Machine)## [,1]
## A 1
## B -1mm <- model.matrix(~ 1 + Machine, Machines2)
dA <- mm[, 2]
var_GM_1 <- 1/4*(A + B + 2*AB)
var_dA_1 <- 1/4*(A + B - 2*AB)
cov_GM_dA_1 <- 1/4*(A - B)
cat("var_GM: ", var_GM_1, "var_dA: ", var_dA_1, "cov_GM_dA: ", cov_GM_dA_1)## var_GM: 30.66955 var_dA: 7.132269 cov_GM_dA: -12.03231x1 <- lmer(score ~ 1 + dA + (1 + dA | Worker), data=Machines2, REML=FALSE)
(vc1 <- as.data.frame(VarCorr(x1)))## grp var1 var2 vcov sdcor
## 1 Worker (Intercept) <NA> 30.669551 5.5380096
## 2 Worker dA <NA> 7.132268 2.6706307
## 3 Worker (Intercept) dA -12.032314 -0.8135451
## 4 Residual <NA> <NA> 1.160278 1.0771619(sigmas_x1 <- setNames(c(sum(vc1[1:2, "vcov"]) + 2*vc1[3, "vcov"],
sum(vc1[1:2, "vcov"]) - 2*vc1[3, "vcov"],
sum(vc1[1, "vcov"]) - vc1[2, "vcov"]),
c("sigma_1^2", "sigma_2^2", "sigma_12")))## sigma_1^2 sigma_2^2 sigma_12
## 13.73719 61.86645 23.53728V1 <- getME_V(x1)
round(V1[seq(1,6,3), seq(1,6,3)], 4)## 2 x 2 sparse Matrix of class "dgCMatrix"
## 1 19
## 1 13.7372 23.5373
## 19 23.5373 61.8664# save
x1_sum_2 <- x1
V1_sum_2 <- V1
# LMM x3
x3 <- lmer(score ~ 1 + dA + (1 + dA || Worker), data=Machines2, REML=FALSE)
(vc3 <- data.frame(VarCorr(x3)))## grp var1 var2 vcov sdcor
## 1 Worker (Intercept) <NA> 30.669553 5.538010
## 2 Worker.1 dA <NA> 7.132268 2.670631
## 3 Residual <NA> <NA> 1.160278 1.077162(sigmas_x3 <- setNames(c(sum(vc3[1:2, "vcov"]), sum(vc3[1:2, "vcov"]),
vc3[1, "vcov"] - vc3[2, "vcov"]),
c("sigma_1^2", "sigma_2^2", "sigma_12")))## sigma_1^2 sigma_2^2 sigma_12
## 37.80182 37.80182 23.53729V3 <- getME_V(x3)
round(V3[seq(1,6,3), seq(1,6,3)], 4)## 2 x 2 sparse Matrix of class "dgCMatrix"
## 1 19
## 1 37.8018 23.5373
## 19 23.5373 37.8018# save
x3_sum_2 <- x3
V3_sum_2 <- V36 Illustration for three levels (k = 3)
6.1 LMM m1a
# LMM x0
x0 <- lmer(score ~ 0 + Machine + (0 + Machine | Worker), data=Machines, REML=FALSE)
(vc0 <- as.data.frame(VarCorr(x0)))## grp var1 var2 vcov sdcor
## 1 Worker MachineA <NA> 13.8157474 3.7169541
## 2 Worker MachineB <NA> 61.9450143 7.8705155
## 3 Worker MachineC <NA> 16.0049079 4.0006134
## 4 Worker MachineA MachineB 23.5372939 0.8045743
## 5 Worker MachineA MachineC 9.2887699 0.6246610
## 6 Worker MachineB MachineC 24.3200687 0.7723871
## 7 Residual <NA> <NA> 0.9246296 0.9615766A <- vc0[1, "vcov"]
B <- vc0[2, "vcov"]
C <- vc0[3, "vcov"]
AB <- vc0[4, "vcov"]
AC <- vc0[5, "vcov"]
BC <- vc0[6, "vcov"]6.2 Sliding contrasts
var(GM) = var(1/9*A + 1/9*B + 1/9*C) = 1/9*var(A+B+C) =
= 1/9*(A^2 + B^2 + C^2 + 2*AB + 2*AC + 2*BC)
var(dBA = var(B - A) = A^2 + B^2 - 2*AB
var(dCB = var(C - B) = B^2 + C^2 - 2*BC
cov(GM, dBA) = 1/3*cov(A+B+C, B-A) = 1/3*(AB - A^2 + B^2 - AB + BC - AC) =
= 1/3*(-A^2 + B^2 + BC - AC)
cov(GM, dCB) = 1/3*cov(A+B+C, C-B) = 1/3*(AC - AB + BC - B^2 + C^2 - BC) =
= 1/3*(-B^2 + C^2 + AC - AB)
cov(dBA, dCB) = cov(B-A, C-B) = BC - B^2 - AC + AB = -B^2 + AB - AC + BCNote. A, B, C are short for A^2, B^2, and C^2.
contrasts(Machines$Machine) <- MASS::contr.sdif(3)
contrasts(Machines$Machine)## 2-1 3-2
## A -0.6666667 -0.3333333
## B 0.3333333 -0.3333333
## C 0.3333333 0.6666667mm <- model.matrix(~ 1 + Machine, Machines)
dBA <- mm[, 2]
dCB <- mm[, 3]
var_GM_1 <- 1/9*(A + B + C + 2*AB + 2*AC + 2*BC)
var_dBA_1 <- A + B - 2*AB
var_dCB_1 <- B + C - 2*BC
cov_GM_dBA_1 <- 1/3*(-A + B + BC - AC)
cov_GM_dCB_1 <- 1/3*(-B + C + AC - AB)
cov_dBA_dCB_1 <- -B + AB - AC + BC
cat("var_GM: ", var_GM_1, "var_dBA: ", var_dBA_1, "var_dCB: ", var_dCB_1, "\n",
"cov_GM_dBA: ", cov_GM_dBA_1, "cov_GM_dCB: ", cov_GM_dCB_1, "cov_dBA_dCB: ", cov_dBA_dCB_1)## var_GM: 22.89533 var_dBA: 28.68617 var_dCB: 29.30978
## cov_GM_dBA: 21.05352 cov_GM_dCB: -20.06288 cov_dBA_dCB: -23.37642# var_GM: 22.89533 var_dBA: 28.68617 var_dCB: 29.30978
# cov_GM_dBA: 21.05352 cov_GM_dCB: -20.06288 cov_dBA_dCB: -23.37642
x1 <- lmer(score ~ 1 + dBA + dCB + (1 + dBA + dCB | Worker), data=Machines, REML=FALSE)
(vc1 <- data.frame(VarCorr(x1)))## grp var1 var2 vcov sdcor
## 1 Worker (Intercept) <NA> 22.8953139 4.7849048
## 2 Worker dBA <NA> 28.6861763 5.3559478
## 3 Worker dCB <NA> 29.3097833 5.4138511
## 4 Worker (Intercept) dBA 21.0535162 0.8215141
## 5 Worker (Intercept) dCB -20.0628670 -0.7744857
## 6 Worker dBA dCB -23.3764229 -0.8061863
## 7 Residual <NA> <NA> 0.9246296 0.9615766(sigmas_x1 <- setNames(c(vc1[1, "vcov"] + 4/9*vc1[2, "vcov"] + 1/9*vc1[3, "vcov"] - 4/3*vc1[4, "vcov"] - 2/3*vc1[5, "vcov"] + 4/9*vc1[6, "vcov"],
vc1[1, "vcov"] + 1/9*vc1[2, "vcov"] + 1/9*vc1[3, "vcov"] + 2/3*vc1[4, "vcov"] - 2/3*vc1[5, "vcov"] - 2/9*vc1[6, "vcov"],
vc1[1, "vcov"] + 1/9*vc1[2, "vcov"] + 4/9*vc1[3, "vcov"] + 2/3*vc1[4, "vcov"] + 4/3*vc1[5, "vcov"] + 4/9*vc1[6, "vcov"],
vc1[1, "vcov"] - 2/9*vc1[2, "vcov"] + 1/9*vc1[3, "vcov"] - 1/3*vc1[4, "vcov"] - 2/3*vc1[5, "vcov"] + 1/9*vc1[6, "vcov"],
vc1[1, "vcov"] - 2/9*vc1[2, "vcov"] - 2/9*vc1[3, "vcov"] - 1/3*vc1[4, "vcov"] + 1/3*vc1[5, "vcov"] - 5/9*vc1[6, "vcov"],
vc1[1, "vcov"] + 1/9*vc1[2, "vcov"] - 2/9*vc1[3, "vcov"] + 2/3*vc1[4, "vcov"] + 1/3*vc1[5, "vcov"] + 1/9*vc1[6, "vcov"]),
c("sigma_1^2", "sigma_2^2", "sigma_3^2",
"sigma_12", "sigma_13", "sigma_23")))## sigma_1^2 sigma_2^2 sigma_3^2 sigma_12 sigma_13 sigma_23
## 13.815737 61.944992 16.004904 23.537276 9.288763 24.320056V1 <- getME_V(x1)
round(V1[seq(1,9,3), seq(1,9,3)], 4)## 3 x 3 sparse Matrix of class "dgCMatrix"
## 1 19 39
## 1 13.8157 23.5373 9.2888
## 19 23.5373 61.9450 24.3201
## 39 9.2888 24.3201 16.0049# save
x1_sld_3 <- x1
V1_sld_3 <- V1
# LMM x3
x3 <- lmer(score ~ 1 + dBA + dCB + (1 + dBA + dCB || Worker), data=Machines, REML=FALSE)
(vc3 <- data.frame(VarCorr(x3)))## grp var1 var2 vcov sdcor
## 1 Worker (Intercept) <NA> 22.8950258 4.7848747
## 2 Worker.1 dBA <NA> 28.1933836 5.3097442
## 3 Worker.2 dCB <NA> 28.8065429 5.3671727
## 4 Residual <NA> <NA> 0.9272635 0.9629452(sigmas_x3 <- setNames(c(vc3[1, "vcov"] + 4/9*vc3[2, "vcov"] + 1/9*vc3[3, "vcov"],
vc3[1, "vcov"] + 1/9*vc3[2, "vcov"] + 1/9*vc3[3, "vcov"],
vc3[1, "vcov"] + 1/9*vc3[2, "vcov"] + 4/9*vc3[3, "vcov"],
vc3[1, "vcov"] - 2/9*vc3[2, "vcov"] + 1/9*vc3[3, "vcov"],
vc3[1, "vcov"] - 2/9*vc3[2, "vcov"] - 2/9*vc3[3, "vcov"],
vc3[1, "vcov"] + 1/9*vc3[2, "vcov"] - 2/9*vc3[3, "vcov"]),
c("sigma_1^2", "sigma_2^2", "sigma_3^2",
"sigma_12", "sigma_13", "sigma_23")))## sigma_1^2 sigma_2^2 sigma_3^2 sigma_12 sigma_13 sigma_23
## 38.62615 29.22835 38.83053 19.83056 10.22838 19.62617V3 <- getME_V(x3)
round(V3[seq(1,9,3), seq(1,9,3)], 4)## 3 x 3 sparse Matrix of class "dgCMatrix"
## 1 19 39
## 1 38.6261 19.8306 10.2284
## 19 19.8306 29.2284 19.6262
## 39 10.2284 19.6262 38.8305# save
x3_sld_3 <- x3
V3_sld_3 <- V3
d_m1_m3 <- V1[seq(1,9,3), seq(1,9,3)] - V3[seq(1,9,3), seq(1,9,3)]6.3 Sum contrasts
var(GM) = var(1/3*A + 1/3*B + 1/3*C) = 1/3*var(A+B+C) =
= 1/9*(A^2 + B^2 + C^2 + 2*AB + 2*AC + 2*BC)
cov(dA, dA) = cov(2/3*A - 1/3*B - 1/3*C, 2/3*A - 1/3*B - 1/3*C) =
= 4/9*A^2 - 2/9*AB - 2/9*AC - 2/9*AB + 1/9*B^2 + 1/9*BC - 2/9*AC + 1/9*BC + 1/9*C^2) =
= 1/9*(4*A^2 + B^2 + C^2 - 4*AB - 4*AC + 2*BC)
cov(dB, dB) = cov(2/3*B - 1/3*A - 1/3*C, 2/3*B - 1/3*A - 1/3*C) =
= 4/9*B^2 - 2/9*AB - 2/9*BC - 2/9*AB + 1/9*A^2 + 1/9*AC - 2/9*BC + 1/9*AC + 1/9*C^2) =
= 1/9*(4*B^2 + A^2 + C^2 - 4*AB - 4*BC + 2*AC)
cov(GM, dA) = cov(1/3*A + 1/3*B + 1/3*C, 2/3*A - 1/3*B - 1/3*C) =
= 2/9*A^2 - 1/9*B^2 - 1/9*C^2 + 1/9*AB + 1/9*AC - 2/9*BC =
= 1/9*(2*A^2 - B^2 - C^2 + AB + AC - 2*BC)
cov(GM, dB) = cov(1/3*B + 1/3*A + 1/3*C, 2/3*B - 1/3*A - 1/3*C) =
= 2/9*B^2 - 1/9*A^2 - 1/9*C^2 + 1/9*AB + 1/9*BC - 2/9*AC =
= 1/9*(2*B^2 - A^2 - C^2 + AB + BC - 2*AC)
cov(dA, dB) = cov(2/3*A - 1/3*B - 1/3*C, 2/3*B - 1/3*A - 1/3*C) =
= 4/9*AB - 2/9*A^2 - 2/9*AC - 2/9*B^2 + 1/9*AB + 1/9*BC - 2/9*BC + 1/9*AC + 1/9*C^2 =
= -2/9*A^2 - 2/9*B^2 + 1/9*C^2 + 5/9*AB -1/9*AC - 1/9*BC =
= 1/9*(-2*A^2 - 2*B^2 + C^2 + 5*AB - AC - BC)contrasts(Machines$Machine) <- contr.sum(3)
contrasts(Machines$Machine)## [,1] [,2]
## A 1 0
## B 0 1
## C -1 -1mm <- model.matrix(~ 1 + Machine, Machines)
dA <- mm[, 2]
dB <- mm[, 3]
var_GM_1 <- 1/9*(A + B + C + 2*AB + 2*AC + 2*BC)
var_dA_1 <- 1/9*(4*A + B + C - 4*AB - 4*AC + 2*BC)
var_dB_1 <- 1/9*(4*B + A + C - 4*AB - 4*BC + 2*AC)
cov_GM_dA_1 <- 1/9*(2*A - B - C + AB + AC - 2*BC)
cov_GM_dB_1 <- 1/9*(2*B - A - C + AB + BC - 2*AC)
cov_dA_dB_1 <- 1/9*(-2*A - 2*B + C + 5*AB - AC - BC)
cat("var_GM: ", var_GM_1, "var_dA: ", var_dA_1, "var_dB: ", var_dB_1, "\n",
"cov_GM_dA: ", cov_GM_dA_1, "cov_GM_dB: ", cov_GM_dB_1, "cov_dA_dB: ", cov_dA_dB_1)## var_GM: 22.89533 var_dA: 5.616533 var_dB: 11.63876
## cov_GM_dA: -7.348056 cov_GM_dB: 13.70547 cov_dA_dB: -5.715443# var_GM: 22.89533 var_dA: 5.61653 var_dB: 11.63876
# cov_GM_dA: -7.34806 cov_GM_dB: 13.70547 cov_dA_dB: -5.71544
x1 <- lmer(score ~ 1 + dA + dB + (1 + dA + dB | Worker), data=Machines, REML=FALSE)
(vc1 <- data.frame(VarCorr(x1)))## grp var1 var2 vcov sdcor
## 1 Worker (Intercept) <NA> 22.8953063 4.7849040
## 2 Worker dA <NA> 5.6165342 2.3699228
## 3 Worker dB <NA> 11.6387526 3.4115616
## 4 Worker (Intercept) dA -7.3480539 -0.6479849
## 5 Worker (Intercept) dB 13.7054559 0.8395896
## 6 Worker dA dB -5.7154435 -0.7069074
## 7 Residual <NA> <NA> 0.9246297 0.9615767(sigmas_x1 <- setNames(c(sum(vc1[1:2, "vcov"]) + 2*vc1[4, "vcov"],
sum(vc1[c(1,3), "vcov"]) + 2*vc1[5, "vcov"],
sum(vc1[1:3, "vcov"]) - 2*vc1[4, "vcov"] - 2*vc1[5, "vcov"] + 2*vc1[6, "vcov"],
vc1[1, "vcov"] + sum(vc1[4:6, "vcov"]),
vc1[1, "vcov"] - vc1[2, "vcov"] - vc1[5, "vcov"] - vc1[6, "vcov"],
vc1[1, "vcov"] - vc1[3, "vcov"] - vc1[4, "vcov"] - vc1[6, "vcov"]),
c("sigma_1^2", "sigma_2^2", "sigma_3^2",
"sigma_12", "sigma_13", "sigma_23")))## sigma_1^2 sigma_2^2 sigma_3^2 sigma_12 sigma_13 sigma_23
## 13.81573 61.94497 16.00490 23.53726 9.28876 24.32005V1 <- getME_V(x1)
round(V1[seq(1,9,3), seq(1,9,3)], 2)## 3 x 3 sparse Matrix of class "dgCMatrix"
## 1 19 39
## 1 13.82 23.54 9.29
## 19 23.54 61.94 24.32
## 39 9.29 24.32 16.00# save
x1_sum_3 <- x1
V1_sum_3 <- V1
# LMM x3
x3 <- lmer(score ~ 1 + dA + dB + (1 + dA + dB || Worker), data=Machines, REML=FALSE)
(vc3 <- data.frame(VarCorr(x3)))## grp var1 var2 vcov sdcor
## 1 Worker (Intercept) <NA> 22.8950150 4.784874
## 2 Worker.1 dA <NA> 5.5147425 2.348349
## 3 Worker.2 dB <NA> 11.4323044 3.381169
## 4 Residual <NA> <NA> 0.9273633 0.962997(sigmas3 <- setNames(c(sum(vc3[c(1,2), "vcov"]), sum(vc3[c(1,3), "vcov"]), sum(vc3[1:3, "vcov"]),
vc3[1, "vcov"], vc3[1, "vcov"] - vc3[2, "vcov"], vc3[1, "vcov"] - vc3[3, "vcov"]),
c("sigma_1^2", "sigma_2^2", "sigma_3^2", "sigma_12", "sigma_13", "sigma_23"))) ## sigma_1^2 sigma_2^2 sigma_3^2 sigma_12 sigma_13 sigma_23
## 28.40976 34.32732 39.84206 22.89501 17.38027 11.46271V3 <- getME_V(x3)
round(V3[seq(1,9,3), seq(1,9,3)], 2)## 3 x 3 sparse Matrix of class "dgCMatrix"
## 1 19 39
## 1 28.41 22.90 17.38
## 19 22.90 34.33 11.46
## 39 17.38 11.46 39.84# save
x3_sum_3 <- x3
V3_sum_3 <- V37 Model comparisons
7.1 Sliding contrasts
suppressMessages(tab_model(x1_sld_2, x3_sld_2, x1_sld_3, x3_sld_3,
dv.labels = c("max_sld_2", "zcp_sld_2", "max_sld_3", "zcp_sld_3"), show.aic=TRUE,
show.ci=NULL, show.p=FALSE, show.r2=FALSE, show.icc=FALSE, show.dev=TRUE))Models with 2 factor levels differ by 1 df; models with 3 factor levels differ by 3 df.
7.2 Sum contrasts
suppressMessages(tab_model(x1_sum_2, x3_sum_2, x1_sum_3, x3_sum_3,
dv.labels = c("max_sum_2", "zcp_sum_2", "max_sum_3", "zcp_sum_3"), show.aic=TRUE,
show.ci=NULL, show.p=FALSE, show.r2=FALSE, show.icc=FALSE, show.dev=TRUE))| max_sum_2 | zcp_sum_2 | max_sum_3 | zcp_sum_3 | |
|---|---|---|---|---|
| Predictors | Estimates | Estimates | Estimates | Estimates |
| (Intercept) | 56.34 | 56.34 | 59.65 | 59.65 |
| d A | -3.98 | -3.98 | -7.29 | -7.29 |
| d B | 0.67 | 0.67 | ||
| Random Effects | ||||
| σ2 | 1.16 | 1.16 | 0.92 | 0.93 |
| τ00 | 30.67 Worker | 30.67 Worker | 22.90 Worker | 22.90 Worker |
| 7.13 Worker.1 | 5.51 Worker.1 | |||
| 11.43 Worker.2 | ||||
| τ11 | 7.13 Worker.dA | 5.62 Worker.dA | ||
| 11.64 Worker.dB | ||||
| ρ01 | -0.81 Worker | -0.65 Worker.dA | ||
| 0.84 Worker.dB | ||||
| Observations | 36 | 36 | 54 | 54 |
| Deviance | 153.622 | 159.758 | 216.418 | 227.491 |
| AIC | 165.622 | 169.758 | 236.418 | 241.491 |
Models with 2 factor levels differ by 1 df; models with 3 factor levels differ by 3 df.
8 Images of random-effect structures
The hiearchical structure of LMMs m1 to m4 can also be seen in comparison of visualization based on estimates for the variance-covariance matrix V based on 3 machines x 6 workers x 3 measures (i.e., 54 x 54).
8.1 Sliding contrasts
Vs <- list(V1_sld_2, V3_sld_2, V1_sld_3, V3_sld_3)
for(i in 1:4) print(image(Vs[[i]], asp=1))
V1_sld_2[seq(1,6,3), seq(1,6,3)]## 2 x 2 sparse Matrix of class "dgCMatrix"
## 1 19
## 1 13.73719 23.53728
## 19 23.53728 61.86645V3_sld_2[seq(1,6,3), seq(1,6,3)]## 2 x 2 sparse Matrix of class "dgCMatrix"
## 1 19
## 1 37.80182 23.53728
## 19 23.53728 37.80182V1_sld_3[seq(1,9,3), seq(1,9,3)]## 3 x 3 sparse Matrix of class "dgCMatrix"
## 1 19 39
## 1 13.815737 23.53728 9.288763
## 19 23.537276 61.94499 24.320056
## 39 9.288763 24.32006 16.004904V3_sld_3[seq(1,9,3), seq(1,9,3)]## 3 x 3 sparse Matrix of class "dgCMatrix"
## 1 19 39
## 1 38.62615 19.83056 10.22838
## 19 19.83056 29.22835 19.62617
## 39 10.22838 19.62617 38.83053The 3-level case shows that at the correlations between observations within groups are not compatible with compound symmetry. This is the case for both the maxLMM m1 and the zcpLMM m3.
8.2 Sum contrasts
Vs <- list(V1_sum_2, V3_sum_2, V1_sum_3, V3_sum_3)
for(i in 1:4) print(image(Vs[[i]], asp=1))
V1_sum_2[seq(1,6,3), seq(1,6,3)]## 2 x 2 sparse Matrix of class "dgCMatrix"
## 1 19
## 1 13.73719 23.53728
## 19 23.53728 61.86645V3_sum_2[seq(1,6,3), seq(1,6,3)]## 2 x 2 sparse Matrix of class "dgCMatrix"
## 1 19
## 1 37.80182 23.53729
## 19 23.53729 37.80182V1_sum_3[seq(1,9,3), seq(1,9,3)]## 3 x 3 sparse Matrix of class "dgCMatrix"
## 1 19 39
## 1 13.81573 23.53726 9.28876
## 19 23.53726 61.94497 24.32005
## 39 9.28876 24.32005 16.00490V3_sum_3[seq(1,9,3), seq(1,9,3)]## 3 x 3 sparse Matrix of class "dgCMatrix"
## 1 19 39
## 1 28.40976 22.89501 17.38027
## 19 22.89501 34.32732 11.46271
## 39 17.38027 11.46271 39.842069 Conclusion
This exercise demonstrates that LMM m3 is nested under LMM m1 (as it was the case for treatment contrasts (Kliegl, 2018f). Thus, LRTs can be used for model comparison. The elements of \(\mathbf\Sigma\) correspond in a transparent way with the model specification in the lmer syntax. The demonstrations also show that at the level of observations within groups (workers) compound symmetry does not hold for 3-level factors, neither for sliding, nor for sum contrasts. Thus, the variance and, in the case of the maxLMM, covariance parameters at the group level may combine to different within-group correlations (i.e., we may end up with six different values.) For treatment contrasts in the zcpLMM, only variance parameters differ; all correlation parameters are equal to the variance parameter of the reference category (see Kliegl, 2018f).
10 References
- Bates, D., Kliegl, R., Vasishth, S., & Baayen, H. (2015a). Parsimonious mixed models. arXiv:1506.04967v2 stat.ME
- Bates, D., Mächler, M., Bolker, B.M., & Walker, S.C. (2015b). Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, 67, 1
- Matuschek et al. (2017). Balancing Type I error and power in linear mixed models. Journal of Memory and Language, 94, 305-315.
Related RPubs:
- Kliegl, R. (2018a). Reduction of complexity of linear mixed models with double-bar syntax (Version 2)
- Kliegl, R. (2018b). Factor contrasts and factor levels in linear mixed models
- Kliegl, R. (2018c). Hierarchical random-effect structures in linear mixed models
- Kliegl, R. (2018d). Images of Contrast-Dependent Random-Effect Structures in Linear Mixed Models
- Kliegl, R. (2018e). Linear Mixed Models with Contrast-Based Measures
- Kliegl, R. (2018f). Decompositon of a Linear Mixed Model with Treatment Contrasts
11 Acknowledgment
This script was motivated and inspired by exchanges with Rune H. B. Christensen. I gratefully acknowledge his insights and comments.
12 Appendix
sessionInfo()## R version 3.5.1 (2018-07-02)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS 10.14
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] bindrcpp_0.2.2 modelr_0.1.2 forcats_0.3.0 stringr_1.3.1
## [5] dplyr_0.7.6 purrr_0.2.5 readr_1.1.1 tidyr_0.8.1
## [9] tibble_1.4.2 ggplot2_3.0.0 tidyverse_1.2.1 sjPlot_2.6.0
## [13] lme4_1.1-18-1 Matrix_1.2-14
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-137 lubridate_1.7.4 httr_1.3.1
## [4] tools_3.5.1 TMB_1.7.14 backports_1.1.2
## [7] R6_2.2.2 sjlabelled_1.0.14 lazyeval_0.2.1
## [10] colorspace_1.3-2 nnet_7.3-12 withr_2.1.2
## [13] tidyselect_0.2.4 mnormt_1.5-5 emmeans_1.2.4
## [16] compiler_3.5.1 cli_1.0.1 rvest_0.3.2
## [19] xml2_1.2.0 sandwich_2.5-0 effects_4.0-3
## [22] scales_1.0.0 mvtnorm_1.0-8 psych_1.8.4
## [25] ggridges_0.5.0 digest_0.6.17 foreign_0.8-71
## [28] minqa_1.2.4 rmarkdown_1.10.11 stringdist_0.9.5.1
## [31] pkgconfig_2.0.2 htmltools_0.3.6 pwr_1.2-2
## [34] rlang_0.2.2 readxl_1.1.0 htmldeps_0.1.1
## [37] rstudioapi_0.7 bindr_0.1.1 zoo_1.8-4
## [40] jsonlite_1.5 magrittr_1.5 modeltools_0.2-22
## [43] bayesplot_1.6.0 Rcpp_0.12.18 munsell_0.5.0
## [46] prediction_0.3.6 stringi_1.2.4 multcomp_1.4-8
## [49] yaml_2.2.0 snakecase_0.9.2 carData_3.0-1
## [52] MASS_7.3-50 plyr_1.8.4 grid_3.5.1
## [55] parallel_3.5.1 sjmisc_2.7.5 crayon_1.3.4
## [58] lattice_0.20-35 ggeffects_0.5.0 haven_1.1.2
## [61] splines_3.5.1 sjstats_0.17.0 hms_0.4.2
## [64] knitr_1.20 pillar_1.3.0.9000 estimability_1.3
## [67] codetools_0.2-15 stats4_3.5.1 glue_1.3.0
## [70] evaluate_0.11 data.table_1.11.4 nloptr_1.0.4
## [73] cellranger_1.1.0 gtable_0.2.0 assertthat_0.2.0
## [76] coin_1.2-2 xtable_1.8-3 broom_0.5.0
## [79] survey_3.33-2 coda_0.19-1 survival_2.42-6
## [82] glmmTMB_0.2.2.0 TH.data_1.0-9