Kesetaraan antara model anova tindakan berulang dan model campuran: lmer vs lme, dan simetri majemuk

Saya memiliki beberapa kesulitan untuk mendapatkan hasil yang setara antara aovmodel tindakan berulang antara-dalam dan lmermodel campuran.

Data dan skrip saya terlihat sebagai berikut

data=read.csv("https://www.dropbox.com/s/zgle45tpyv5t781/fitness.csv?dl=1")
data$id=factor(data$id)
data
   id  FITNESS      TEST PULSE
1   1  pilates   CYCLING    91
2   2  pilates   CYCLING    82
3   3  pilates   CYCLING    65
4   4  pilates   CYCLING    90
5   5  pilates   CYCLING    79
6   6  pilates   CYCLING    84
7   7 aerobics   CYCLING    84
8   8 aerobics   CYCLING    77
9   9 aerobics   CYCLING    71
10 10 aerobics   CYCLING    91
11 11 aerobics   CYCLING    72
12 12 aerobics   CYCLING    93
13 13    zumba   CYCLING    63
14 14    zumba   CYCLING    87
15 15    zumba   CYCLING    67
16 16    zumba   CYCLING    98
17 17    zumba   CYCLING    63
18 18    zumba   CYCLING    72
19  1  pilates   JOGGING   136
20  2  pilates   JOGGING   119
21  3  pilates   JOGGING   126
22  4  pilates   JOGGING   108
23  5  pilates   JOGGING   122
24  6  pilates   JOGGING   101
25  7 aerobics   JOGGING   116
26  8 aerobics   JOGGING   142
27  9 aerobics   JOGGING   137
28 10 aerobics   JOGGING   134
29 11 aerobics   JOGGING   131
30 12 aerobics   JOGGING   120
31 13    zumba   JOGGING    99
32 14    zumba   JOGGING    99
33 15    zumba   JOGGING    98
34 16    zumba   JOGGING    99
35 17    zumba   JOGGING    87
36 18    zumba   JOGGING    89
37  1  pilates SPRINTING   179
38  2  pilates SPRINTING   195
39  3  pilates SPRINTING   188
40  4  pilates SPRINTING   189
41  5  pilates SPRINTING   173
42  6  pilates SPRINTING   193
43  7 aerobics SPRINTING   184
44  8 aerobics SPRINTING   179
45  9 aerobics SPRINTING   179
46 10 aerobics SPRINTING   174
47 11 aerobics SPRINTING   164
48 12 aerobics SPRINTING   182
49 13    zumba SPRINTING   111
50 14    zumba SPRINTING   103
51 15    zumba SPRINTING   113
52 16    zumba SPRINTING   118
53 17    zumba SPRINTING   127
54 18    zumba SPRINTING   113

Pada dasarnya, 3 x 6 subjek ( id) menjadi sasaran tiga FITNESSskema latihan yang berbeda masing-masing dan mereka PULSEdiukur setelah melakukan tiga jenis ketahanan yang berbeda TEST.

Saya kemudian memasang aovmodel berikut :

library(afex)
library(car)
set_sum_contrasts()
fit1 = aov(PULSE ~ FITNESS*TEST + Error(id/TEST),data=data)
summary(fit1)
Error: id
          Df Sum Sq Mean Sq F value   Pr(>F)    
FITNESS    2  14194    7097   115.1 7.92e-10 ***
Residuals 15    925      62                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: id:TEST
             Df Sum Sq Mean Sq F value   Pr(>F)    
TEST          2  57459   28729   253.7  < 2e-16 ***
FITNESS:TEST  4   8200    2050    18.1 1.16e-07 ***
Residuals    30   3397     113                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Hasilnya saya dapatkan menggunakan

set_sum_contrasts()
fit2=aov.car(PULSE ~ FITNESS*TEST+Error(id/TEST),data=data,type=3,return="Anova")
summary(fit2)

identik dengan ini.

Model campuran dijalankan menggunakan nlmememberikan hasil yang setara langsung, misalnya menggunakan lme:

library(lmerTest)    
lme1=lme(PULSE ~ FITNESS*TEST, random=~1|id, correlation=corCompSymm(form=~1|id),data=data)
anova(lme1)
             numDF denDF   F-value p-value
(Intercept)      1    30 12136.126  <.0001
FITNESS          2    15   115.127  <.0001
TEST             2    30   253.694  <.0001
FITNESS:TEST     4    30    18.103  <.0001


