1 Overview

The zero-correlation parameter linear mixed model (zcpLMM) serves as a reference point for various model comparisons of experimental data. It is of intermediate model complexity because, relative to maximal LMM (maxLMM) it estimates no correlation parameters and relative to the minimal LMM (minLMM) it retains all the variance parameters of the maxLMM. Given its role as a pivot model, it is important to have a clear understanding of its mathematical foundation.

In an analysis with contr.treatment(), on CrossValidated Rune H. B. Christensen and Dmitry Kobak (aka amoeba) in response to a question by Maarten Jung (aka statmerkur, 2018) observed that in zcpLMM the special structure on \(\mathbf\Sigma\) contains a constant covariance so the specification does not correspond to uncorrelated random effects, counter to the reported random-effect structure (RES) which lists no estimates of correlation parameters (CPs), but only estimates of variance parameters (VPs) that correspond directly to orthogonal variance components (VCs) in the model formula, triggered, for example, by the || specification. Here I check whether this result holds also for other contrasts (e.g., contr.sum, MASS::contr.sdif).

2 Function: `getME_V()`

V <- getME_V <- function(model=model){
  Z <- getME(model, "Z")
  Lambda <- getME(model, "Lambda") 
  V <- sigma(model)^2 * Z %*% Lambda %*% t(Lambda) %*% t(Z)
  V
}

Returns variance-covariance matrix V for observations. Residual variance which would add \(\sigma^2 \mathbf I\) to V is ignored.

3 Level-based random-effect structures

3.1 Maximal model `m1a`

Machine data are from six workers (random factor) who operated three machines (fixed factor) each for three times; total number of observations is 6 x 3 x 3 = 54.

# Indicator covariates for RE0
mm0 <- model.matrix(~ 0 + Machine, Machines)
A <- mm0[, 1]
B <- mm0[, 2]
C <- mm0[, 3]

# m1a
m1a <- lmer(score ~ 1 + Machine + (0 + Machine | Worker), data=Machines, REML=FALSE)
VarCorr(m1a)

##  Groups   Name     Std.Dev. Corr       
##  Worker   MachineA 3.71695             
##           MachineB 7.87051  0.805      
##           MachineC 4.00061  0.625 0.772
##  Residual          0.96158

V1a <- getME_V(m1a)
(V1a_W1 <- round(V1a[1:9, 1:9], 2))

## 9 x 9 sparse Matrix of class "dgCMatrix"
##        1     2     3    19    20    21    39    38    37
## 1  13.82 13.82 13.82 23.54 23.54 23.54  9.29  9.29  9.29
## 2  13.82 13.82 13.82 23.54 23.54 23.54  9.29  9.29  9.29
## 3  13.82 13.82 13.82 23.54 23.54 23.54  9.29  9.29  9.29
## 19 23.54 23.54 23.54 61.94 61.94 61.94 24.32 24.32 24.32
## 20 23.54 23.54 23.54 61.94 61.94 61.94 24.32 24.32 24.32
## 21 23.54 23.54 23.54 61.94 61.94 61.94 24.32 24.32 24.32
## 39  9.29  9.29  9.29 24.32 24.32 24.32 16.00 16.00 16.00
## 38  9.29  9.29  9.29 24.32 24.32 24.32 16.00 16.00 16.00
## 37  9.29  9.29  9.29 24.32 24.32 24.32 16.00 16.00 16.00

# variance-covariance for observations from first worker
image(V1a)

For the structure of V of m1a the observations on different workers are independent, but the (co)variances of observations within a worker (across the three machines) are:

equal for the (co)variances from the same machine
different for the (co)variances between different machines

3.2 Models `m2a`, `m3a`, and `m4a`

# chunk-option: out.width='33%' reduces image size, but keeps them in row
## m2a
m2a <- lmer(score ~ 1 + Machine + (0 + A | Worker) + (0 + B + C | Worker), Machines, REML=FALSE)  
V2a <- getME_V(m2a)
(V2a_W1 <- round(V2a[1:9, 1:9], 2))

