Kesetaraan spesifikasi efek acak (0 + faktor | kelompok) dan (1 | kelompok) + (1 | kelompok: faktor) jika simetri senyawa

13

Douglas Bates menyatakan bahwa model-model berikut ini setara "jika matriks varians-kovarians untuk efek acak bernilai vektor memiliki bentuk khusus, yang disebut simetri gabungan" ( slide 91 dalam presentasi ini ):

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

Khususnya Bates menggunakan contoh ini:

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

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

dengan output yang sesuai:

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

Adakah yang bisa menjelaskan perbedaan antara model dan bagaimana m1mengurangi ke m2(diberikan simetri gabungan) dengan cara yang intuitif?

statmerkur
sumber
6
+1 dan, ya, ini benar-benar topik. Pilih untuk membuka kembali.
Amuba kata Reinstate Monica
2
@Peter Flom mengapa Anda menganggap pertanyaan ini sebagai di luar topik?
statmerkur
3
Mungkin tidak jelas bahwa Anda bertanya tentang model daripada tentang lme4sintaks. Akan sangat membantu - & memperlebar kumpulan penjawab potensial - jika Anda menjelaskannya untuk orang yang belum terbiasa lme4.
Scortchi
Sepertinya ini tentang coding.
Peter Flom - Reinstate Monica
1
Jika ini berguna, berikut adalah dua posting bagus tentang apa yang dilakukan sintaks lme4, dan apa simetri majemuk dalam konteks model campuran (lihat jawaban yang diterima pada kedua pertanyaan). stats.stackexchange.com/questions/13166/rs-lmer-cheat-sheet dan stats.stackexchange.com/questions/15102/…
Jacob Socolar

Jawaban:

11

Inn1nnyμγγy

yN(μ+γ,σy2I54)

σy2 adalah varian "residual".

γ

yN(μ,σy2I54+Σ)

ΣγΣ

m1

γ=Zθ

Z=I1813θT=[θ1,A,θ1,B,θ1,Cθ6,A,θ6,B,θ6,C]

θN(0,I6Λ)

ΛΛσθτ

Λ=[σθ2+τ2τ2τ2τ2σθ2+τ2τ2τ2τ2σθ2+τ2]

Λ memiliki semua elemen pada set offdiagonal dengan nilai yang sama.)

Σ=Z(I6Λ)ZTσθ2+τ2+σy2i,ju,v

cov(yi,u,yj,v)={0if ijτ2if i=j,uvσθ2+τ2if i=j,u=v

For your m2, the random effects decompose into:

γ=Zω+Xη

Where Z is as before, X=I619 is a design matrix that maps the random intercepts per worker onto observations, ωT=[ω1,A,ω1,B,ω1,C,,ω6,A,ω6,B,ω6,C] is the 18-dimensional vector of random intercepts for every combination of machine and worker; and ηT=[η1,,η6] is the 6-dimensional vector of random intercepts for worker. These are distributed as,

ηN(0,ση2I6)
ωN(0,σω2I18)
Where ση2,σω2 are the variances of these random intercepts.

The marginal covariance structure of m2 is Σ=σω2ZZT+ση2XXT, so that the variance of a given observation is σω2+ση2+σy2, and the covariance between two observations from workers i,j and machines u,v is:

cov(yi,u,yj,v)={0if ijση2if i=j,uvσω2+ση2if i=j,u=v

So ... σθ2σω2 and τ2ση2. 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
Nate Pope
sumber
1
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.
statmerkur
@statmerkur Thanks, I've tried to clarify this line. Let me know if you have another suggestion.
Nate Pope
1
Should X be defined as X=I619?
S. Catterall Reinstate Monica
1
@S.Catterall Yup, that's a typo -- thanks for catching it! I've fixed in my answer.
Nate Pope
2
@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.
Nate Pope