Equivalencia de las especificaciones de efectos aleatorios (0 + factor|grupo) y (1|grupo) + (1|grupo:factor) en caso de simetría compuesta

Question

Douglas Bates afirma que los siguientes modelos son equivalentes "si la matriz de varianza-covarianza para los efectos aleatorios vectoriales tiene una forma especial, llamada simetría compuesta" ( diapositiva 91 de esta presentación ):

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

En concreto, Bates utiliza este ejemplo:

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

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

con las correspondientes salidas:

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

¿Puede alguien explicar la diferencia entre los modelos y cómo m1 se reduce a m2 (dada la simetría compuesta) de forma intuitiva?

Answer 1

En este ejemplo, hay tres observaciones para cada combinación de las tres máquinas (A, B, C) y los seis trabajadores. Utilizaré $I_n$ para denotar un $n$ -y la matriz de identidad de las dimensiones $1_n$ para denotar un $n$ -vector de unos. Digamos que $y$ es el vector de observaciones, que supondré que está ordenado por trabajador, luego máquina y luego réplica. Sea $\mu$ sean los valores esperados correspondientes (por ejemplo, los efectos fijos), y que $\gamma$ sea un vector de desviaciones específicas del grupo respecto a los valores esperados (por ejemplo, los efectos aleatorios). Condicionado a $\gamma$ el modelo de $y$ se puede escribir:

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

donde $\sigma^2_y$ es la varianza "residual".

Para entender cómo la estructura de covarianza de los efectos aleatorios induce una estructura de covarianza entre las observaciones, es más intuitivo trabajar con la representación "marginal" equivalente que se integra sobre los efectos aleatorios $\gamma$ . La forma marginal de este modelo es,

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

Aquí, $\Sigma$ es una matriz de covarianza que depende de la estructura de $\gamma$ (por ejemplo, los "componentes de la varianza" subyacentes a los efectos aleatorios). Me referiré a $\Sigma$ como la covarianza "marginal".

En su m1 los efectos aleatorios se descomponen como $$\gamma = Z \theta$$

Donde $Z = I_{18} \otimes 1_3$ es una matriz de diseño que asigna los coeficientes aleatorios a las observaciones, y $\theta^T = [\theta_{1,A}, \theta_{1,B}, \theta_{1,C} \dots \theta_{6,A}, \theta_{6,B}, \theta_{6,C}]$ es el vector de 18 dimensiones de coeficientes aleatorios ordenados por trabajador y luego por máquina, y se distribuye como

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

Aquí $\Lambda$ es la covarianza de los coeficientes aleatorios. La suposición de simetría compuesta significa que $\Lambda$ tiene dos parámetros, que llamaré $\sigma_\theta$ y $\tau$ y la estructura:

$$\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]$$

(En otras palabras, la matriz de correlación subyacente $\Lambda$ tiene todos los elementos de la offdiagonal con el mismo valor).

La estructura de covarianza marginal inducida por estos efectos aleatorios es $\Sigma = Z (I_6 \otimes \Lambda) Z^T$ para que la varianza de una observación dada sea $\sigma^2_\theta + \tau^2 + \sigma^2_y$ y la covarianza entre dos observaciones (separadas) de los trabajadores $i, j$ y máquinas $u, v$ es: $$\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}$$

Para su m2 los efectos aleatorios se descomponen en

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

Donde Z es como antes, $X = I_6 \otimes 1_9$ es una matriz de diseño que asigna los interceptos aleatorios por trabajador a las observaciones, $\omega^T = [\omega_{1,A}, \omega_{1,B}, \omega_{1,C}, \dots, \omega_{6,A}, \omega_{6,B}, \omega_{6,C}]$ es el vector de 18 dimensiones de intercepciones aleatorias para cada combinación de máquina y trabajador; y $\eta^T = [\eta_{1}, \dots, \eta_{6}]$ es el vector de 6 dimensiones de interceptos aleatorios para el trabajador. Estos se distribuyen como, $$\eta \sim \mathcal{N}(0, \sigma^2_\eta I_6)$$ $$\omega \sim \mathcal{N}(0, \sigma^2_\omega I_{18})$$ Donde $\sigma_\eta^2, \sigma_\omega^2$ son las varianzas de estos interceptos aleatorios.

La estructura de covarianza marginal de m2 es $\Sigma = \sigma^2_\omega Z Z^T + \sigma^2_\eta X X^T$ para que la varianza de una observación dada sea $\sigma^2_\omega + \sigma^2_\eta + \sigma^2_y$ y la covarianza entre dos observaciones de los trabajadores $i, j$ y máquinas $u, v$ es: $$\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}$$

Así que... $\sigma^2_\theta \equiv \sigma^2_\omega$ y $\tau^2 \equiv \sigma^2_\eta$ . Si m1 asumió una simetría compuesta (que no tiene con su llamada a lmer, porque la covarianza de los efectos aleatorios no está estructurada).

La brevedad no es mi punto fuerte: todo esto no es más que una forma larga y enrevesada de decir que cada modelo tiene dos parámetros de varianza para los efectos aleatorios, y no son más que dos formas diferentes de escribir el mismo modelo "marginal".

En código...

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