Saya menganalisis hasil percobaan waktu reaksi dalam R.

Saya menjalankan pengukuran ANOVA berulang (1 faktor dalam subjek dengan 2 level dan 1 antara faktor subjek dengan 2 level). Saya menjalankan model campuran linear yang sama dan saya ingin meringkas hasil lmer dalam bentuk tabel ANOVA menggunakan lmerTest::anova.

Jangan salah paham: Saya tidak mengharapkan hasil yang sama, namun saya tidak yakin tentang tingkat kebebasan dalam lmerTest::anovahasil. Menurut saya ini agak mencerminkan ANOVA tanpa agregasi pada tingkat subjek.

Saya menyadari fakta bahwa menghitung derajat kebebasan dalam model efek campuran itu rumit, tetapi lmerTest::anovadisebutkan sebagai salah satu solusi yang mungkin dalam ?pvaluestopik yang diperbarui ( lme4paket).

Apakah perhitungan ini benar? Apakah hasil lmerTest::anovamencerminkan model yang ditentukan dengan benar?

Pembaruan: Saya membuat perbedaan individu lebih besar. Tingkat kebebasan dalam lmerTest::anovalebih berbeda dari anova sederhana, tetapi saya masih tidak yakin, mengapa mereka begitu besar untuk faktor / interaksi dalam subjek.

# mini example with ANT dataset from ez package
library(ez); library(lme4); library(lmerTest)

# repeated measures ANOVA with ez package
data(ANT)
ANT.2 <- subset(ANT, !error)
# update: make individual differences larger
baseline.shift <- rnorm(length(unique(ANT.2$subnum)), 0, 50)
ANT.2$rt <- ANT.2$rt + baseline.shift[as.numeric(ANT.2$subnum)]

anova.ez <- ezANOVA(data = ANT.2, dv = .(rt), wid = .(subnum), 
  within = .(direction), between = .(group))
anova.ez

# similarly with lmer and lmerTest::anova
model <- lmer(rt ~ group * direction + (1 | subnum), data = ANT.2)
lmerTest::anova(model)

# simple ANOVA on all available data
m <- lm(rt ~ group * direction, data = ANT.2)
anova(m)

Hasil dari kode di atas [ diperbarui ]:

anova.ez

$ ANOVA

           Effect DFn DFd         F          p p<.05          ges
2           group   1  18 2.6854464 0.11862957       0.1294475137
3       direction   1  18 0.9160571 0.35119193       0.0001690471
4 group:direction   1  18 4.9169156 0.03970473     * 0.0009066868

lmerTest :: anova (model)

Analysis of Variance Table of type 3  with  Satterthwaite 
approximation for degrees of freedom
                Df Sum Sq Mean Sq F value Denom Pr(>F)
group            1  13293   13293  2.6830    18 0.1188
direction        1   1946    1946  0.3935  5169 0.5305
group:direction  1  11563   11563  2.3321  5169 0.1268

anova (m)

Analysis of Variance Table

Response: rt
                  Df   Sum Sq Mean Sq  F value Pr(>F)    
group              1  1791568 1791568 242.3094 <2e-16 ***
direction          1      728     728   0.0985 0.7537    
group:direction    1    12024   12024   1.6262 0.2023    
Residuals       5187 38351225    7394                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Saya berpikir bahwa lmerTestsemakin benar dan ezanovasemakin salah dalam hal ini.

• hasil dari lmerTestsetuju dengan intuisi / pemahaman saya
• dua perhitungan berbeda dalam lmerTest(Satterthwaite dan Kenward-Roger) setuju
• mereka juga setuju nlme::lme
• ketika saya menjalankannya, ezanovamemberikan peringatan, yang saya tidak sepenuhnya mengerti, tetapi yang tidak boleh diabaikan ...

Menjalankan kembali contoh:

library(ez); library(lmerTest); library(nlme)
data(ANT)
ANT.2 <- subset(ANT, !error)
set.seed(101)  ## for reproducibility
baseline.shift <- rnorm(length(unique(ANT.2$subnum)), 0, 50)
ANT.2$rt <- ANT.2$rt + baseline.shift[as.numeric(ANT.2$subnum)]

