Perbedaan antara PROC Mixed dan lme / lmer dalam R - derajat kebebasan

Catatan: pertanyaan ini adalah repost, karena pertanyaan saya sebelumnya harus dihapus karena alasan hukum.

Sambil membandingkan PROC CAMPURAN dari SAS dengan fungsi lmedari nlmepaket di R, saya menemukan beberapa perbedaan yang agak membingungkan. Lebih khusus lagi, derajat kebebasan dalam berbagai tes berbeda antara PROC MIXEDdan lme, dan saya bertanya-tanya mengapa.

Mulai dari dataset berikut (kode R yang diberikan di bawah):

• ind: faktor yang menunjukkan individu tempat pengukuran dilakukan
• fac: organ tempat pengukuran dilakukan
• trt: faktor yang mengindikasikan perawatan
• y: beberapa variabel respons kontinu

Idenya adalah untuk membangun model-model sederhana berikut:

y ~ trt + (ind): indsebagai faktor acak y ~ trt + (fac(ind)): facbersarang indsebagai faktor acak

Perhatikan bahwa model terakhir harus menyebabkan singularitas, karena hanya ada 1 nilai yuntuk setiap kombinasi inddan fac.

Model Pertama

Di SAS, saya membuat model berikut:

PROC MIXED data=Data;
    CLASS ind fac trt;
    MODEL y = trt /s;
    RANDOM ind /s;
run;

Menurut tutorial, model yang sama dalam R menggunakan nlmeharus:

> require(nlme)
> options(contrasts=c(factor="contr.SAS",ordered="contr.poly"))
> m2<-lme(y~trt,random=~1|ind,data=Data)

Kedua model memberikan estimasi yang sama untuk koefisien dan SE mereka, tetapi ketika melakukan uji F untuk efek trt, mereka menggunakan jumlah derajat kebebasan yang berbeda:

SAS : 
Type 3 Tests of Fixed Effects 
Effect Num DF Den DF     F  Value Pr > F 
trt         1      8  0.89        0.3724 

R : 
> anova(m2)
            numDF denDF  F-value p-value
(Intercept)     1     8 70.96836  <.0001
trt             1     6  0.89272  0.3812

Pertanyaan1: Apa perbedaan antara kedua tes? Keduanya dipasang menggunakan REML, dan menggunakan kontras yang sama.

CATATAN: Saya mencoba nilai yang berbeda untuk opsi DDFM = (termasuk BETWITHIN, yang secara teoritis akan memberikan hasil yang sama dengan lme)

Model Kedua

Dalam SAS:

PROC MIXED data=Data;
    CLASS ind fac trt;
    MODEL y = trt /s;
    RANDOM fac(ind) /s;
run;

Model setara dalam R harus:

> m4<-lme(y~trt,random=~1|ind/fac,data=Data)

Dalam hal ini, ada beberapa perbedaan yang sangat aneh:

• R cocok tanpa mengeluh, sedangkan SAS mencatat bahwa goni akhir tidak pasti positif (yang tidak sedikit mengejutkan saya, lihat di atas)
• SE pada koefisien berbeda (lebih kecil dalam SAS)
• Sekali lagi, uji F menggunakan jumlah DF yang berbeda (pada kenyataannya, dalam SAS jumlah itu = 0)

Output SAS:

Effect     trt Estimate Std Error  DF t Value Pr > |t| 
Intercept        0.8863    0.1192  14    7.43 <.0001 
trt       Cont  -0.1788    0.1686   0   -1.06 . 

R Output:

> summary(m4)
...
Fixed effects: y ~ trt 
               Value Std.Error DF   t-value p-value
(Intercept)  0.88625 0.1337743  8  6.624963  0.0002
trtCont     -0.17875 0.1891855  6 -0.944840  0.3812
...

(Perhatikan bahwa dalam kasus ini, uji F dan T setara dan gunakan DF yang sama.)

Menariknya, ketika menggunakan lme4dalam R model bahkan tidak cocok:

> require(lme4)
> m4r <- lmer(y~trt+(1|ind/fac),data=Data)
Error in function (fr, FL, start, REML, verbose)  : 
  Number of levels of a grouping factor for the random effects
must be less than the number of observations

Pertanyaan 2 : Apa perbedaan antara model ini dengan faktor bersarang? Apakah mereka ditentukan dengan benar dan jika demikian, bagaimana hasilnya sangat berbeda?

Data Simulasi dalam R:

Data <- structure(list(y = c(1.05, 0.86, 1.02, 1.14, 0.68, 1.05, 0.22, 
1.07, 0.46, 0.65, 0.41, 0.82, 0.6, 0.49, 0.68, 1.55), ind = structure(c(1L, 
2L, 3L, 1L, 3L, 4L, 4L, 2L, 5L, 6L, 7L, 8L, 6L, 5L, 7L, 8L), .Label = c("1", 
"2", "3", "4", "5", "6", "7", "8"), class = "factor"), fac = structure(c(1L, 
1L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 1L, 2L, 2L), .Label = c("l", 
"r"), class = "factor"), trt = structure(c(2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("Cont", 
"Treat"), class = "factor")), .Names = c("y", "ind", "fac", "trt"
), row.names = c(NA, -16L), class = "data.frame")

Data Simulasi:

   y ind fac   trt
1.05   1   l Treat
0.86   2   l Treat
1.02   3   l Treat
1.14   1   r Treat
0.68   3   r Treat
1.05   4   l Treat
0.22   4   r Treat
1.07   2   r Treat
0.46   5   r  Cont
0.65   6   l  Cont
0.41   7   l  Cont
0.82   8   l  Cont
0.60   6   r  Cont
0.49   5   l  Cont
0.68   7   r  Cont
1.55   8   r  Cont

Untuk pertanyaan pertama, metode default di SAS untuk menemukan df tidak terlalu pintar; ia mencari istilah dalam efek acak yang secara sintaksis menyertakan efek tetap, dan menggunakannya. Dalam hal ini, karena trttidak ditemukan ind, itu tidak melakukan hal yang benar. Saya belum pernah mencoba BETWITHINdan tidak tahu detailnya, tetapi baik opsi Satterthwaite ( satterth) atau menggunakan ind*trtsebagai efek acak memberikan hasil yang benar.

PROC MIXED data=Data;
    CLASS ind fac trt;
    MODEL y = trt /s ddfm=satterth;
    RANDOM ind /s;
run;

PROC MIXED data=Data;
    CLASS ind fac trt;
    MODEL y = trt /s;
    RANDOM ind*trt /s;
run;

Adapun pertanyaan kedua, kode SAS Anda tidak cukup cocok dengan kode R Anda; hanya memiliki istilah untuk fac*ind, sedangkan kode R memiliki istilah untuk keduanya inddan fac*ind. (Lihat output Variance Components untuk melihat ini.) Menambahkan ini memberikan SE yang sama untuk trtsemua model di Q1 dan Q2 (0.1892).

Seperti yang Anda perhatikan, ini adalah model aneh yang cocok karena fac*indistilah tersebut memiliki satu pengamatan untuk setiap level, sehingga setara dengan istilah kesalahan. Ini tercermin dalam output SAS, di mana fac*indistilah tersebut memiliki nol varians. Ini juga yang diberitahukan pesan kesalahan dari lme4 kepada Anda; alasan kesalahan adalah bahwa Anda kemungkinan besar salah menentukan sesuatu karena Anda memasukkan istilah kesalahan dalam model dalam dua cara yang berbeda. Menariknya, ada satu sedikit perbedaan dalam model nlme; itu entah bagaimana menemukan istilah varians untuk fac*indistilah di samping istilah kesalahan, tetapi Anda akan melihat bahwa jumlah kedua varians ini sama dengan istilah kesalahan dari SAS dan nlme tanpa fac*indistilah. Namun, SE untuk trttetap sama (0,1892)trt yang bersarang diind, jadi istilah varian yang lebih rendah ini tidak memengaruhinya.

Akhirnya, catatan umum tentang derajat kebebasan dalam model ini: Mereka dihitung setelah model cocok, dan perbedaan derajat kebebasan antara berbagai program atau opsi program tidak selalu berarti bahwa model tersebut sedang fit berbeda. Untuk itu, kita harus melihat pada perkiraan parameter, baik parameter efek tetap dan parameter kovarian.

Juga, menggunakan perkiraan t dan F dengan jumlah derajat kebebasan tertentu cukup kontroversial. Tidak hanya ada beberapa cara untuk memperkirakan df, beberapa percaya praktik melakukannya juga bukan ide yang baik. Beberapa kata nasihat:

1. Jika semuanya seimbang, bandingkan hasilnya dengan metode kuadrat terkecil tradisional, seperti yang seharusnya disepakati. Jika hampir seimbang, hitung sendiri (dengan asumsi keseimbangan) sehingga Anda dapat memastikan yang Anda gunakan berada di stadion baseball yang tepat.

2. Jika Anda memiliki ukuran sampel yang besar, tingkat kebebasan tidak terlalu penting karena distribusinya mendekati normal dan chi-square.

3. Periksa metode inferensi Doug Bates. Metode yang lebih lama didasarkan pada simulasi MCMC; metode yang lebih baru didasarkan pada profil kemungkinan.