14
$\begingroup$

Douglas Bates states that the following models are equivalent "if the variance-covariance matrix for the vector-valued random effects has a special form, called compound symmetry" (slide 91 in this presentation):

m1 <- lmer(y ~ factor + (0 + factor|group), data)
m2 <- lmer(y ~ factor + (1|group) + (1|group:factor), data)

Specifically Bates uses this example:

library(lme4)
data("Machines", package = "MEMSS")

m1a <- lmer(score ~ Machine + (0 + Machine|Worker), Machines)
m2a <- lmer(score ~ Machine + (1|Worker) + (1|Worker:Machine), Machines)

with the corresponding outputs:

print(m1a, corr = FALSE)

Linear mixed model fit by REML ['lmerMod']
Formula: score ~ Machine + (0 + Machine | Worker)
   Data: Machines
REML criterion at convergence: 208.3112
Random effects:
 Groups   Name     Std.Dev. Corr     
 Worker   MachineA 4.0793            
          MachineB 8.6253   0.80     
          MachineC 4.3895   0.62 0.77
 Residual          0.9616            
Number of obs: 54, groups:  Worker, 6
Fixed Effects:
(Intercept)     MachineB     MachineC  
     52.356        7.967       13.917  

print(m2a, corr = FALSE)

Linear mixed model fit by REML ['lmerMod']
Formula: score ~ Machine + (1 | Worker) + (1 | Worker:Machine)
   Data: Machines
REML criterion at convergence: 215.6876
Random effects:
 Groups         Name        Std.Dev.
 Worker:Machine (Intercept) 3.7295  
 Worker         (Intercept) 4.7811  
 Residual                   0.9616  
Number of obs: 54, groups:  Worker:Machine, 18; Worker, 6
Fixed Effects:
(Intercept)     MachineB     MachineC  
     52.356        7.967       13.917

Can anyone explain the difference between the models and how m1 reduces to m2 (given compound symmetry) in an intuitive way?

$\endgroup$
7
  • 7
    $\begingroup$ +1 and, imho, this is absolutely on topic. Vote to reopen. $\endgroup$
    – amoeba
    Commented Sep 21, 2017 at 20:36
  • 2
    $\begingroup$ @Peter Flom why do you consider this question as off-topic? $\endgroup$
    – statmerkur
    Commented Sep 21, 2017 at 21:13
  • 4
    $\begingroup$ It perhaps wasn't clear that you were asking about the models rather than about the lme4 syntax. It'd be helpful - & widen the pool of potential answerers - if you explained them for people unfamiliar with lme4. $\endgroup$ Commented Sep 21, 2017 at 22:22
  • 1
    $\begingroup$ I think this is almost literally the definition of compound symmetry. Is there some specific definition of compound symmetry that an answer should use, and then prove that that definition implies the equivalence of the models m1 and m2? $\endgroup$ Commented Sep 28, 2017 at 15:56
  • 2
    $\begingroup$ If it's useful, here are two good posts about what the lme4 syntax is doing, and what compound symmetry is in the context of mixed models (see accepted answers on both questions). stats.stackexchange.com/questions/13166/rs-lmer-cheat-sheet and stats.stackexchange.com/questions/15102/… $\endgroup$ Commented Sep 28, 2017 at 15:58

1 Answer 1

11
+75
$\begingroup$

In this example, there are three observations for each combination of the three machines (A, B, C) and the six workers. I'll use $I_n$ to denote a $n$-dimensional identity matrix and $1_n$ to denote an $n$-dimensional vector of ones. Let's say $y$ is the vector of observations, that I will assume is ordered by worker then machine then replicate. Let $\mu$ be the corresponding expected values (e.g. the fixed effects), and let $\gamma$ be a vector of group-specific deviations from the expected values (e.g. the random effects). Conditional on $\gamma$, the model for $y$ can be written:

$$y \sim \mathcal{N}(\mu + \gamma, \sigma^2_y I_{54})$$

where $\sigma^2_y$ is the "residual" variance.

To understand how the covariance structure of the random effects induces a covariance structure among observations, it is more intuitive to work with the equivalent "marginal" representation, which integrates over the random effects $\gamma$. The marginal form of this model is,

$$y \sim \mathcal{N}(\mu, \sigma^2_y I_{54} + \Sigma)$$

Here, $\Sigma$ is a covariance matrix that depends on the structure of $\gamma$ (e.g. the "variance components" underlying the random effects). I'll refer to $\Sigma$ as the "marginal" covariance.

In your m1, the random effects decompose as: $$\gamma = Z \theta$$

Where $Z = I_{18} \otimes 1_3$ is a design matrix that maps the random coefficients onto observations, and $\theta^T = [\theta_{1,A}, \theta_{1,B}, \theta_{1,C} \dots \theta_{6,A}, \theta_{6,B}, \theta_{6,C}]$ is the 18-dimensional vector of random coefficients ordered by worker then machine, and is distributed as:

$$\theta \sim \mathcal{N}(0, I_6 \otimes \Lambda)$$

Here $\Lambda$ is the covariance of the random coefficients. The assumption of compound symmetry means that $\Lambda$ has two parameters, that I'll call $\sigma_\theta$ and $\tau$, and the structure:

$$\Lambda = \left[\begin{matrix} \sigma^2_\theta + \tau^2 & \tau^2 & \tau^2 \\ \tau^2 & \sigma^2_\theta + \tau^2 & \tau^2 \\ \tau^2 & \tau^2 & \sigma^2_\theta + \tau^2 \end{matrix}\right]$$

(In other words, the correlation matrix underlying $\Lambda$ has all the elements on the offdiagonal set to the same value.)

The marginal covariance structure induced by these random effects is $\Sigma = Z (I_6 \otimes \Lambda) Z^T$, so that the variance of a given observation is $\sigma^2_\theta + \tau^2 + \sigma^2_y$ and the covariance between two (separate) observations from workers $i, j$ and machines $u, v$ is: $$\mathrm{cov}(y_{i,u}, y_{j,v}) = \begin{cases} 0 & \text{if } i\neq j \\ \tau^2 & \text{if } i=j, u\neq v \\ \sigma^2_\theta + \tau^2 & \text{if } i=j, u=v \end{cases}$$

For your m2, the random effects decompose into:

$$\gamma = Z \omega + X \eta$$

Where Z is as before, $X = I_6 \otimes 1_9$ is a design matrix that maps the random intercepts per worker onto observations, $\omega^T = [\omega_{1,A}, \omega_{1,B}, \omega_{1,C}, \dots, \omega_{6,A}, \omega_{6,B}, \omega_{6,C}]$ is the 18-dimensional vector of random intercepts for every combination of machine and worker; and $\eta^T = [\eta_{1}, \dots, \eta_{6}]$ is the 6-dimensional vector of random intercepts for worker. These are distributed as, $$\eta \sim \mathcal{N}(0, \sigma^2_\eta I_6)$$ $$\omega \sim \mathcal{N}(0, \sigma^2_\omega I_{18})$$ Where $\sigma_\eta^2, \sigma_\omega^2$ are the variances of these random intercepts.

The marginal covariance structure of m2 is $\Sigma = \sigma^2_\omega Z Z^T + \sigma^2_\eta X X^T$, so that the variance of a given observation is $\sigma^2_\omega + \sigma^2_\eta + \sigma^2_y$, and the covariance between two observations from workers $i, j$ and machines $u, v$ is: $$\mathrm{cov}(y_{i,u}, y_{j,v}) = \begin{cases} 0 & \text{if } i\neq j \\ \sigma_\eta^2 & \text{if } i=j,u\neq v \\ \sigma^2_\omega + \sigma^2_\eta & \text{if } i=j,u=v \end{cases}$$

So ... $\sigma^2_\theta \equiv \sigma^2_\omega$ and $\tau^2 \equiv \sigma^2_\eta$. If m1 assumed compound symmetry (which it doesn't with your call to lmer, because the random effects covariance is unstructured).

Brevity is not my strong point: this is all just a long, convoluted way of saying that each model has two variance parameters for the random effects, and are just two different ways of writing of the same "marginal" model.

In code ...

sigma_theta <- 1.8
tau         <- 0.5
sigma_eta   <- tau
sigma_omega <- sigma_theta
Z <- kronecker(diag(18), rep(1,3))
rownames(Z) <- paste(paste0("worker", rep(1:6, each=9)), 
                     rep(paste0("machine", rep(1:3, each=3)),6))
X <- kronecker(diag(6), rep(1,9))
rownames(X) <- rownames(Z)
Lambda <- diag(3)*sigma_theta^2 + tau^2

# marginal covariance for m1:
Z%*%kronecker(diag(6), Lambda)%*%t(Z)
# for m2:
X%*%t(X)*sigma_eta^2 + Z%*%t(Z)*sigma_omega^2
$\endgroup$
8
  • 1
    $\begingroup$ Very nice answer! But I think the phrase "machine nested within worker" could be misleading as the same three machines appear in more than one (in fact every) level of worker. $\endgroup$
    – statmerkur
    Commented Oct 10, 2017 at 9:44
  • $\begingroup$ @statmerkur Thanks, I've tried to clarify this line. Let me know if you have another suggestion. $\endgroup$
    – Nate Pope
    Commented Oct 10, 2017 at 18:01
  • 1
    $\begingroup$ Should $X$ be defined as $X = I_6 \otimes 1_9$? $\endgroup$ Commented Oct 18, 2017 at 16:57
  • 1
    $\begingroup$ @S.Catterall Yup, that's a typo -- thanks for catching it! I've fixed in my answer. $\endgroup$
    – Nate Pope
    Commented Oct 18, 2017 at 18:05
  • 2
    $\begingroup$ @statmerkur can you clarify what you mean? There's no continuous covariates here, so not sure what you mean by "slope". The way I think of the model is that there are systematic differences in the mean of the response between machines (the fixed effects); then a random deviation for each worker (random intercepts/worker); then a random deviation for each machine-worker combination; and finally a random deviation per observation. The greater the variance of the random deviations per worker, the more correlated observations from a given worker would be, etc. $\endgroup$
    – Nate Pope
    Commented Oct 18, 2017 at 18:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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