## 9 x 9 sparse Matrix of class "dgCMatrix"
##        1     2     3    19    20    21    39    38    37
## 1  13.82 13.82 13.82  .     .     .     .     .     .   
## 2  13.82 13.82 13.82  .     .     .     .     .     .   
## 3  13.82 13.82 13.82  .     .     .     .     .     .   
## 19  .     .     .    61.94 61.94 61.94 24.32 24.32 24.32
## 20  .     .     .    61.94 61.94 61.94 24.32 24.32 24.32
## 21  .     .     .    61.94 61.94 61.94 24.32 24.32 24.32
## 39  .     .     .    24.32 24.32 24.32 16.00 16.00 16.00
## 38  .     .     .    24.32 24.32 24.32 16.00 16.00 16.00
## 37  .     .     .    24.32 24.32 24.32 16.00 16.00 16.00

# variance-covariance for observations from first worker

## m3a
m3a <- lmer(score ~ 1 + Machine + (0 + A + B + C || Worker), Machines, REML=FALSE) 
V3a <- getME_V(m3a)
round(V3a[1:9, 1:9], 2)

## 9 x 9 sparse Matrix of class "dgCMatrix"
##        1     2     3    19    20    21 39 38 37
## 1  13.82 13.82 13.82  .     .     .     .  .  .
## 2  13.82 13.82 13.82  .     .     .     .  .  .
## 3  13.82 13.82 13.82  .     .     .     .  .  .
## 19  .     .     .    61.95 61.95 61.95  .  .  .
## 20  .     .     .    61.95 61.95 61.95  .  .  .
## 21  .     .     .    61.95 61.95 61.95  .  .  .
## 39  .     .     .     .     .     .    16 16 16
## 38  .     .     .     .     .     .    16 16 16
## 37  .     .     .     .     .     .    16 16 16

## m4a
m4a <- lmer(score ~ 1 + Machine + (0 + B + C || Worker), Machines, REML=FALSE)
V4a <- getME_V(m4a)
round(V4a[1:9, 1:9], 2)

## 9 x 9 sparse Matrix of class "dgCMatrix"
##    1 2 3    19    20    21    39    38    37
## 1  . . .  .     .     .     .     .     .   
## 2  . . .  .     .     .     .     .     .   
## 3  . . .  .     .     .     .     .     .   
## 19 . . . 59.97 59.97 59.97  .     .     .   
## 20 . . . 59.97 59.97 59.97  .     .     .   
## 21 . . . 59.97 59.97 59.97  .     .     .   
## 39 . . .  .     .     .    14.03 14.03 14.03
## 38 . . .  .     .     .    14.03 14.03 14.03
## 37 . . .  .     .     .    14.03 14.03 14.03

## m5
m5 <- lmer(score ~ 1 + Machine + (1 | Worker), data=Machines, REML=FALSE)
V5<- getME_V(m5)
(V5 <- round(V5[1:9, 1:9],2))

## 9 x 9 sparse Matrix of class "dgCMatrix"
##        1     2     3    19    20    21    39    38    37
## 1  21.93 21.93 21.93 21.93 21.93 21.93 21.93 21.93 21.93
## 2  21.93 21.93 21.93 21.93 21.93 21.93 21.93 21.93 21.93
## 3  21.93 21.93 21.93 21.93 21.93 21.93 21.93 21.93 21.93
## 19 21.93 21.93 21.93 21.93 21.93 21.93 21.93 21.93 21.93
## 20 21.93 21.93 21.93 21.93 21.93 21.93 21.93 21.93 21.93
## 21 21.93 21.93 21.93 21.93 21.93 21.93 21.93 21.93 21.93
## 39 21.93 21.93 21.93 21.93 21.93 21.93 21.93 21.93 21.93
## 38 21.93 21.93 21.93 21.93 21.93 21.93 21.93 21.93 21.93
## 37 21.93 21.93 21.93 21.93 21.93 21.93 21.93 21.93 21.93

Vs <- list(V2a, V3a, V4a)
for(i in 1:3) print(image(Vs[[i]], asp=1))

Model complexity decreases from left to right

For LMM m3a (with a diagonal random-effect variance-covariance matrix), observations for the same Worker on different Machines are independent with covariance zero. We can see this also for the \(9\times 9\) block of observations for the first worker.

This is not the case for LMM m2a due to the additional correlation parameter between machines B and C. Correlations with machine A are zero, as specified in the formula. A comparison between LMM m3a and LMM m2a of the worker’s values suggest that the former is nested under the latter.

