Tindakan berulang anova: lm vs lmer

10

Saya mencoba mereproduksi beberapa uji interaksi antara keduanya lmdan lmerpada tindakan berulang (2x2x2). Alasan saya ingin membandingkan kedua metode ini adalah karena GLM SPSS untuk tindakan berulang menghasilkan hasil yang sama persis dengan lmpendekatan yang disajikan di sini, jadi pada akhirnya saya ingin membandingkan SPSS vs R-lmer. Sejauh ini, saya hanya berhasil mereproduksi (dekat) beberapa interaksi ini.

Anda akan menemukan skrip di bawah ini untuk menggambarkan maksud saya dengan lebih baik:

library(data.table)
library(tidyr)
library(lmerTest)
library(MASS)

set.seed(1)

N     <- 100 # number of subjects
sigma <- 1   # popuplation sd
rho   <- .6  # correlation between variables

# X1:   a  a  a  a  b  b  b  b
# X2:   a  a  b  b  a  a  b  b
# X3:   a  b  a  b  a  b  a  b
mu <- c(5, 3, 3, 5, 3, 5, 5, 3) # means

# Simulate the data
sigma.mat <- rep(sigma, length(mu))
S <- matrix(sigma.mat, ncol = length(sigma.mat), nrow = length(sigma.mat))
Sigma <- t(S) * S * rho  
diag(Sigma) <- sigma**2
X <- data.table( mvrnorm(N, mu, Sigma) )
setnames(X, names(X), c("aaa", "aab", "aba", "abb", "baa", "bab", "bba", "bbb"))
X[, id := 1:.N]

# Long format
XL <- data.table( gather(X, key, Y, aaa:bbb) )
XL[, X1 := substr(key, 1, 1)]
XL[, X2 := substr(key, 2, 2)]
XL[, X3 := substr(key, 3, 3)]

# Recode long format (a = +1; b = -1)
XL[, X1c := ifelse(X1 == "a", +1, -1)]
XL[, X2c := ifelse(X2 == "a", +1, -1)]
XL[, X3c := ifelse(X3 == "a", +1, -1)]


### Composite scores to be used with lm
# X2:X3 2-way interaction (for half the data; i.e. when X1 == "a")
X[, X1a_X2.X3 := (aaa - aab) - (aba - abb)]

# X2:X3 2-way interaction (for all the data)
X[, aa := (aaa + baa) / 2]
X[, ab := (aab + bab) / 2]
X[, ba := (aba + bba) / 2]
X[, bb := (abb + bbb) / 2]
X[, X2.X3 := (aa - ab) - (ba - bb)]

# X1:X2:X3 3-way interaction (for all the data)
X[, X1.X2.X3 := ( (aaa - aab) - (aba - abb) ) - ( (baa - bab) - (bba - bbb) )]


### Fit models
# X2:X3 2-way interaction (for half the data; i.e. when X1 == "a")
summary( lm(X1a_X2.X3 ~ 1, X) ) # t = 34.13303
summary( lmer(Y ~ X2c*X3c + (X2c+X3c|id), XL[X1 == "a"]) ) # t = 34.132846  close match
summary( lmer(Y ~ X2c*X3c + (X2c+X3c||id), XL[X1 == "a"]) ) # t = 34.134624  close match

# X2:X3 2-way interaction (for all the data) 
summary( lm(X2.X3 ~ 1, X) ) # t = 0.3075025
summary( lmer(Y ~ X2c*X3c + (X2c+X3c|id), XL) ) # t = 0.1641932
summary( lmer(Y ~ X2c*X3c + (X2c+X3c||id), XL) ) # t = 0.1640710
summary( lmer(Y ~ X2c*X3c + (X2c*X3c|id), XL) ) # t = 0.1641765
anova(   lmer(Y ~ X2c*X3c + (X2c*X3c|id), XL), ddf = "Kenward-Roger" ) # t = 0.1643168
summary( lmer(Y ~ X2c*X3c + (X2c*X3c|id), XL, REML = FALSE) ) # t = 0.1645303
summary( lmer(Y ~ X2c*X3c + (X2c*X3c||id), XL) ) # t = 0.1640704

# X1:X2:X3 3-way interaction (for all the data)
summary( lm(X1.X2.X3 ~ 1, X) ) # t = 46.50177
summary( lmer(Y ~ X1c*X2c*X3c + (X1c*X2c*X3c - X1c:X2c:X3c|id), XL) ) # t = 49.0317599
anova(   lmer(Y ~ X1c*X2c*X3c + (X1c*X2c*X3c - X1c:X2c:X3c|id), XL), ddf = "Kenward-Roger" ) # t = 49.03176
summary( lmer(Y ~ X1c*X2c*X3c + (X1c*X2c*X3c - X1c:X2c:X3c|id), XL, REML = FALSE) ) # t = 49.2677606
summary( lmer(Y ~ X1c*X2c*X3c + (X1c*X2c*X3c - X1c:X2c:X3c||id), XL) ) # t = 46.5193774 close match
summary( lmer(Y ~ X1c*X2c*X3c + (X1c*X2c*X3c|id), XL) ) # unidentifiable
summary( lmer(Y ~ X1c*X2c*X3c + (X1c*X2c*X3c|id), XL,
              control = lmerControl(check.nobs.vs.nRE="ignore")) ) # t = 46.5148684 close match

