Multilevel regression form for lmm with different lmer factor representations in fixed and random part

Question

Consider a linear mixed model with a categorical predictor (factor) with levels A, B, and C, the corresponding custom contrasts c1 and c2, and the grouping factor group. Is there a way to express models like m3a/m3b and m4a/m4b (which are discussed here), where the representation of the factor in the fixed and random part "don't match", in hierarchical/multilevel regression form?

Update: In what follows c1 refers to the difference A - (B + C)/2 and c2 to B - C. The distinction between version a and b is made for illustration purposes, these models are always equivalent.

contrasts(data$factor) <- MASS::ginv(rbind(c(1, -0.5, -0.5),  # A - (B + C)/2
                                           c(0, 1, -1)))      # B - C

mm1 <- model.matrix(~ 1 + factor, data)  # depends on the contrasts
c1 <- mm1[, 2]
c2 <- mm1[, 3] 

mm0 <- model.matrix(~ 0 + factor, data)  # does NOT depend on the contrasts
A <- mm0[, 1]
B <- mm0[, 2]
C <- mm0[, 3] 

m1a <- lmer(y ~ 1 + factor + (1 + c1 + c2 | group), data)
m1b <- lmer(y ~ 1 + c1 + c2 + (1 + c1 + c2 | group), data) 
m2a <- lmer(y ~ 1 + factor + (1 + c1 + c2 || group), data) 
m2b <- lmer(y ~ 1 + c1 + c2 + (1 + c1 + c2 || group), data)
m3a <- lmer(y ~ 1 + factor + (0 + A + B + C | group), data)
m3b <- lmer(y ~ 1 + c1 + c2 + (0 + A + B + C | group), data) 
m4a <- lmer(y ~ 1 + factor + (1 + A + B + C || group), data)
m4b <- lmer(y ~ 1 + c1 + c2  + (1 + A + B + C || group), data)

m1a/m1bcorrespond to:

\begin{aligned} &Y_{c,g,i}\ =\ \beta_0 + G_{0,g} + \left(\beta_1 + G_{1,g}\right)X_{1,c} +\left(\beta_2 + G_{2,g}\right)X_{2,c} + \epsilon_{c,g,i},\\ &\begin{pmatrix} G_{0,g}\\ G_{1,g}\\ G_{2,g} \end{pmatrix} \sim N \begin{bmatrix} \begin{pmatrix} 0\\ 0\\ 0 \end{pmatrix}, \begin{pmatrix} \tau_{00}^2 & \rho_{g_{01}}\tau_{00}\tau_{11} & \rho_{g_{02}}\tau_{00}\tau_{22}\\ \rho_{g_{01}}\tau_{00}\tau_{11} & \tau_{11}^2 & \rho_{g_{12}}\tau_{11}\tau_{22}\\ \rho_{g_{02}}\tau_{00}\tau_{22} & \rho_{g_{12}}\tau_{11}\tau_{22} & \tau_{22}^2 \end{pmatrix} \end{bmatrix},\\ &\epsilon_{c,g,i}\ \sim\ N \left(0,\ \sigma^2\right). \end{aligned}

$Y_{c,g,i}$ refers to the dependent variable, the subscripts stand for condition (c) and group (g).The by-group intercepts and slopes are $G_{0,g}$, $G_{1,g}$ and $G_{2,g}$, respectively. The fixed effects are $\beta_0$, $\beta_1$ and $\beta_2$ where $\beta_1$ and $\beta_2$ represents the coefficients/slopes estimated for c1 and c2, respectively. The residual error is $\epsilon_{c,g,i}$.

m2a/m2b correspond to:

\begin{aligned} &Y_{c,g,i}\ =\ \beta_0 + G_{0,g} + \left(\beta_1 + G_{1,g}\right)X_{1,c} +\left(\beta_2 + G_{2,g}\right)X_{2,c} + \epsilon_{c,g,i},\\ &\begin{pmatrix} G_{0,g}\\ G_{1,g}\\ G_{2,g} \end{pmatrix} \sim N \begin{bmatrix} \begin{pmatrix} 0\\ 0\\ 0 \end{pmatrix}, \begin{pmatrix} \tau_{00}^2 & 0 & 0\\ 0 & \tau_{11}^2 & 0\\ 0 & 0 & \tau_{22}^2 \end{pmatrix} \end{bmatrix},\\ &\epsilon_{c,g,i}\ \sim\ N \left(0,\ \sigma^2\right). \end{aligned}

But what about m3a/m3b and m4a/m4b?