LMM m4a has a structure similar to LMM m3a, except that the variance/covariances for machine A, which was not estimated in this model, are zero. The variance/covariance values estimated in both models, however, are different. This was not the case for LMM m2a.

Question: Is LMM m4a nested under LMM m3a?

4 Contrast-based random-effect structures

4.1 Treatment contrast

Now we look at the same models but with a treatment contrast with machine A as base. Of course, the same qualitative pattern of results was obtained with machines B or C as base.

# Extract indicator covariates for treatment contrast with base A
contrasts(Machines$Machine) <- contr.treatment(3, base=1)
mm <- model.matrix(~ 1 + Machine, Machines)
c1 <- mm[, 2]
c2 <- mm[, 3]

# Models

## m1
m1_trt <- lmer(score ~ 1 + Machine + (1 + c1 + c2 | Worker), data=Machines, REML=FALSE)
V1_trt <- getME_V(m1_trt)
(V1_trt_W1 <- round(V1_trt[1:9, 1:9], 2))

## 9 x 9 sparse Matrix of class "dgCMatrix"
##        1     2     3    19    20    21    39    38    37
## 1  13.82 13.82 13.82 23.54 23.54 23.54  9.29  9.29  9.29
## 2  13.82 13.82 13.82 23.54 23.54 23.54  9.29  9.29  9.29
## 3  13.82 13.82 13.82 23.54 23.54 23.54  9.29  9.29  9.29
## 19 23.54 23.54 23.54 61.95 61.95 61.95 24.32 24.32 24.32
## 20 23.54 23.54 23.54 61.95 61.95 61.95 24.32 24.32 24.32
## 21 23.54 23.54 23.54 61.95 61.95 61.95 24.32 24.32 24.32
## 39  9.29  9.29  9.29 24.32 24.32 24.32 16.00 16.00 16.00
## 38  9.29  9.29  9.29 24.32 24.32 24.32 16.00 16.00 16.00
## 37  9.29  9.29  9.29 24.32 24.32 24.32 16.00 16.00 16.00

## m2
m2_trt <- lmer(score ~ 1 + Machine + (1 | Worker) + (0 + c1 + c2 | Worker), data=Machines, REML=FALSE)
V2_trt <- getME_V(m2_trt)
(V2_trt_W1 <- round(V2_trt[1:9, 1:9], 2))

## 9 x 9 sparse Matrix of class "dgCMatrix"
##        1     2     3    19    20    21    39    38    37
## 1  13.72 13.72 13.72 13.72 13.72 13.72 13.72 13.72 13.72
## 2  13.72 13.72 13.72 13.72 13.72 13.72 13.72 13.72 13.72
## 3  13.72 13.72 13.72 13.72 13.72 13.72 13.72 13.72 13.72
## 19 13.72 13.72 13.72 42.83 42.83 42.83 19.14 19.14 19.14
## 20 13.72 13.72 13.72 42.83 42.83 42.83 19.14 19.14 19.14
## 21 13.72 13.72 13.72 42.83 42.83 42.83 19.14 19.14 19.14
## 39 13.72 13.72 13.72 19.14 19.14 19.14 24.76 24.76 24.76
## 38 13.72 13.72 13.72 19.14 19.14 19.14 24.76 24.76 24.76
## 37 13.72 13.72 13.72 19.14 19.14 19.14 24.76 24.76 24.76

## m3
m3_trt <- lmer(score ~ 1 + Machine + (1 + c1 + c2 || Worker), data=Machines, REML=FALSE)
V3_trt <- getME_V(m3_trt)
(V3_trt_W1 <- round(V3_trt[1:9, 1:9], 2))

## 9 x 9 sparse Matrix of class "dgCMatrix"
##        1     2     3    19    20    21    39    38    37
## 1  13.78 13.78 13.78 13.78 13.78 13.78 13.78 13.78 13.78
## 2  13.78 13.78 13.78 13.78 13.78 13.78 13.78 13.78 13.78
## 3  13.78 13.78 13.78 13.78 13.78 13.78 13.78 13.78 13.78
## 19 13.78 13.78 13.78 42.60 42.60 42.60 13.78 13.78 13.78
## 20 13.78 13.78 13.78 42.60 42.60 42.60 13.78 13.78 13.78
## 21 13.78 13.78 13.78 42.60 42.60 42.60 13.78 13.78 13.78
## 39 13.78 13.78 13.78 13.78 13.78 13.78 24.71 24.71 24.71
## 38 13.78 13.78 13.78 13.78 13.78 13.78 24.71 24.71 24.71
## 37 13.78 13.78 13.78 13.78 13.78 13.78 24.71 24.71 24.71

