Bagaimana memilih struktur efek acak dan tetap dalam model linier campuran?

Pertimbangkan data berikut dari dua arah dalam desain subjek:

df <- "http://personality-project.org/r/datasets/R.appendix4.data"
df <- read.table(df,header=T)
head(df)

Observation Subject Task Valence Recall
1           1     Jim Free     Neg      8
2           2     Jim Free     Neu      9
3           3     Jim Free     Pos      5
4           4     Jim Cued     Neg      7
5           5     Jim Cued     Neu      9
6           6     Jim Cued     Pos     10

Saya ingin menganalisis ini menggunakan model linier campuran. Mempertimbangkan semua kemungkinan efek tetap dan acak ada beberapa model yang mungkin:

# different fixed effects with random-intercept
a0 <- lmer(Recall~1 + (1|Subject), REML=F,df)
a1 <- lmer(Recall~Task + (1|Subject), REML=F,df)
a2 <- lmer(Recall~Valence + (1|Subject), REML=F,df)
a3 <- lmer(Recall~Task+Valence + (1|Subject), REML=F,df)
a4 <- lmer(Recall~Task*Valence + (1|Subject), REML=F,df)

# different fixed effects with random-intercept-random-slope
b0 <- lmer(Recall~1 + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)
b1 <- lmer(Recall~Task + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)
b2 <- lmer(Recall~Valence + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)
b3 <- lmer(Recall~Task+Valence + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)
b4 <- lmer(Recall~Task*Valence + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)

# different fixed effects with random-intercept-random-slope including variance-covariance matrix
c0 <- lmer(Recall~1 + (1 + Valence + Task|Subject), REML=F,df)
c1 <- lmer(Recall~Task + (1 + Valence + Task|Subject), REML=F,df)
c2 <- lmer(Recall~Valence + (1 + Valence + Task|Subject), REML=F,df)
c3 <- lmer(Recall~Task+Valence + (1 + Valence + Task|Subject), REML=F,df)
c4 <- lmer(Recall~Task*Valence + (1 + Valence + Task|Subject), REML=F,df)

1. Apa cara yang disarankan untuk memilih model pemasangan terbaik dalam konteks ini? Saat menggunakan tes rasio log-likelihood apa prosedur yang direkomendasikan? Membuat model ke atas (dari model nol ke model paling kompleks) atau ke bawah (dari model paling kompleks ke model nol)? Inklusi atau pengecualian bertahap? Atau direkomendasikan untuk menempatkan semua model dalam satu uji rasio log-likelihood dan memilih model dengan nilai p terendah? Bagaimana membandingkan model yang tidak bersarang?

2. Apakah disarankan untuk pertama-tama menemukan struktur efek-efek yang sesuai dan kemudian struktur efek-acak yang sesuai atau sebaliknya (saya telah menemukan referensi untuk kedua opsi ...)?

3. Apa cara yang disarankan untuk melaporkan hasil? Melaporkan nilai-p dari uji rasio log-likelihood yang membandingkan model campuran lengkap (dengan efek yang dipertanyakan) dengan model yang dikurangi (tanpa efek yang dipertanyakan). Atau lebih baik menggunakan uji rasio log-likelihood untuk menemukan model pemasangan terbaik dan kemudian menggunakan lmerTest untuk melaporkan nilai-p dari efek dalam model pemasangan terbaik?

Saya tidak yakin benar-benar ada jawaban kanonik untuk ini, tetapi saya akan mencobanya.

Apa cara yang disarankan untuk memilih model pemasangan terbaik dalam konteks ini? Saat menggunakan tes rasio log-likelihood apa prosedur yang direkomendasikan? Membuat model ke atas (dari model nol ke model paling kompleks) atau ke bawah (dari model paling kompleks ke model nol)? Inklusi atau pengecualian bertahap? Atau direkomendasikan untuk menempatkan semua model dalam satu uji rasio log-likelihood dan memilih model dengan nilai p terendah? Bagaimana membandingkan model yang tidak bersarang?