Mencari tahu desain eksperimental

with(ANT.2,table(subnum,group,direction))

Jadi sepertinya individu ( subnum) ditempatkan dalam kelompok kontrol atau perlakuan, dan masing-masing diuji untuk kedua arah - yaitu arah dapat diuji dalam individu (penyebut df besar), tetapi kelompok dan kelompok: arah hanya dapat diuji di antara individu

(anova.ez <- ezANOVA(data = ANT.2, dv = .(rt), wid = .(subnum), 
    within = .(direction), between = .(group)))
## $ANOVA
##            Effect DFn DFd         F          p p<.05          ges
## 2           group   1  18 2.4290721 0.13651174       0.1183150147
## 3       direction   1  18 0.9160571 0.35119193       0.0002852171
## 4 group:direction   1  18 4.9169156 0.03970473     * 0.0015289914

Di sini saya dapatkan Warning: collapsing data to cell means. *IF* the requested effects are a subset of the full design, you must use the "within_full" argument, else results may be inaccurate. . Denominator DF terlihat sedikit funky (semua sama dengan 18): Saya pikir mereka harus lebih besar untuk arah dan grup: arah, yang dapat diuji secara independen (tetapi akan lebih kecil jika Anda ditambahkan (direction|subnum)ke model)?

# similarly with lmer and lmerTest::anova
model <- lmer(rt ~ group * direction + (1 | subnum), data = ANT.2)
lmerTest::anova(model)
##                 Df  Sum Sq Mean Sq F value Denom Pr(>F)
## group            1 12065.7 12065.7  2.4310    18 0.1364
## direction        1  1952.2  1952.2  0.3948  5169 0.5298
## group:direction  1 11552.2 11552.2  2.3299  5169 0.1270

yang Dfkolom di sini merujuk pada df pembilang, Denom(kedua-untuk-terakhir) memberikan perkiraan denominator df tersebut; mereka setuju dengan intuisi klasik. Lebih penting lagi, kami juga mendapatkan jawaban berbeda untuk nilai F ...

Kami juga dapat memeriksa ulang dengan Kenward-Roger ( sangat lambat karena melibatkan pengubahan model beberapa kali)

lmerTest::anova(model,ddf="Kenward-Roger")

Hasilnya identik.

Untuk contoh ini lme(dari nlmepaket) sebenarnya melakukan pekerjaan yang sangat baik dengan menebak penyebut yang sesuai df (nilai F dan p sangat sedikit berbeda):

model3 <- lme(rt ~ group * direction, random=~1|subnum, data = ANT.2)
anova(model3)[-1,]
##                 numDF denDF   F-value p-value
## group               1    18 2.4334314  0.1362
## direction           1  5169 0.3937316  0.5304
## group:direction     1  5169 2.3298847  0.1270

Jika saya cocok dengan interaksi antara directiondan subnumdf untuk directiondan group:directionjauh lebih kecil (saya akan berpikir mereka akan berusia 18, tapi mungkin saya mendapatkan sesuatu yang salah):

model2 <- lmer(rt ~ group * direction + (direction | subnum), data = ANT.2)
lmerTest::anova(model2)
##                 Df  Sum Sq Mean Sq F value   Denom Pr(>F)
## group            1 20334.7 20334.7  2.4302  17.995 0.1364
## direction        1  1804.3  1804.3  0.3649 124.784 0.5469
## group:direction  1 10616.6 10616.6  2.1418 124.784 0.1459

Saya biasanya setuju dengan analisis Ben, tetapi izinkan saya menambahkan beberapa komentar dan sedikit intuisi.

Pertama, hasil keseluruhan:

1. Hasil lmerTest menggunakan metode Satterthwaite sudah benar
2. Metode Kenward-Roger juga benar dan setuju dengan Satterthwaite

Ben menguraikan desain yang subnumbersarang groupsementara direction dan group:directiondisilangkan subnum. Ini berarti bahwa istilah kesalahan alami (yaitu yang disebut "melampirkan kesalahan strata") untuk groupadalah subnumsedangkan strata kesalahan melampirkan untuk istilah lain (termasuk subnum) adalah residual.