## m4
m4_trt <- lmer(score ~ 1 + Machine + (1 | Worker) + (0 + c2 || Worker), data=Machines, REML=FALSE)
V4_trt <- getME_V(m4_trt)
(V4_trt_W1 <- round(V4_trt[1:9, 1:9], 2))

## 9 x 9 sparse Matrix of class "dgCMatrix"
##        1     2     3    19    20    21    39    38    37
## 1  27.22 27.22 27.22 27.22 27.22 27.22 27.22 27.22 27.22
## 2  27.22 27.22 27.22 27.22 27.22 27.22 27.22 27.22 27.22
## 3  27.22 27.22 27.22 27.22 27.22 27.22 27.22 27.22 27.22
## 19 27.22 27.22 27.22 27.22 27.22 27.22 27.22 27.22 27.22
## 20 27.22 27.22 27.22 27.22 27.22 27.22 27.22 27.22 27.22
## 21 27.22 27.22 27.22 27.22 27.22 27.22 27.22 27.22 27.22
## 39 27.22 27.22 27.22 27.22 27.22 27.22 36.13 36.13 36.13
## 38 27.22 27.22 27.22 27.22 27.22 27.22 36.13 36.13 36.13
## 37 27.22 27.22 27.22 27.22 27.22 27.22 36.13 36.13 36.13

Vs <- list(V2_trt, V3_trt, V4_trt)
for(i in 1:3) print(image(Vs[[i]], asp=1))

Model complexity decreases from left to right

LMM m3 model looks like LMM m1 (which is mathematically identical withm1a): Observations on the same Worker but different Machines have a non-zero covariance in V.

Question: For same structure, different values: Is LMM m3 nested under LMM m1?

The structure of V has the peculiar structure that the covariance of observations on different Machines is the same, and the same as the covariance of observations on the base Machine (in this case: A).

For LMM m4 (without VP c1) the structure of V represents a further simplification. Now all but the covariance of three observations on Machine C are equal. Again, the structure of m4 appears to be simpler at the level of covariances. Is it nested under m3?

For LMM m2 (with CP for c1 and c2) the structure is a bit more complicated relative to m3. At the level of values of covariances it looks like m3 could be nested under m2.

4.2 Sum contrast

I use the default contr.sum().

# Extract indicator covariates for sum contrast
contrasts(Machines$Machine) <- contr.sum(3)
mm <- model.matrix(~ 1 + Machine, Machines)
c1 <- mm[, 2]
c2 <- mm[, 3]

# Models

## m1
m1_sum <- lmer(score ~ 1 + Machine + (1 + c1 + c2 | Worker), data=Machines, REML=FALSE)
V1_sum <- getME_V(m1_sum)
(V1_sum_W1 <- round(V1_sum[1:9, 1:9], 2))

## 9 x 9 sparse Matrix of class "dgCMatrix"
##        1     2     3    19    20    21    39    38    37
## 1  13.82 13.82 13.82 23.54 23.54 23.54  9.29  9.29  9.29
## 2  13.82 13.82 13.82 23.54 23.54 23.54  9.29  9.29  9.29
## 3  13.82 13.82 13.82 23.54 23.54 23.54  9.29  9.29  9.29
## 19 23.54 23.54 23.54 61.94 61.94 61.94 24.32 24.32 24.32
## 20 23.54 23.54 23.54 61.94 61.94 61.94 24.32 24.32 24.32
## 21 23.54 23.54 23.54 61.94 61.94 61.94 24.32 24.32 24.32
## 39  9.29  9.29  9.29 24.32 24.32 24.32 16.00 16.00 16.00
## 38  9.29  9.29  9.29 24.32 24.32 24.32 16.00 16.00 16.00
## 37  9.29  9.29  9.29 24.32 24.32 24.32 16.00 16.00 16.00

