1
$\begingroup$

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?

$\endgroup$
0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Thanks for this clear depiction! I was thinking of custom contrasts and made the corresponding edit (should have made this clear earlier as it seems also relevant to our email exchange; probably using "standard contrasts" wasn't the best choice). Indeed I'm thinking of the contrasts as fixed slopes and their corresponding random slopes. I got this idea and the corresponding mathematical notation from Matuschek et al. (including Reinhold Kliegl) (2017) and Barr et al. (2013). $\endgroup$ – statmerkur Jun 1 '18 at 9:19
  • $\begingroup$ I don't think which contrasts are used really affects my answer. I have added an illustration of custom contrasts but any other contrasts should be easy to work out. $\endgroup$ – Rune H Christensen Jun 1 '18 at 10:19
  • $\begingroup$ I see how V for m2a (with contr.treatment) corresponds to print(VarCorr(m2a), comp="Variance"). But what is the relation between V for m2aX and print(VarCorr(m2aX), comp="Variance")? $\endgroup$ – statmerkur Jun 2 '18 at 14:33
  • $\begingroup$ If I'm interested in specific contrasts and their associated (co)variances I think m1a and m2a are the models that would tell me something about that and I do not see how models that use (0+A+B+C) in the lmer random parts can achieve the same. Am I missing something here? $\endgroup$ – statmerkur Jun 2 '18 at 14:41
  • $\begingroup$ Besides that, a fundamental Q on your "remark on mathematical notation": When categorical covariates are transformed into numerical covariates why should these numerical covariates and their corresponding random coefficients/slopes be expressed in a different way than continuous variable are? $\endgroup$ – statmerkur Jun 2 '18 at 14:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.