summary(lme1)
Linear mixed-effects model fit by REML
 Data: data 
       AIC      BIC    logLik
  371.5375 393.2175 -173.7688

Random effects:
 Formula: ~1 | id
        (Intercept) Residual
StdDev:    1.699959 9.651662

Correlation Structure: Compound symmetry
 Formula: ~1 | id 
 Parameter estimate(s):
       Rho 
-0.2156615 
Fixed effects: PULSE ~ FITNESS * TEST 
                                 Value Std.Error DF   t-value p-value
(Intercept)                   81.33333  4.000926 30 20.328628  0.0000
FITNESSpilates                 0.50000  5.658164 15  0.088368  0.9308
FITNESSzumba                  -6.33333  5.658164 15 -1.119327  0.2806
TESTJOGGING                   48.66667  6.143952 30  7.921069  0.0000
TESTSPRINTING                 95.66667  6.143952 30 15.570868  0.0000
FITNESSpilates:TESTJOGGING   -11.83333  8.688861 30 -1.361897  0.1834
FITNESSzumba:TESTJOGGING     -28.50000  8.688861 30 -3.280062  0.0026
FITNESSpilates:TESTSPRINTING   8.66667  8.688861 30  0.997446  0.3265
FITNESSzumba:TESTSPRINTING   -56.50000  8.688861 30 -6.502579  0.0000

Atau menggunakan gls:

library(lmerTest)    
gls1=gls(PULSE ~ FITNESS*TEST, correlation=corCompSymm(form=~1|id),data=data)
anova(gls1)

Namun, hasil yang saya memperoleh menggunakan lme4's lmerberbeda:

set_sum_contrasts()
fit3=lmer(PULSE ~ FITNESS*TEST+(1|id),data=data)
summary(fit3)
Linear mixed model fit by REML ['lmerMod']
Formula: PULSE ~ FITNESS * TEST + (1 | id)
   Data: data

REML criterion at convergence: 362.4

Random effects:
 Groups   Name        Variance Std.Dev.
 id       (Intercept)  0.00    0.0     
 Residual             96.04    9.8     
...

Anova(fit3,test.statistic="F",type=3)
Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)

Response: PULSE
                    F Df Df.res    Pr(>F)    
(Intercept)  7789.360  1     15 < 2.2e-16 ***
FITNESS        73.892  2     15 1.712e-08 ***
TEST          299.127  2     30 < 2.2e-16 ***
FITNESS:TEST   21.345  4     30 2.030e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Adakah yang salah dengan lmermodel yang saya lakukan ? Atau dari mana perbedaan itu berasal? Mungkinkah ia melakukan sesuatu dengan lmertidak mengizinkan corelations intraclass negatif atau sesuatu seperti itu? Mengingat bahwa nlme's glsdan lmejangan mengembalikan hasil yang benar, meskipun, saya bertanya-tanya bagaimana ini berbeda dalam glsdan lme? Apakah pilihan itu correlation=corCompSymm(form=~1|id)menyebabkan mereka untuk secara langsung memperkirakan korelasi intraclass, yang bisa positif atau negatif, sedangkan lmermemperkirakan komponen varians, yang tidak boleh negatif (dan akhirnya diperkirakan nol dalam kasus ini)?

Seperti yang telah disebutkan oleh Ben Bolker dalam komentar, masalahnya adalah seperti yang Anda duga: lmer()Model menjadi tersandung karena mencoba untuk memperkirakan model komponen varians, dengan estimasi komponen varians dibatasi menjadi non-negatif. Apa yang akan saya coba lakukan adalah memberikan pemahaman yang agak intuitif tentang apa itu dataset Anda yang mengarah ke ini, dan mengapa ini menyebabkan masalah untuk model komponen varians.

Ini adalah plot dataset Anda. Titik-titik putih adalah pengamatan aktual dan titik-titik hitam adalah sarana subjek.

Untuk membuat segalanya lebih sederhana, tetapi tanpa mengubah semangat masalah, saya akan mengurangi efek tetap (yaitu, FITNESSdan TESTefek, serta grand mean) dan menangani data residual sebagai masalah efek acak satu arah . Jadi inilah tampilan dataset baru:

Perhatikan baik-baik pola di plot ini. Pikirkan tentang bagaimana pengamatan yang diambil dari subjek yang sama berbeda dari pengamatan yang diambil dari subyek yang berbeda. Secara khusus, perhatikan pola berikut: Karena salah satu pengamatan untuk subjek lebih tinggi (atau lebih rendah) di atas (atau di bawah) rata-rata subjek, pengamatan lain dari subjek tersebut cenderung berada di sisi berlawanan dari rata-rata subjek. Dan semakin jauh observasi itu dari mean subjek, semakin jauh observasi lain cenderung dari mean subjek di sisi yang berlawanan. Ini menunjukkan korelasi intra kelas yang negatif. Dua pengamatan yang diambil dari subjek yang sama sebenarnya cenderung kurang serupa, rata-rata, dari dua pengamatan yang diambil secara acak dari dataset.

Cara lain untuk berpikir tentang pola ini adalah dalam hal besaran relatif antara varians antara subjek dan dalam subjek. Tampaknya ada varians dalam-subjek yang jauh lebih besar dibandingkan dengan varians antar-subjek. Tentu saja, kami berharap ini terjadi sampai batas tertentu. Setelah semua, varians dalam subjek didasarkan pada variasi dalam poin data individu, sedangkan varians antara subjek didasarkan pada variasi dalam cara poin data individu (yaitu, berarti subjek), dan kita tahu bahwa varians dari rata-rata akan cenderung menurun karena jumlah hal yang dirata-rata meningkat. Namun dalam dataset ini perbedaannya cukup mencolok: Ada caralebih banyak di dalam subjek daripada di antara variasi subjek. Sebenarnya perbedaan ini adalah alasan mengapa korelasi intra-kelas negatif muncul.

Oke, jadi inilah masalahnya. Model komponen varians mengasumsikan bahwa setiap titik data adalah jumlah dari efek subjek dan kesalahan:$y_{ij}=u_j+e_{ij}$dimana $u_j$ adalah efek dari $j$subjek ke-5. Jadi mari kita berpikir tentang apa yang akan terjadi jika benar-benar ada 0 varians dalam efek subjek - dengan kata lain, jika komponen varians antar-subyek yang benar adalah 0. Dengan dataset aktual yang dihasilkan dalam model ini, jika kita menghitung sampel berarti untuk setiap data yang diamati setiap subjek, sampel tersebut berarti masih memiliki beberapa varian tidak nol, tetapi mereka hanya mencerminkan varians kesalahan, dan bukan varians subjek "benar" (karena kami berasumsi tidak ada varians).

Jadi bagaimana variabel yang kita harapkan dari subjek ini? Nah, pada dasarnya setiap perkiraan efek subjek adalah rata-rata, dan kami tahu rumus untuk varians rata-rata:$\text{var}(\bar{X})=\text{var}(X_i)/n$dimana $n$adalah jumlah hal yang dirata-rata. Sekarang mari kita terapkan formula ini untuk dataset Anda dan melihat berapa banyak varians kita akan mengharapkan untuk melihat efek subjek diperkirakan jika benar antara-subyek komponen varians persis 0.

Varians dalam subjek berhasil $348$, dan setiap efek subjek dihitung sebagai rata-rata dari 3 pengamatan. Jadi deviasi standar yang diharapkan dalam subjek berarti - dengan asumsi varians antara-subjek yang sebenarnya adalah 0 - berhasil menjadi sekitar$10.8$. Sekarang bandingkan ini dengan standar deviasi pada subjek yang berarti kita benar-benar mengamati:$4.3$! Variasi yang diamati secara substansial kurang dari variasi yang diharapkan ketika kami mengasumsikan 0 antar-subjek varians. Untuk model komponen varians, satu-satunya cara bahwa variasi yang diamati dapat diharapkan serendah apa yang sebenarnya kami amati adalah jika varians antar-subjek yang sebenarnya entah bagaimana negatif . Dan di situlah letak masalahnya. Data menyiratkan bahwa ada entah bagaimana komponen varians negatif, tetapi perangkat lunak (masuk akal) tidak akan memungkinkan perkiraan negatif komponen varians, karena varians sebenarnya tidak pernah bisa menjadi negatif. Model-model lain yang Anda cocok menghindari masalah ini dengan langsung memperkirakan korelasi intra-kelas daripada mengasumsikan model komponen varians sederhana.

Jika Anda ingin melihat bagaimana Anda benar-benar bisa mendapatkan estimasi komponen varians negatif yang tersirat oleh dataset Anda, Anda dapat menggunakan prosedur yang saya ilustrasikan (dengan Rkode yang menyertainya ) dalam jawaban saya terbaru ini . Prosedur itu tidak sepenuhnya sepele, tetapi juga tidak terlalu sulit (untuk desain yang seimbang seperti ini).