## m2
m2_sum <- lmer(score ~ 1 + Machine + (1 | Worker) + (0 + c1 + c2 | Worker), data=Machines, REML=FALSE)
V2_sum <- getME_V(m2_sum)
(V2_sum_W1 <- round(V2_sum[1:9, 1:9], 2))

## 9 x 9 sparse Matrix of class "dgCMatrix"
##        1     2     3    19    20    21    39    38    37
## 1  28.51 28.51 28.51 17.18 17.18 17.18 22.99 22.99 22.99
## 2  28.51 28.51 28.51 17.18 17.18 17.18 22.99 22.99 22.99
## 3  28.51 28.51 28.51 17.18 17.18 17.18 22.99 22.99 22.99
## 19 17.18 17.18 17.18 34.53 34.53 34.53 16.97 16.97 16.97
## 20 17.18 17.18 17.18 34.53 34.53 34.53 16.97 16.97 16.97
## 21 17.18 17.18 17.18 34.53 34.53 34.53 16.97 16.97 16.97
## 39 22.99 22.99 22.99 16.97 16.97 16.97 28.72 28.72 28.72
## 38 22.99 22.99 22.99 16.97 16.97 16.97 28.72 28.72 28.72
## 37 22.99 22.99 22.99 16.97 16.97 16.97 28.72 28.72 28.72

## m3
m3_sum <- lmer(score ~ 1 + Machine + (1 + c1 + c2 || Worker), data=Machines, REML=FALSE)
V3_sum <- getME_V(m3_sum)
(V3_sum_W1 <- round(V3_sum[1:9, 1:9], 2))

## 9 x 9 sparse Matrix of class "dgCMatrix"
##        1     2     3    19    20    21    39    38    37
## 1  28.41 28.41 28.41 22.90 22.90 22.90 17.38 17.38 17.38
## 2  28.41 28.41 28.41 22.90 22.90 22.90 17.38 17.38 17.38
## 3  28.41 28.41 28.41 22.90 22.90 22.90 17.38 17.38 17.38
## 19 22.90 22.90 22.90 34.33 34.33 34.33 11.46 11.46 11.46
## 20 22.90 22.90 22.90 34.33 34.33 34.33 11.46 11.46 11.46
## 21 22.90 22.90 22.90 34.33 34.33 34.33 11.46 11.46 11.46
## 39 17.38 17.38 17.38 11.46 11.46 11.46 39.84 39.84 39.84
## 38 17.38 17.38 17.38 11.46 11.46 11.46 39.84 39.84 39.84
## 37 17.38 17.38 17.38 11.46 11.46 11.46 39.84 39.84 39.84

## m4
m4_sum <- lmer(score ~ 1 + Machine + (1 | Worker) + (0 + c2 || Worker), data=Machines, REML=FALSE)
V4_sum <- getME_V(m4_sum)
(V4_sum_W1 <- round(V4_sum[1:9, 1:9], 2))

## 9 x 9 sparse Matrix of class "dgCMatrix"
##        1     2     3    19    20    21    39    38    37
## 1  22.49 22.49 22.49 22.49 22.49 22.49 22.49 22.49 22.49
## 2  22.49 22.49 22.49 22.49 22.49 22.49 22.49 22.49 22.49
## 3  22.49 22.49 22.49 22.49 22.49 22.49 22.49 22.49 22.49
## 19 22.49 22.49 22.49 29.22 29.22 29.22 15.77 15.77 15.77
## 20 22.49 22.49 22.49 29.22 29.22 29.22 15.77 15.77 15.77
## 21 22.49 22.49 22.49 29.22 29.22 29.22 15.77 15.77 15.77
## 39 22.49 22.49 22.49 15.77 15.77 15.77 29.22 29.22 29.22
## 38 22.49 22.49 22.49 15.77 15.77 15.77 29.22 29.22 29.22
## 37 22.49 22.49 22.49 15.77 15.77 15.77 29.22 29.22 29.22

Vs <- list(V2_sum, V3_sum, V4_sum)
for(i in 1:3) print(image(Vs[[i]], asp=1))

Model complexity decreases from left to right

The pattern for contr.sum is similar to the one for contr.treatment. Covariances are zero between subjects, but not between machines within subjects. The hierarchy between the three models was not as clearly visible as in the last chunk.

4.3 Sliding-difference contrast

I use the default MASS::contr.sdif().