Answer 1

The statistical model for all your models (m1a $\ldots$ m4b) can be written on the following "standard" form: $$y_{ijk} = \mu + \alpha_j + d_{ij} + e_{ijk}, \quad d_{i} \sim N(0, \Sigma), \quad e_{ijk} \sim N(0, \sigma^2)$$ where $\alpha_j$ is the fixed-effect of the $j$'th treatment (factor), and $d_i = (d_{i1}, \ldots, d_{iJ})^\top$ is the vector-random-effect for the $i$'th subject (group) possibly depending on the $j$'th treatment.

The only difference between your models is how $\Sigma$ is parameterised.

Let's first consider the fixed-effects structure:

When fitting the model it is brought on the general form: $$y = X\beta + Z b + e$$ where the fixed-effects coefficients $\beta$ are functions of the parameters $\{\mu, \alpha_i\}$. The standard treatment contrast coding in R means that $X$ is set up such that $\beta$ represents the following function of the parameters: $\beta = (\mu + \alpha_1, \alpha_2 - \alpha_1, \alpha_3 - \alpha_1)^\top$ assuming $J = 3$ treatments but there are many other contrast coding options. From the coefficient vector, $\beta$ you can extract all estimable contrasts (linear functions) of the parameters that you would like or you can, as you suggest in the code, parameterise $X$ differently such that $\beta$ represent different estimable functions of the parameters. To uniquely estimate a vector $\beta$ it is required that $X$ has full rank - this is why it only has 3 columns even though there are 4 parameters.

Now let's turn to the random-effect structure:

A key point here is that $Z$ is just an indicator matrix - no "coding" (treatment contrasts or otherwise) is required and it is not required to be of full rank: $Z$ has a column for each $d_{ij}$. If there are 3 treatments and 5 subjects, $Z$ will have 15 columns. Thus, using lmer syntax such as (1 + c1 + c2 || group) may affect the way that lmer sets up its computations, but in terms of the model expression above $Z$ remains the same indicator matrix, instead it is the structure of $\Sigma$ that changes.

A couple of examples are

1. (factor | group), (1 + c1 + c2 | group) and (0 + A + B + C | group) all parameterise a symmetric but otherwise unconstrained $\Sigma=\begin{bmatrix}\sigma^2_A & \sigma_{AB} & \sigma_{AC} \\ \sigma_{AB} & \sigma^2_{B} & \sigma^2_{BC} \\\sigma_{AC} & \sigma_{BC} & \sigma^2_{C} \end{bmatrix}$ using 6 parameters.

2. (0 + A + B + C || group) constrains $\Sigma$ to be diagonal and use 3 parameters for that; $\Sigma=\begin{bmatrix}\sigma^2_A & 0 & 0 \\ 0 & \sigma^2_{B} & 0 \\ 0 & 0 & \sigma^2_{C} \end{bmatrix}$.

3. (1 + A + B + C || group) which is the same as (1 | group) + (0 + A + B + C || group) constrains the off-diagonal entries in $\Sigma$ to be the same but with different variances in the diagonal. 4 parameters are used for that, $\Sigma=\begin{bmatrix}\sigma^2_A + \sigma^2_X & \sigma^2_X & \sigma^2_X \\ \sigma^2_X & \sigma^2_{B} + \sigma^2_X & \sigma^2_X \\ \sigma^2_X & \sigma^2_X & \sigma^2_{C} + \sigma^2_X \end{bmatrix}$.

Working out the structure and parameterisation of $\Sigma$ for, say, (1 + c1 + c2 || group) is more difficult but with a treatment contrast coding of c1 and c2 it can be worked out that 2nd and 3rd diagonal entries are free parameters, but that all off-diagonal entries are restricted to equal the 1st diagonal entry; $\Sigma=\begin{bmatrix}\sigma^2_X & \sigma^2_X & \sigma^2_X \\ \sigma^2_X & \sigma^2_{B} + \sigma^2_X & \sigma^2_X \\ \sigma^2_X & \sigma^2_X & \sigma^2_{C} + \sigma^2_X \end{bmatrix}$. I cannot think of a situation where this is a natural model but it would be interesting if you know of one.

We can work that last structure out analytically, but we can also simply demonstrate it using lmer. Here I use the Machines dataset as an example:

library("lme4")
data("Machines", package = "MEMSS")
Machines <- Machines[do.call(order, Machines), ] # reorder by Worker, then Machine
mm <- model.matrix(~ 1 + Machine, Machines)
c1 <- mm[, 2]
c2 <- mm[, 3]
m2a <- lmer(score ~ 1 + Machine + (1 + c1 + c2 || Worker), data=Machines, REML=FALSE)
Z <- getME(m2a, "Z")
Lambda <- getME(m2a, "Lambda")
V <- sigma(m2a)^2 * Z %*% Lambda %*% t(Lambda) %*% t(Z)
round(V[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

This matrix shows the random-effects contribution to the marginal variance-covariance matrix, i.e., the first term in $\mathsf{Cov}(y_i) = Z_i \Sigma Z_i^\top + \sigma^2 I$. This is the same matrix for all Workers (in this nicely balanced example) and the code just extracts it for the first Worker. Because there are 3 observations on each Machine-Worker combination each element of $\Sigma$ is scaled up to a ($3\times 3$) sub-matrix.

Using alternative contrasts

Using alternative contrasts changes the model, but only the variance-covariance structure of the random-effects (as long as $X$ is set up with 3 linearly independent columns). Using your custom contrasts we get

Machines$Machine2 <- Machines$Machine
contrasts(Machines$Machine2) <- MASS::ginv(rbind(c(1, -0.5, -0.5),  # A - (B + C)/2
                                                c(0, 1, -1))) 
mm <- model.matrix(~ 1 + Machine2, Machines)
c1 <- mm[, 2]
c2 <- mm[, 3]
m2aX <- lmer(score ~ 1 + Machine + (1 + c1 + c2 || Worker), data=Machines, REML=FALSE)
m2aZ <- lmer(score ~ 1 + c1 + c2 + (1 + c1 + c2 || Worker), data=Machines, REML=FALSE)
anova(m2a, m2aX, m2aZ)

Data: Machines
Models:
m2a: score ~ 1 + Machine + ((1 | Worker) + (0 + c1 | Worker) + (0 + 
m2a:     c2 | Worker))
m2aX: score ~ 1 + Machine + ((1 | Worker) + (0 + c1 | Worker) + (0 + 
m2aX:     c2 | Worker))
m2aZ: score ~ 1 + c1 + c2 + ((1 | Worker) + (0 + c1 | Worker) + (0 + 
m2aZ:     c2 | Worker))
     Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
m2a   7 235.11 249.04 -110.56   221.11                        
m2aX  7 238.84 252.76 -112.42   224.84     0      0          1
m2aZ  7 238.84 252.76 -112.42   224.84     0      0          1

Z <- getME(m2aX, "Z")
Lambda <- getME(m2aX, "Lambda")
V <- sigma(m2aX)^2 * Z %*% Lambda %*% t(Lambda) %*% t(Z)
round(V[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 20.09 20.09 20.09 20.09 20.09 20.09
2  28.51 28.51 28.51 20.09 20.09 20.09 20.09 20.09 20.09
3  28.51 28.51 28.51 20.09 20.09 20.09 20.09 20.09 20.09
19 20.09 20.09 20.09 31.63 31.63 31.63 16.97 16.97 16.97
20 20.09 20.09 20.09 31.63 31.63 31.63 16.97 16.97 16.97
21 20.09 20.09 20.09 31.63 31.63 31.63 16.97 16.97 16.97
39 20.09 20.09 20.09 16.97 16.97 16.97 31.63 31.63 31.63
38 20.09 20.09 20.09 16.97 16.97 16.97 31.63 31.63 31.63
37 20.09 20.09 20.09 16.97 16.97 16.97 31.63 31.63 31.63

1. the anova output illustrates that changing the contrasts changes the model because the random-effects are set up differently - not because the coding of the fixed-effects has changed.

2. The last bit of output shows that the variance-covariance matrix for the random-effects has the following parameterisation: $\Sigma=\begin{bmatrix}\sigma^2_1 & \sigma^2_2 & \sigma^2_2 \\ \sigma^2_2 & \sigma^2_{3} & \sigma^2_4 \\ \sigma^2_2 & \sigma^2_4 & \sigma^2_{3}\end{bmatrix}$ with some dependence among the parameters as there are only 3 independent parameters. Again a form that I find somewhat peculiar.

A remark on mathematical notation

The way you have written up your models mathematically doesn't make sense to me and so I'm not sure how to reconcile your notation with mine. I may just be displaying my ignorance here, but it seems to me that your notation implies that the random-effects design matrix should be the same as the fixed-effects design matrix? The notation has similarities with a notation that I would use if there were a continuous variable for which we could estimate random slopes. But there are no continuous variables in this example and so it does not make sense to me to talk about random slopes.