Struktur ini dapat direpresentasikan dalam diagram struktur-faktor:

names <- c(expression("[I]"[5169]^{5191}),
           expression("[subnum]"[18]^{20}), expression(grp:dir[1]^{4}),
           expression(dir[1]^{2}), expression(grp[1]^{2}), expression(0[1]^{1}))
x <- c(2, 4, 4, 6, 6, 8)
y <- c(5, 7, 5, 3, 7, 5)
plot(NA, NA, xlim=c(2, 8), ylim=c(2, 8), type="n", axes=F, xlab="", ylab="")
text(x, y, names) # Add text according to ’names’ vector
# Define coordinates for start (x0, y0) and end (x1, y1) of arrows:
x0 <- c(1.8, 1.8, 4.2, 4.2, 4.2, 6, 6) + .5
y0 <- c(5, 5, 7, 5, 5, 3, 7)
x1 <- c(2.7, 2.7, 5, 5, 5, 7.2, 7.2) + .5
y1 <- c(5, 7, 7, 3, 7, 5, 5)
arrows(x0, y0, x1, y1, length=0.1)

Di sini istilah acak terlampir dalam tanda kurung, 0mewakili keseluruhan rata-rata (atau mencegat), [I]mewakili istilah kesalahan, angka super-skrip adalah jumlah level dan nomor sub-skrip adalah jumlah derajat kebebasan dengan asumsi desain seimbang. Diagram menunjukkan bahwa istilah kesalahan alami (melampirkan strata kesalahan) untuk groupadalah subnumdan pembilang df untuk subnum, yang sama dengan penyebut df untuk group, adalah 18: 20 dikurangi 1 df untuk groupdan 1 df untuk rata-rata keseluruhan. Pengantar yang lebih komprehensif untuk diagram struktur faktor tersedia di bab 2 di sini: https://02429.compute.dtu.dk/eBook .

Jika data benar-benar seimbang, kami akan dapat membuat uji-F dari dekomposisi SSQ seperti yang disediakan oleh anova.lm. Karena dataset sangat seimbang, kami dapat memperoleh perkiraan F-tes sebagai berikut:

ANT.2 <- subset(ANT, !error)
set.seed(101)
baseline.shift <- rnorm(length(unique(ANT.2$subnum)), 0, 50)
ANT.2$rt <- ANT.2$rt + baseline.shift[as.numeric(ANT.2$subnum)]
fm <- lm(rt ~ group * direction + subnum, data=ANT.2)
(an <- anova(fm))
Analysis of Variance Table

Response: rt
                  Df   Sum Sq Mean Sq  F value Pr(>F)    
group              1   994365  994365 200.5461 <2e-16 ***
direction          1     1568    1568   0.3163 0.5739    
subnum            18  7576606  420923  84.8927 <2e-16 ***
group:direction    1    11561   11561   2.3316 0.1268    
Residuals       5169 25629383    4958                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Di sini semua nilai F dan p dihitung dengan asumsi bahwa semua istilah memiliki residual sebagai strata kesalahan terlampirnya, dan itu berlaku untuk semua kecuali 'grup'. Sebaliknya, F -test 'seimbang-benar' untuk grup adalah:

F_group <- an["group", "Mean Sq"] / an["subnum", "Mean Sq"]
c(Fvalue=F_group, pvalue=pf(F_group, 1, 18, lower.tail = FALSE))
   Fvalue    pvalue 
2.3623466 0.1416875 

di mana kita menggunakan subnumMS sebagai ganti ResidualsMS dalam penyebut nilai- F .

Perhatikan bahwa nilai-nilai ini sangat cocok dengan hasil Satterthwaite:

model <- lmer(rt ~ group * direction + (1 | subnum), data = ANT.2)
anova(model, type=1)
Type I Analysis of Variance Table with Satterthwaite's method
                 Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
group           12065.3 12065.3     1    18  2.4334 0.1362
direction        1951.8  1951.8     1  5169  0.3936 0.5304
group:direction 11552.2 11552.2     1  5169  2.3299 0.1270

Perbedaan yang tersisa disebabkan oleh data yang tidak seimbang.