# Extract indicator covariates for sliding difference con
contrasts(Machines$Machine) <- MASS::contr.sdif(3)
mm <- model.matrix(~ 1 + Machine, Machines)
c1 <- mm[, 2]
c2 <- mm[, 3]

# Models

## m1
m1_sld <- lmer(score ~ 1 + Machine + (1 + c1 + c2 | Worker), data=Machines, REML=FALSE)
V1_sld <- getME_V(m1_sld)
(V1_sld_W1 <- round(V1_sld[1:9, 1:9], 2))

## 9 x 9 sparse Matrix of class "dgCMatrix"
##        1     2     3    19    20    21    39    38    37
## 1  13.82 13.82 13.82 23.54 23.54 23.54  9.29  9.29  9.29
## 2  13.82 13.82 13.82 23.54 23.54 23.54  9.29  9.29  9.29
## 3  13.82 13.82 13.82 23.54 23.54 23.54  9.29  9.29  9.29
## 19 23.54 23.54 23.54 61.94 61.94 61.94 24.32 24.32 24.32
## 20 23.54 23.54 23.54 61.94 61.94 61.94 24.32 24.32 24.32
## 21 23.54 23.54 23.54 61.94 61.94 61.94 24.32 24.32 24.32
## 39  9.29  9.29  9.29 24.32 24.32 24.32 16.00 16.00 16.00
## 38  9.29  9.29  9.29 24.32 24.32 24.32 16.00 16.00 16.00
## 37  9.29  9.29  9.29 24.32 24.32 24.32 16.00 16.00 16.00

## m2
m2_sld <- lmer(score ~ 1 + Machine + (1 | Worker) + (0 + c1 + c2 | Worker), data=Machines, REML=FALSE)
V2_sld <- getME_V(m2_sld)
(V2_sld_W1 <- round(V2_sld[1:9, 1:9],2))

## 9 x 9 sparse Matrix of class "dgCMatrix"
##        1     2     3    19    20    21    39    38    37
## 1  28.51 28.51 28.51 17.18 17.18 17.18 22.99 22.99 22.99
## 2  28.51 28.51 28.51 17.18 17.18 17.18 22.99 22.99 22.99
## 3  28.51 28.51 28.51 17.18 17.18 17.18 22.99 22.99 22.99
## 19 17.18 17.18 17.18 34.53 34.53 34.53 16.97 16.97 16.97
## 20 17.18 17.18 17.18 34.53 34.53 34.53 16.97 16.97 16.97
## 21 17.18 17.18 17.18 34.53 34.53 34.53 16.97 16.97 16.97
## 39 22.99 22.99 22.99 16.97 16.97 16.97 28.72 28.72 28.72
## 38 22.99 22.99 22.99 16.97 16.97 16.97 28.72 28.72 28.72
## 37 22.99 22.99 22.99 16.97 16.97 16.97 28.72 28.72 28.72

## m3
m3_sld <- lmer(score ~ 1 + Machine + (1 + c1 + c2 || Worker), data=Machines, REML=FALSE)
V3_sld <- getME_V(m3_sld)
(V3_sld_W1 <- round(V3_sld[1:9, 1:9],2))

## 9 x 9 sparse Matrix of class "dgCMatrix"
##        1     2     3    19    20    21    39    38    37
## 1  38.63 38.63 38.63 19.83 19.83 19.83 10.23 10.23 10.23
## 2  38.63 38.63 38.63 19.83 19.83 19.83 10.23 10.23 10.23
## 3  38.63 38.63 38.63 19.83 19.83 19.83 10.23 10.23 10.23
## 19 19.83 19.83 19.83 29.23 29.23 29.23 19.63 19.63 19.63
## 20 19.83 19.83 19.83 29.23 29.23 29.23 19.63 19.63 19.63
## 21 19.83 19.83 19.83 29.23 29.23 29.23 19.63 19.63 19.63
## 39 10.23 10.23 10.23 19.63 19.63 19.63 38.83 38.83 38.83
## 38 10.23 10.23 10.23 19.63 19.63 19.63 38.83 38.83 38.83
## 37 10.23 10.23 10.23 19.63 19.63 19.63 38.83 38.83 38.83

## m4
m4_sld <- lmer(score ~ 1 + Machine + (1 | Worker) + (0 + c2 || Worker), data=Machines, REML=FALSE)
V4_sld <- getME_V(m4_sld)
(V4_sld_W1 <- round(V4_sld[1:9, 1:9],2))