Itu tergantung apa tujuan Anda.

• Secara umum Anda harus sangat , sangat berhati-hati dalam pemilihan model (lihat misalnya jawaban ini , atau posting ini , atau hanya Google "Harrell stepwise" ...).
• Jika Anda tertarik untuk memiliki nilai-p Anda menjadi bermakna (yaitu Anda melakukan pengujian hipotesis konfirmasi), Anda tidak harus melakukan pemilihan model. Namun : tidak begitu jelas bagi saya apakah prosedur pemilihan model cukup buruk jika Anda melakukan pemilihan model pada bagian non-fokus dari model , misalnya melakukan pemilihan model pada efek acak jika minat utama Anda adalah kesimpulan tentang efek tetap.
• Tidak ada yang namanya "menempatkan semua model dalam satu uji rasio kemungkinan" - pengujian rasio kemungkinan adalah prosedur berpasangan. Jika Anda ingin melakukan pemilihan model (misalnya) pada efek acak, saya mungkin akan merekomendasikan pendekatan "sekaligus" menggunakan kriteria informasi seperti dalam contoh ini - yang setidaknya menghindari beberapa masalah pendekatan bertahap (tapi bukan dari pemilihan model lebih umum).
• Barr et al. 2013 "Keep it maximal" Jurnal Memori dan Bahasa (doi: 10.1016 / j.jml.2012.11.001) akan merekomendasikan menggunakan model maksimal (hanya).
• Shravan Vasishth tidak setuju , dengan alasan bahwa model seperti itu akan kurang bertenaga dan karenanya bermasalah kecuali kumpulan data sangat besar (dan rasio sinyal-ke-noise tinggi)
• Pendekatan lain yang cukup masuk akal adalah mencocokkan model yang besar tetapi masuk akal dan kemudian, jika cocok adalah tunggal, hapus persyaratan sampai tidak lagi
• Dengan beberapa peringatan (disebutkan dalam FAQ GLMM ), Anda dapat menggunakan kriteria informasi untuk membandingkan model yang tidak bersarang dengan berbagai efek acak (walaupun Brian Ripley tidak setuju: lihat bagian bawah hal. 6 di sini )

Apakah disarankan untuk pertama-tama menemukan struktur efek-efek yang sesuai dan kemudian struktur efek-acak yang sesuai atau sebaliknya (saya telah menemukan referensi untuk kedua opsi ...)?

