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

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

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

Específicamente, Bates usa 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 salidas correspondientes:

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

¿Alguien puede explicar la diferencia entre los modelos y cómo se m1reduce a m2(simetría compuesta dada) de una manera intuitiva?

En este ejemplo, hay tres observaciones para cada combinación de las tres máquinas (A, B, C) y los seis trabajadores. Usaré para denotar una matriz de identidad n- dimensional y 1 n para denotar un vector n- dimensional de unos. Digamos que y es el vector de observaciones, que supongo que está ordenado por trabajador, luego máquina y luego replicar. Supongamos que μ son los valores esperados correspondientes (por ejemplo, los efectos fijos) y que γ sea ​​un vector de desviaciones específicas del grupo de los valores esperados (por ejemplo, los efectos aleatorios). Condicional a γ , el modelo para y se puede escribir:$I_n$$n$$1_n$$n$$y$$\mu$$\gamma$$\gamma$$y$

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

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

Para comprender 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 . La forma marginal de este modelo es,$\gamma$

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

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

En tu m1, los efectos aleatorios se descomponen como:

$$\gamma = Z \theta$$

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

$$\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 tiene dos parámetros, que llamaré σ θ y τ , y la estructura:$\Lambda$$\sigma_\theta$$\tau$

$$\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 tiene todos los elementos en el conjunto fuera del diagnóstico con el mismo valor).$\Lambda$

La estructura de covarianza marginal inducida por estos efectos aleatorios es , de modo que la varianza de una observación dada es σ 2 θ + τ 2 + σ 2 y y la covarianza entre dos observaciones (separadas) de los trabajadores i , j y máquinas u , v es: c o v ( y i , u , y j , v ) $\Sigma = Z (I_6 \otimes \Lambda) Z^T$$\sigma^2_\theta + \tau^2 + \sigma^2_y$$i, j$$u, v$

$$\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, es una matriz de diseño que mapea las intersecciones aleatorias por trabajador en observaciones, ω T = [ ω 1 , A , ω 1 , B , ω 1 , C , ... , ω 6 , A , ω 6 T = [ η 1 , … , η 6$X = I_6 \otimes 1_9$$\omega^T = [\omega_{1,A}, \omega_{1,B}, \omega_{1,C}, \dots, \omega_{6,A}, \omega_{6,B}, \omega_{6,C}]$$\eta^T = [\eta_{1}, \dots, \eta_{6}]$

$$\eta \sim \mathcal{N}(0, \sigma^2_\eta I_6)$$ $$\omega \sim \mathcal{N}(0, \sigma^2_\omega I_{18})$$ $\sigma_\eta^2, \sigma_\omega^2$

m2$\Sigma = \sigma^2_\omega Z Z^T + \sigma^2_\eta X X^T$$\sigma^2_\omega + \sigma^2_\eta + \sigma^2_y$$i, j$$u, v$

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

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