## 9 x 9 sparse Matrix of class "dgCMatrix"
##        1     2     3    19    20    21    39    38    37
## 1  23.33 23.33 23.33 23.33 23.33 23.33 19.98 19.98 19.98
## 2  23.33 23.33 23.33 23.33 23.33 23.33 19.98 19.98 19.98
## 3  23.33 23.33 23.33 23.33 23.33 23.33 19.98 19.98 19.98
## 19 23.33 23.33 23.33 23.33 23.33 23.33 19.98 19.98 19.98
## 20 23.33 23.33 23.33 23.33 23.33 23.33 19.98 19.98 19.98
## 21 23.33 23.33 23.33 23.33 23.33 23.33 19.98 19.98 19.98
## 39 19.98 19.98 19.98 19.98 19.98 19.98 26.67 26.67 26.67
## 38 19.98 19.98 19.98 19.98 19.98 19.98 26.67 26.67 26.67
## 37 19.98 19.98 19.98 19.98 19.98 19.98 26.67 26.67 26.67

Vs <- list(V2_sld, V3_sld, V4_sld)
for(i in 1:3) print(image(Vs[[i]], asp=1))

Model complexity decreases from left to right

5 Models for Machine x Worker interaction

## m6 / fm5
m6 <- lmer(score ~ 1 + Machine + (1 | Worker) + (0 + A + B + C || Worker), data=Machines, REML=FALSE)
V6 <- getME_V(m6)
(V6_W1 <- round(V6[1:9, 1:9],2))

## 9 x 9 sparse Matrix of class "dgCMatrix"
##        1     2     3    19    20    21    39    38    37
## 1  15.03 15.03 15.03 11.94 11.94 11.94 11.94 11.94 11.94
## 2  15.03 15.03 15.03 11.94 11.94 11.94 11.94 11.94 11.94
## 3  15.03 15.03 15.03 11.94 11.94 11.94 11.94 11.94 11.94
## 19 11.94 11.94 11.94 40.65 40.65 40.65 11.94 11.94 11.94
## 20 11.94 11.94 11.94 40.65 40.65 40.65 11.94 11.94 11.94
## 21 11.94 11.94 11.94 40.65 40.65 40.65 11.94 11.94 11.94
## 39 11.94 11.94 11.94 11.94 11.94 11.94 18.64 18.64 18.64
## 38 11.94 11.94 11.94 11.94 11.94 11.94 18.64 18.64 18.64
## 37 11.94 11.94 11.94 11.94 11.94 11.94 18.64 18.64 18.64

## m7 / fm3
m7 <- lmer(score ~ 1 + Machine + (1 | Worker) + (1 | Machine:Worker), data=Machines, REML=FALSE)
V7 <- getME_V(m7)
(V7_W1 <- round(V7[1:9, 1:9],2))

## 9 x 9 sparse Matrix of class "dgCMatrix"
##        1     2     3    19    20    21    39    38    37
## 1  30.59 30.59 30.59 19.05 19.05 19.05 19.05 19.05 19.05
## 2  30.59 30.59 30.59 19.05 19.05 19.05 19.05 19.05 19.05
## 3  30.59 30.59 30.59 19.05 19.05 19.05 19.05 19.05 19.05
## 19 19.05 19.05 19.05 30.59 30.59 30.59 19.05 19.05 19.05
## 20 19.05 19.05 19.05 30.59 30.59 30.59 19.05 19.05 19.05
## 21 19.05 19.05 19.05 30.59 30.59 30.59 19.05 19.05 19.05
## 39 19.05 19.05 19.05 19.05 19.05 19.05 30.59 30.59 30.59
## 38 19.05 19.05 19.05 19.05 19.05 19.05 30.59 30.59 30.59
## 37 19.05 19.05 19.05 19.05 19.05 19.05 30.59 30.59 30.59

## m8 / fm2
m8 <- lmer(score ~ 1 + Machine + (1 | Machine:Worker), data=Machines, REML=FALSE)
V8 <- getME_V(m8)
(V8_W1 <- round(V8[1:9, 1:9],2))