Saya tidak berpikir ada yang tahu. Lihat jawaban sebelumnya tentang pemilihan model secara lebih umum. Jika Anda dapat mendefinisikan tujuan Anda dengan cukup jelas (yang dilakukan beberapa orang), pertanyaannya mungkin dapat dijawab. Jika Anda memiliki referensi untuk kedua opsi, akan berguna untuk mengedit pertanyaan Anda untuk memasukkannya ... (Untuk apa nilainya, contoh ini (sudah dikutip di atas) menggunakan kriteria informasi untuk memilih bagian efek acak, kemudian menghindari seleksi pada bagian efek tetap dari model.

Apa cara yang disarankan untuk melaporkan hasil? Melaporkan nilai-p dari uji rasio log-likelihood yang membandingkan model campuran lengkap (dengan efek yang dipertanyakan) dengan model yang dikurangi (tanpa efek yang dipertanyakan). Atau lebih baik menggunakan uji rasio log-likelihood untuk menemukan model pemasangan terbaik dan kemudian menggunakan lmerTest untuk melaporkan nilai-p dari efek dalam model pemasangan terbaik?

Ini (sayangnya) pertanyaan sulit lainnya. Jika Anda melaporkan efek marginal seperti dilansir lmerTest, Anda harus khawatir tentang keterpinggiran (misalnya, apakah perkiraan efek utama Adan Bbermakna bila ada A-by- Binteraksi dalam model); ini adalah kaleng besar cacing, tetapi agak dikurangi jika Anda menggunakan contrasts="sum"seperti yang disarankan oleh afex::mixed(). Desain seimbang membantu sedikit juga. Jika Anda benar-benar ingin mengatasi semua retakan ini, saya pikir saya akan merekomendasikan afex::mixed, yang memberi Anda hasil yang serupa lmerTest, tetapi mencoba menangani masalah ini.

Pembaruan Mei 2017 : Ternyata, sebagian dari apa yang saya tulis di sini agak salah . Beberapa pembaruan dilakukan di seluruh pos.

Saya setuju banyak dengan apa yang telah dikatakan oleh Ben Bolker (terima kasih atas teriakannya afex::mixed()), tetapi izinkan saya menambahkan beberapa pemikiran yang lebih umum dan spesifik tentang masalah ini.

Fokus pada efek tetap versus acak dan cara melaporkan hasil

Untuk jenis penelitian eksperimental yang direpresentasikan dalam contoh kumpulan data dari Jonathan Baron Anda menggunakan pertanyaan penting biasanya apakah faktor yang dimanipulasi memiliki efek keseluruhan. Sebagai contoh, apakah kita menemukan efek utama atau interaksi keseluruhan Task? Poin penting adalah bahwa dalam set data biasanya semua faktor berada di bawah kontrol eksperimental lengkap dan ditetapkan secara acak. Akibatnya, fokus bunga biasanya pada efek tetap.
Sebaliknya, komponen efek acak dapat dilihat sebagai parameter "gangguan" yang menangkap varians sistematis (yaitu, perbedaan antar individu dalam ukuran efek) yang tidak selalu penting untuk pertanyaan utama. Dari sudut pandang ini saran untuk menggunakan struktur efek acak maksimal seperti yang dianjurkan oleh Barr et al. mengikuti agak alami. Mudah dibayangkan bahwa manipulasi eksperimental tidak memengaruhi semua individu dengan cara yang persis sama dan kami ingin mengendalikannya. Di sisi lain, jumlah faktor atau level biasanya tidak terlalu besar sehingga bahaya overfitting tampaknya relatif kecil.

Akibatnya, saya akan mengikuti saran dari Barr et al. dan tentukan struktur efek acak maksimal dan laporkan pengujian efek tetap sebagai hasil utama saya. Untuk menguji efek tetap, saya juga menyarankan untuk menggunakannya afex::mixed()karena melaporkan pengujian efek atau faktor (alih-alih uji parameter) dan menghitung tes tersebut dengan cara yang agak masuk akal (misalnya, menggunakan struktur efek acak yang sama untuk semua model di mana suatu efek tunggal dihapus, menggunakan jumlah-ke-nol-kontras, menawarkan metode yang berbeda untuk menghitung nilai- p , ...).

Bagaimana dengan contoh data

Masalah dengan contoh data yang Anda berikan adalah bahwa untuk dataset ini struktur efek acak maksimal mengarah ke model jenuh karena hanya ada satu titik data per sel desain:

> with(df, table(Valence, Subject, Task))
, , Task = Cued

       Subject
Valence Faye Jason Jim Ron Victor
    Neg    1     1   1   1      1
    Neu    1     1   1   1      1
    Pos    1     1   1   1      1

, , Task = Free

       Subject
Valence Faye Jason Jim Ron Victor
    Neg    1     1   1   1      1
    Neu    1     1   1   1      1
    Pos    1     1   1   1      1

Akibatnya, lmertersedak pada struktur efek acak maksimal:

> lmer(Recall~Task*Valence + (Valence*Task|Subject), df)
Error: number of observations (=30) <= number of random effects (=30) for term
(Valence * Task | Subject); the random-effects parameters and the residual variance
(or scale parameter) are probably unidentifiable

Sayangnya, setahu saya tidak ada cara yang disepakati untuk menangani masalah ini. Tapi izinkan saya membuat sketsa dan membahas beberapa:

1. Solusi pertama bisa dengan menghapus kemiringan acak tertinggi dan menguji efek untuk model ini:

   require(afex)
   mixed(Recall~Task*Valence + (Valence+Task|Subject), df)
           Effect    F ndf  ddf F.scaling p.value
   1         Task 6.56   1 4.00      1.00     .06
   2      Valence 0.80   2 3.00      0.75     .53
   3 Task:Valence 0.42   2 8.00      1.00     .67

   Namun, solusi ini sedikit ad-hoc dan tidak terlalu termotivasi.

   Pembaruan Mei 2017: Ini adalah pendekatan yang saya setujui saat ini. Lihat posting blog ini dan draf bab saya penulis bersama , bagian "Struktur Efek Acak untuk Desain ANOVA Tradisional".

2. Solusi alternatif (dan yang bisa dilihat sebagai advokasi oleh diskusi Barr et al.) Adalah untuk selalu menghilangkan lereng acak untuk efek terkecil. Ini memiliki dua masalah: (1) Struktur efek acak mana yang kita gunakan untuk mencari tahu apa efek terkecil dan (2) R enggan untuk menghapus efek orde rendah seperti efek utama jika efek orde tinggi seperti interaksi efek ini hadir (lihat di sini ). Sebagai akibatnya, seseorang perlu mengatur struktur efek acak ini dengan tangan dan meneruskan matriks model yang dikonstruksi ke lmer call.

3. Solusi ketiga bisa menggunakan parametrization alternatif dari bagian efek acak, yaitu yang sesuai dengan model RM-ANOVA untuk data ini. Sayangnya (?), lmerTidak mengizinkan "varians negatif" sehingga parameterisasi ini tidak persis sesuai dengan RM-ANOVA untuk semua set data , lihat diskusi di sini dan di tempat lain (misalnya di sini dan di sini ). "Lmer-ANOVA" untuk data ini adalah:

   > mixed(Recall~Task*Valence + (1|Subject) + (1|Task:Subject) + (1|Valence:Subject), df)
           Effect    F ndf  ddf F.scaling p.value
   1         Task 7.35   1 4.00      1.00     .05
   2      Valence 1.46   2 8.00      1.00     .29
   3 Task:Valence 0.29   2 8.00      1.00     .76

Mengingat semua masalah ini, saya tidak akan menggunakan lmerset data yang hanya memiliki satu titik data per sel desain kecuali ada solusi yang lebih disepakati untuk masalah struktur efek acak maksimal yang tersedia.

4. Sebaliknya, saya akan One juga masih bisa menggunakan ANOVA klasik. Menggunakan salah satu pembungkus untuk car::Anova()dalam afexhasilnya adalah:

   > aov4(Recall~Task*Valence + (Valence*Task|Subject), df)
           Effect         df  MSE      F  ges   p
   1      Valence 1.44, 5.75 4.67   1.46  .02 .29
   2         Task       1, 4 4.08 7.35 +  .07 .05
   3 Valence:Task 1.63, 6.52 2.96   0.29 .003 .71

   Perhatikan bahwa afexsekarang juga memungkinkan untuk mengembalikan model yang cocok dengan aovyang dapat diteruskan ke lsmeansuntuk tes post-hoc (tetapi untuk uji efek yang dilaporkan car::Anovamasih lebih masuk akal):

   > require(lsmeans)
   > m <- aov4(Recall~Task*Valence + (Valence*Task|Subject), df, return = "aov")
   > lsmeans(m, ~Task+Valence)
    Task Valence lsmean       SE   df lower.CL upper.CL
    Cued Neg       11.8 1.852026 5.52  7.17157 16.42843
    Free Neg       10.2 1.852026 5.52  5.57157 14.82843
    Cued Neu       13.0 1.852026 5.52  8.37157 17.62843
    Free Neu       11.2 1.852026 5.52  6.57157 15.82843
    Cued Pos       13.6 1.852026 5.52  8.97157 18.22843
    Free Pos       11.0 1.852026 5.52  6.37157 15.62843
   
   Confidence level used: 0.95