Seperti yang dapat Anda lihat dari atas, tidak ada lmestimasi yang sama persis dengan yang lmerdiestimasi. Meskipun beberapa hasil sangat mirip dan mungkin berbeda hanya karena alasan numerik / komputasi. Kesenjangan antara kedua metode estimasi khusus besar untuk X2:X3 2-way interaction (for all the data).

Pertanyaan saya adalah apakah ada cara untuk mendapatkan hasil yang sama persis dengan kedua metode, dan jika ada cara yang benar untuk melakukan analisis dengan lmer(meskipun mungkin tidak cocok dengan lmhasilnya).


Bonus:

Saya perhatikan bahwa yang t valueterkait dengan interaksi 3-arah dipengaruhi oleh faktor-faktor yang dikodekan, yang tampaknya sangat aneh bagi saya:

summary( lmer(Y ~ X1*X2*X3 + (X1*X2*X3 - X1:X2:X3||id), XL) ) # t = 48.36
summary( lmer(Y ~ X1c*X2c*X3c + (X1c*X2c*X3c - X1c:X2c:X3c||id), XL) ) # t = 56.52
tikar
sumber
1
+1 karena terlihat menarik tetapi saya tidak tahu apa yang Anda lakukan di sini :) Dapatkah Anda menjelaskan dalam kata-kata atau matematika mengapa panggilan lm dan lmer ini menghasilkan koefisien yang sama? Dan apa logika di balik seluruh latihan ini?
amoeba
@amoeba Saya memperbarui posting saya untuk memperjelas tujuan posting ini. Pada dasarnya, saya ingin mereproduksi hasil dari SPSS (yang dapat diterjemahkan ke dalam lmmodel) dengan lmer, dan juga tahu apa analisis yang benar lmer untuk data jenis ini.
mat
Alasan perbedaan besar dalam hal interaksi dua arah untuk data lengkap adalah bahwa Anda memiliki 2 titik data per kombinasi parameter. Intuisi adalah bahwa ukuran sampel efektif untuk model campuran 2x lebih kecil dari untuk lm; Saya menduga itu sebabnya statistik-t kira-kira dua kali lebih kecil di lmer. Anda mungkin akan dapat mengamati fenomena yang sama menggunakan desain 2x2 yang lebih sederhana dan melihat efek utama, tanpa repot dengan 2x2x2 dan interaksi yang rumit.
amoeba

Jawaban:

3

Aneh, ketika saya menggunakan model terakhir Anda, saya menemukan pasangan yang cocok, bukan pasangan yang cocok:

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)  3.91221    0.07242 99.00001  54.025   <2e-16 ***
X1c          0.03277    0.05006 99.00000   0.655    0.514    
X2c         -0.04836    0.04644 99.00000  -1.042    0.300    
X3c          0.04248    0.05009 99.00001   0.848    0.398    
X1c:X2c      0.08370    0.08747 98.99998   0.957    0.341    
X1c:X3c     -0.07025    0.08895 98.99994  -0.790    0.432    
X2c:X3c     -0.02957    0.09616 99.00000  -0.308    0.759    
X1c:X2c:X3c -8.14099    0.17507 99.00003 -46.502   <2e-16 ***
pengguna244839
sumber
1
Untuk lebih jelasnya, model mana yang Anda maksud?
mat
ringkasan (lmer (Y ~ X1c X2c X3c + (X1c X2c X3c | id), XL, control = lmerControl (check.nobs.vs.nRE = "abaikan")))
user244839
Ini memang sangat aneh! summary( lmer(Y ~ X1c*X2c*X3c + (X1c*X2c*X3c|id), XL, control=lmerControl(check.nobs.vs.nRE="ignore")) )$coefficientskembali t = 46.5148684untukku. Mungkinkah ada masalah versi? Saya menggunakan R version 3.5.3 (2019-03-11)dan lmerTest 3.1-0.
mat
Saya memiliki versi R & lmerTest yang sama dengan @mat dan mendapatkan hasil yang sama seperti mereka (walaupun dengan banyak peringatan - kegagalan untuk bertemu, dll).
mkt - Pasang kembali Monica
1
@ ah Mungkin saya tidak jelas - saya mendapatkan hasil yang sama seperti Anda! Saya pikir Anda mungkin benar bahwa user244839 menggunakan versi yang berbeda dari kami.
mkt - Pasang kembali Monica