## 9 x 9 sparse Matrix of class "dgCMatrix"
##        1     2     3    19    20    21    39    38    37
## 1  30.59 30.59 30.59  .     .     .     .     .     .   
## 2  30.59 30.59 30.59  .     .     .     .     .     .   
## 3  30.59 30.59 30.59  .     .     .     .     .     .   
## 19  .     .     .    30.59 30.59 30.59  .     .     .   
## 20  .     .     .    30.59 30.59 30.59  .     .     .   
## 21  .     .     .    30.59 30.59 30.59  .     .     .   
## 39  .     .     .     .     .     .    30.59 30.59 30.59
## 38  .     .     .     .     .     .    30.59 30.59 30.59
## 37  .     .     .     .     .     .    30.59 30.59 30.59

Vs <- list(V6, V7, V8)
for(i in 1:3) print(image(Vs[[i]], asp=1))

Model complexity decreases from left to right

6 Images

6.1 Across models for factor levels (when estimated)

Row 1: m1a (identical with m1), m2a m5

Row 2: m3a, m4a, and m8

6.2 Across models and contrasts

Rows: models m1, m2 , m3, m4, m5

Columns: contr.treatment(), contr.sum, MASS::contr.sdif()

Model complexity decreases across rows within columns

7 Summary

There is always independence of observations betwen workers.
For models with RE0 specification, the structure of Vis (often) simpler, but, in a way, that is not surprising because the machines are all treated equally.
For models with RE1 specification, there is never independence of observations between machines from the same worker. This, again, is not that surprsing given that effects estimate differences between machines, that is, there are always two of the machines involved. At least for the treatment contrast it looks like the hiearchy is reflected in the patterning of covariance values. For sum and sliding contrasts this was not visible in the patterns.
The interaction LMM m7 also “spreads covariance” within subjects across machines.

In summary, the absence of a correlation parameter in the RE structure of an LMM does not imply that there are no correlations between random effects predicted by model parameters (and observstions). Indeed, whenever the LMM specification estimates parameters for differences between levels rather than levels, there will be correlations between observations across the levels of the fixed factor within the levels of the random factor.

8 Acknowledgment

This script was motivated and inspired by exchanges with Rune H. B. Christensen. I gratefully acknowledge his insights and comments.

9 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] lme4_1.1-18.9000 Matrix_1.2-14   
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.18      lattice_0.20-35   digest_0.6.17    
##  [4] MASS_7.3-50       grid_3.5.1        nlme_3.1-137     
##  [7] magrittr_1.5      evaluate_0.11     stringi_1.2.4    
## [10] minqa_1.2.4       nloptr_1.0.4      rmarkdown_1.10.11
## [13] splines_3.5.1     htmldeps_0.1.1    tools_3.5.1      
## [16] stringr_1.3.1     yaml_2.2.0        compiler_3.5.1   
## [19] htmltools_0.3.6   knitr_1.20

Images of Contrast-Dependent Random-Effect Structures in Linear Mixed Models

Reinhold Kliegl

2018-06-02 (revised: 2018-10-02 23:15:18)

1 Overview

2 Function: `getME_V()`

3 Level-based random-effect structures

3.1 Maximal model `m1a`

3.2 Models `m2a`, `m3a`, and `m4a`

4 Contrast-based random-effect structures

4.1 Treatment contrast

4.2 Sum contrast

4.3 Sliding-difference contrast

5 Models for Machine x Worker interaction

6 Images

6.1 Across models for factor levels (when estimated)

6.2 Across models and contrasts

7 Summary

8 Acknowledgment

9 Appendix

Images of Contrast-Dependent Random-Effect Structures in Linear Mixed Models

Reinhold Kliegl

2018-06-02 (revised: 2018-10-02 23:15:18)

1 Overview

2 Function: getME_V()

3 Level-based random-effect structures

3.1 Maximal model m1a

3.2 Models m2a, m3a, and m4a

4 Contrast-based random-effect structures

4.1 Treatment contrast

4.2 Sum contrast

4.3 Sliding-difference contrast

5 Models for Machine x Worker interaction

6 Images

6.1 Across models for factor levels (when estimated)

6.2 Across models and contrasts

7 Summary

8 Acknowledgment

9 Appendix

2 Function: `getME_V()`

3.1 Maximal model `m1a`

3.2 Models `m2a`, `m3a`, and `m4a`