OP membandingkan anova.lmdengan anova.lmerModLmerTest, yang ok, tetapi untuk membandingkan suka dengan seperti kita harus menggunakan kontras yang sama. Dalam hal ini ada perbedaan antara anova.lmdan anova.lmerModLmerTestkarena mereka menghasilkan tes Tipe I dan III secara default masing-masing, dan untuk dataset ini terdapat perbedaan (kecil) antara kontras Tipe I dan III:

show_tests(anova(model, type=1))$group
               (Intercept) groupTreatment directionright groupTreatment:directionright
groupTreatment           0              1    0.005202759                     0.5013477

show_tests(anova(model, type=3))$group # type=3 is default
               (Intercept) groupTreatment directionright groupTreatment:directionright
groupTreatment           0              1              0                           0.5

Jika set data benar-benar seimbang, kontras tipe I akan sama dengan kontras tipe III (yang tidak terpengaruh oleh jumlah sampel yang diamati).

Satu komentar terakhir adalah bahwa 'kelambatan' dari metode Kenward-Roger bukan karena model pemasangan kembali, tetapi karena melibatkan perhitungan dengan matriks varians-kovarians marginal dari pengamatan / residu (5191x5191 dalam kasus ini) yang tidak kasus untuk metode Satterthwaite.

Mengenai model2

Adapun model2 situasinya menjadi lebih kompleks dan saya pikir lebih mudah untuk memulai diskusi dengan model lain di mana saya telah memasukkan interaksi 'klasik' antara subnumdan direction:

model3 <- lmer(rt ~ group * direction + (1 | subnum) +
                 (1 | subnum:direction), data = ANT.2)
VarCorr(model3)
 Groups           Name        Std.Dev.  
 subnum:direction (Intercept) 1.7008e-06
 subnum           (Intercept) 4.0100e+01
 Residual                     7.0415e+01

Karena varians yang terkait dengan interaksi pada dasarnya nol (dengan adanya subnumefek utama acak), istilah interaksi tidak berpengaruh pada perhitungan derajat kebebasan, nilai- F dan nilai- p :

anova(model3, type=1)
Type I Analysis of Variance Table with Satterthwaite's method
                 Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
group           12065.3 12065.3     1    18  2.4334 0.1362
direction        1951.8  1951.8     1  5169  0.3936 0.5304
group:direction 11552.2 11552.2     1  5169  2.3299 0.1270

Namun, subnum:directionadalah strata kesalahan melampirkan untuk subnumjadi jika kita menghapus subnumsemua SSQ terkait kembali kesubnum:direction

model4 <- lmer(rt ~ group * direction +
                 (1 | subnum:direction), data = ANT.2)

Sekarang istilah kesalahan alami untuk group, directiondan group:directionadalah subnum:directiondan dengan nlevels(with(ANT.2, subnum:direction))= 40 dan empat parameter, derajat kebebasan penyebut untuk istilah-istilah tersebut seharusnya sekitar 36:

anova(model4, type=1)
Type I Analysis of Variance Table with Satterthwaite's method
                 Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
group           24004.5 24004.5     1 35.994  4.8325 0.03444 *
direction          50.6    50.6     1 35.994  0.0102 0.92020  
group:direction   273.4   273.4     1 35.994  0.0551 0.81583  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Ini F -tests juga dapat didekati dengan 'seimbang-benar' F -tests:

an4 <- anova(lm(rt ~ group*direction + subnum:direction, data=ANT.2))
an4[1:3, "F value"] <- an4[1:3, "Mean Sq"] / an4[4, "Mean Sq"]
an4[1:3, "Pr(>F)"] <- pf(an4[1:3, "F value"], 1, 36, lower.tail = FALSE)
an4
Analysis of Variance Table

Response: rt
                   Df   Sum Sq Mean Sq F value Pr(>F)    
group               1   994365  994365  4.6976 0.0369 *  
direction           1     1568    1568  0.0074 0.9319    
group:direction     1    10795   10795  0.0510 0.8226    
direction:subnum   36  7620271  211674 42.6137 <2e-16 ***
Residuals        5151 25586484    4967                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

sekarang beralih ke model2:

model2 <- lmer(rt ~ group * direction + (direction | subnum), data = ANT.2)

Model ini menggambarkan struktur kovarians efek-acak yang agak rumit dengan matriks varians-kovarians 2x2. Parameterisasi default tidak mudah untuk ditangani dan kami lebih baik dengan parameterisasi ulang model:

model2 <- lmer(rt ~ group * direction + (0 + direction | subnum), data = ANT.2)

Jika kita membandingkan model2dengan model4, mereka memiliki sama banyak acak-efek; 2 untuk masing-masing subnum, yaitu 2 * 20 = 40 total. Sementara model4menetapkan parameter varians tunggal untuk semua 40 efek acak, model2menetapkan bahwa setiap subnumpasangan efek acak memiliki distribusi normal bi-variate dengan matriks varians-kovarians 2x2 yang parameternya diberikan oleh

VarCorr(model2)
 Groups   Name           Std.Dev. Corr 
 subnum   directionleft  38.880        
          directionright 41.324   1.000
 Residual                70.405        

Ini menunjukkan terlalu pas, tapi mari kita simpan itu untuk hari lain. Yang penting di sini adalah bahwa model4adalah kasus khusus model2 dan yang modeladalah juga kasus khusus dari model2. Berbicara secara longgar (dan secara intuitif) (direction | subnum)mengandung atau menangkap variasi yang terkait dengan efek utama subnum serta interaksi direction:subnum. Dalam hal efek acak kita dapat menganggap dua efek atau struktur ini sebagai menangkap variasi antara baris dan baris-demi-kolom masing-masing:

head(ranef(model2)$subnum)
  directionleft directionright
1    -25.453576     -27.053697
2     16.446105      17.479977
3    -47.828568     -50.835277
4     -1.980433      -2.104932
5      5.647213       6.002221
6     41.493591      44.102056

Dalam hal ini, perkiraan efek acak ini serta estimasi parameter varians keduanya menunjukkan bahwa kami benar-benar hanya memiliki efek utama acak subnum(variasi antara baris) yang ada di sini. Apa yang menyebabkan semua ini adalah bahwa tingkat kebebasan Denominate Satterthwaite masuk

anova(model2, type=1)
Type I Analysis of Variance Table with Satterthwaite's method
                 Sum Sq Mean Sq NumDF   DenDF F value Pr(>F)
group           12059.8 12059.8     1  17.998  2.4329 0.1362
direction        1803.6  1803.6     1 125.135  0.3638 0.5475
group:direction 10616.6 10616.6     1 125.136  2.1418 0.1458

adalah kompromi antara efek-utama dan struktur interaksi ini: Kelompok DenDF tetap di 18 (bersarang subnumoleh desain) tetapi directiondan group:directionDenDF adalah kompromi antara 36 ( model4) dan 5169 ( model).

Saya tidak berpikir apa pun di sini menunjukkan bahwa perkiraan Satterthwaite (atau implementasinya di lmerTest ) salah.

Tabel setara dengan metode Kenward-Roger memberi

anova(model2, type=1, ddf="Ken")
Type I Analysis of Variance Table with Kenward-Roger's method
                 Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
group           12059.8 12059.8     1 18.000  2.4329 0.1362
direction        1803.2  1803.2     1 17.987  0.3638 0.5539
group:direction 10614.7 10614.7     1 17.987  2.1414 0.1606

Tidak mengherankan bahwa KR dan Satterthwaite dapat berbeda tetapi untuk semua tujuan praktis perbedaan dalam nilai- p adalah kecil. Analisis saya di atas menunjukkan bahwa DenDFuntuk directiondan group:directiontidak boleh lebih kecil dari ~ 36 dan mungkin lebih besar dari yang diberikan bahwa kita pada dasarnya hanya memiliki efek utama acak directionsaat ini, jadi jika saya pikir ini merupakan indikasi bahwa metode KR mendapatkan DenDFterlalu rendah pada kasus ini. Tetapi perlu diingat bahwa data tidak benar-benar mendukung (group | direction)struktur sehingga perbandingannya sedikit buatan - akan lebih menarik jika model tersebut benar-benar didukung.