Apa yang harus dilakukan dengan korelasi efek acak yang sama dengan 1 atau -1?

Kejadian yang tidak biasa ketika berhadapan dengan model campuran maksimal yang kompleks (memperkirakan semua efek acak yang mungkin untuk data dan model yang diberikan) adalah sempurna (+1 atau -1) atau korelasi yang hampir sempurna di antara beberapa efek acak. Untuk tujuan diskusi, mari kita amati model dan ringkasan model berikut ini

Model: Y ~ X*Cond + (X*Cond|subj)

# Y = logit variable  
# X = continuous variable  
# Condition = values A and B, dummy coded; the design is repeated 
#             so all participants go through both Conditions  
# subject = random effects for different subjects  

Random effects:
 Groups  Name             Variance Std.Dev. Corr             
 subject (Intercept)      0.85052  0.9222                    
         X                0.08427  0.2903   -1.00            
         CondB            0.54367  0.7373   -0.37  0.37      
         X:CondB          0.14812  0.3849    0.26 -0.26 -0.56
Number of obs: 39401, groups:  subject, 219

Fixed effects:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)       2.49686    0.06909   36.14  < 2e-16 ***
X                -1.03854    0.03812  -27.24  < 2e-16 ***
CondB            -0.19707    0.06382   -3.09  0.00202 ** 
X:CondB           0.22809    0.05356    4.26 2.06e-05 ***

Alasan yang diduga di balik korelasi sempurna ini adalah bahwa kami telah menciptakan model yang terlalu rumit untuk data yang kami miliki. Saran umum yang diberikan dalam situasi ini adalah (misalnya, Matuschek et al., 2017; kertas ) untuk memperbaiki koefisien overparameter ke 0, karena model degenerasi seperti itu cenderung menurunkan daya. Jika kita mengamati perubahan nyata pada efek tetap dalam model yang diperkecil, kita harus menerimanya; jika tidak ada perubahan, maka tidak ada masalah dalam menerima yang asli.

Namun, mari kita asumsikan bahwa kita tidak hanya tertarik pada efek tetap yang dikendalikan untuk RE (efek acak), tetapi juga dalam struktur RE. Dalam kasus yang diberikan, secara teoritis akan masuk akal untuk mengasumsikan bahwa Interceptdan lereng Xmemiliki korelasi negatif bukan nol. Beberapa pertanyaan berikut:

1. Apa yang harus dilakukan dalam situasi seperti itu? Haruskah kita melaporkan korelasi sempurna dan mengatakan bahwa data kami tidak "cukup baik" untuk memperkirakan korelasi "nyata"? Atau haruskah kita melaporkan model korelasi 0? Atau haruskah kita mencoba menetapkan korelasi lain ke 0 dengan harapan bahwa yang "penting" tidak akan sempurna lagi? Saya rasa tidak ada 100% jawaban yang benar di sini, saya lebih suka mendengar pendapat Anda.

2. Bagaimana cara menulis kode yang memperbaiki korelasi 2 efek acak spesifik ke 0, tanpa mempengaruhi korelasi antara parameter lain?

Matriks kovarians efek-acak tunggal

Memperoleh perkiraan korelasi efek acak +1 atau -1 berarti algoritma optimasi mencapai "batas": korelasi tidak boleh lebih tinggi dari +1 atau lebih rendah dari -1. Bahkan jika tidak ada kesalahan atau peringatan konvergensi eksplisit, ini berpotensi mengindikasikan beberapa masalah dengan konvergensi karena kami tidak mengharapkan korelasi sebenarnya terletak pada batas. Seperti yang Anda katakan, ini biasanya berarti bahwa tidak ada cukup data untuk memperkirakan semua parameter dengan andal. Matuschek et al. 2017 mengatakan bahwa dalam situasi ini kekuatan dapat dikompromikan.

Cara lain untuk mencapai batas adalah untuk mendapatkan estimasi varians 0: Mengapa saya mendapatkan nol varians dari efek acak dalam model campuran saya, meskipun ada beberapa variasi dalam data?

$4\times 4$$4\times 4$(pracetak yang tidak diterbitkan) merekomendasikan menggunakan analisis komponen utama (PCA) untuk memeriksa apakah matriks kovarians yang diperoleh adalah singular. Jika ya, mereka menyarankan untuk memperlakukan situasi ini dengan cara yang sama seperti situasi tunggal di atas.

Jadi apa yang harus dilakukan?

Jika tidak ada cukup data untuk memperkirakan semua parameter model secara andal, maka kita harus mempertimbangkan untuk menyederhanakan model. Mengambil contoh model Anda,, X*Cond + (X*Cond|subj)ada berbagai cara yang memungkinkan untuk menyederhanakannya:

1. Hapus salah satu efek acak, biasanya korelasi tingkat tertinggi:

   X*Cond + (X+Cond|subj)

2. Singkirkan semua parameter korelasi:

   X*Cond + (X*Cond||subj)

   Pembaruan: sebagai catatan @Henrik, ||sintaks hanya akan menghapus korelasi jika semua variabel di sebelahnya adalah numerik. Jika variabel kategorikal (seperti Cond) terlibat, orang harusnya menggunakan afexpaket praktisnya (atau solusi manual yang rumit). Lihat jawabannya untuk lebih jelasnya.

3. Singkirkan beberapa parameter korelasi dengan memecah istilah menjadi beberapa, misalnya:

   X*Cond + (X+Cond|subj) + (0+X:Cond|subj)

4. Batasi matriks kovarians dengan cara tertentu, misalnya dengan menetapkan satu korelasi spesifik (yang mengenai batas) menjadi nol, seperti yang Anda sarankan. Tidak ada cara lme4bawaan untuk mencapai ini. Lihat jawaban @ BenBolker di SO untuk demonstrasi cara mencapai ini melalui beberapa peretasan cerdas.

Bertentangan dengan apa yang Anda katakan, saya tidak berpikir Matuschek et al. 2017 secara khusus merekomendasikan # 4. Inti dari Matuschek et al. 2017 dan Bates et al. 2015 tampaknya dimulai dengan model maksimal ala Barr et al. 2013 dan kemudian mengurangi kompleksitas sampai matriks kovarians peringkat penuh. (Selain itu, mereka akan sering merekomendasikan untuk mengurangi kompleksitas lebih jauh, untuk meningkatkan daya.) Pembaruan: Sebaliknya, Barr et al. merekomendasikan untuk mengurangi kompleksitas SAJA jika model tidak bertemu; mereka bersedia mentoleransi matriks kovarian tunggal. Lihat jawaban @ Henrik.

Jika seseorang setuju dengan Bates / Matuschek, maka saya pikir tidak masalah untuk mencoba berbagai cara mengurangi kompleksitas untuk menemukan orang yang melakukan pekerjaan sambil melakukan "kerusakan paling tidak". Melihat daftar saya di atas, matriks kovarians asli memiliki 10 parameter; # 1 memiliki 6 parameter, # 2 memiliki 4 parameter, # 3 memiliki 7 parameter. Model mana yang akan menghilangkan korelasi sempurna tidak mungkin dikatakan tanpa mencocokkannya.

Tetapi bagaimana jika Anda tertarik pada parameter ini?

Diskusi di atas memperlakukan matriks efek kovarians acak sebagai parameter gangguan. Anda mengajukan pertanyaan menarik tentang apa yang harus dilakukan jika Anda secara khusus tertarik pada parameter korelasi yang harus Anda "menyerah" untuk mendapatkan solusi peringkat penuh yang berarti.

Perhatikan bahwa memperbaiki parameter korelasi pada nol tidak harus menghasilkan BLUP ( ranef) yang tidak berkorelasi; bahkan, mereka mungkin tidak terlalu terpengaruh sama sekali (lihat jawaban @ Placidia untuk demonstrasi ). Jadi salah satu pilihan adalah melihat korelasi BLUP dan melaporkannya.

Pilihan lain, yang mungkin kurang menarik, adalah menggunakan suguhan subjectsebagai efek tetap Y~X*cond*subj, mendapatkan perkiraan untuk setiap mata pelajaran dan menghitung korelasi di antara mereka. Ini sama dengan menjalankan Y~X*condregresi terpisah untuk setiap subjek secara terpisah dan mendapatkan estimasi korelasi dari mereka.

Lihat juga bagian tentang model tunggal dalam model campuran Ben Bolker FAQ:

$\theta$

Saya setuju dengan semua yang dikatakan dalam jawaban amuba yang memberikan ringkasan bagus dari diskusi saat ini tentang masalah ini. Saya akan mencoba menambahkan beberapa poin tambahan dan merujuk pada handout kursus model campuran saya yang baru-baru ini yang juga merangkum poin-poin ini.

Menekan parameter korelasi (opsi 2 dan 3 dalam jawaban amuba) melalui ||karya hanya untuk kovariat numerik lmerdan bukan untuk faktor. Ini dibahas secara rinci dengan kode oleh Reinhold Kliegl .

Namun, afexpaket saya menyediakan fungsionalitas untuk menekan korelasi juga di antara faktor-faktor jika argumen expand_re = TRUEdalam panggilan untuk mixed()(lihat juga berfungsi lmer_alt()). Ini pada dasarnya melakukannya dengan menerapkan pendekatan yang dibahas oleh Reinhold Kliegl (yaitu, mentransformasikan faktor menjadi kovariat numerik dan menentukan struktur efek-acak pada faktor-faktor tersebut).

Contoh sederhana:

library("afex")
data("Machines", package = "MEMSS") # same data as in Kliegl code

# with correlation:
summary(lmer(score ~ Machine + (Machine  | Worker), data=Machines))
# Random effects:
#  Groups   Name        Variance Std.Dev. Corr       
#  Worker   (Intercept) 16.6405  4.0793              
#           MachineB    34.5467  5.8776    0.48      
#           MachineC    13.6150  3.6899   -0.37  0.30
#  Residual              0.9246  0.9616              
# Number of obs: 54, groups:  Worker, 6

## crazy results:
summary(lmer(score ~ Machine + (Machine  || Worker), data=Machines))
# Random effects:
#  Groups   Name        Variance Std.Dev. Corr     
#  Worker   (Intercept)  0.2576  0.5076            
#  Worker.1 MachineA    16.3829  4.0476            
#           MachineB    74.1381  8.6103   0.80     
#           MachineC    19.0099  4.3600   0.62 0.77
#  Residual              0.9246  0.9616            
# Number of obs: 54, groups:  Worker, 6

## as expected:
summary(lmer_alt(score ~ Machine + (Machine  || Worker), data=Machines))
# Random effects:
#  Groups   Name         Variance Std.Dev.
#  Worker   (Intercept)  16.600   4.0743  
#  Worker.1 re1.MachineB 34.684   5.8894  
#  Worker.2 re1.MachineC 13.301   3.6471  
#  Residual               0.926   0.9623  
# Number of obs: 54, groups:  Worker, 6

Bagi mereka yang tidak mengetahui afex, fungsi utama untuk model campuran adalah untuk memberikan nilai-p untuk efek tetap, misalnya ,:

(m1 <- mixed(score ~ Machine + (Machine  || Worker), data=Machines, expand_re = TRUE))
# Mixed Model Anova Table (Type 3 tests, KR-method)
# 
# Model: score ~ Machine + (Machine || Worker)
# Data: Machines
#    Effect      df        F p.value
# 1 Machine 2, 5.98 20.96 **    .002
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1

summary(m1)  
# [...]
# Random effects:
#  Groups   Name         Variance Std.Dev.
#  Worker   (Intercept)  27.4947  5.2435  
#  Worker.1 re1.Machine1  6.6794  2.5845  
#  Worker.2 re1.Machine2 13.8015  3.7150  
#  Residual               0.9265  0.9626  
# Number of obs: 54, groups:  Worker, 6
# [...]

Dale Barr dari Barr et al. (2013) makalah lebih berhati-hati dalam merekomendasikan pengurangan struktur efek-acak daripada disajikan dalam jawaban amuba. Dalam pertukaran twitter baru - baru ini dia menulis:

• "Mengurangi model memperkenalkan risiko antikonservativitas yang tidak diketahui, dan harus dilakukan dengan hati-hati, jika sama sekali." dan
• "Kekhawatiran utama saya adalah bahwa orang memahami risiko yang terkait dengan pengurangan model dan bahwa meminimalkan risiko ini memerlukan pendekatan yang lebih konservatif daripada yang biasanya diadopsi (misalnya setiap lereng diuji pada 0,05)."

Jadi hati-hati disarankan.

Sebagai salah satu pengulas saya juga dapat memberikan beberapa wawasan tentang mengapa kami the Bates et al. (2015) kertas tetap tidak dipublikasikan. Saya dan dua pengulas lainnya (yang menandatangani, tetapi akan tetap tidak disebutkan namanya di sini) memiliki beberapa kritik dengan pendekatan PCA (tampaknya tidak berprinsip dan tidak ada bukti bahwa itu lebih unggul dalam hal kekuasaan). Lebih lanjut, saya percaya ketiganya mengkritik bahwa makalah ini tidak fokus pada masalah bagaimana menentukan struktur efek-acak, tetapi juga mencoba untuk memperkenalkan GAMM. Dengan demikian, kertas Bates et al (2015) berubah menjadi Matuschek et al. (2017) makalah yang membahas masalah struktur efek acak dengan simulasi dan Baayen et al. (2017) makalah yang memperkenalkan GAMM.

Ulasan lengkap saya tentang Bates et al. draft dapat ditemukan di sini . IIRC, ulasan lain memiliki jenis poin utama yang serupa.

Saya juga memiliki masalah ini ketika menggunakan estimasi kemungkinan maksimum - hanya saya menggunakan algoritma Goldstein IGLS yang diimplementasikan melalui perangkat lunak MLwiN dan bukan LME4 di R. Namun dalam setiap kasus masalah telah teratasi ketika saya beralih ke estimasi MCMC menggunakan hal yang sama perangkat lunak. Saya bahkan memiliki korelasi lebih dari 3 yang diselesaikan ketika saya mengubah estimasi. Menggunakan IGLS, korelasi dihitung estimasi paska kovarians dibagi dengan produk dari akar kuadrat dari produk dari varians terkait - dan ini tidak memperhitungkan ketidakpastian dalam masing-masing estimasi konstituen.

Perangkat lunak IGLS tidak 'tahu' bahwa kovarian menyiratkan korelasi dan hanya menghitung perkiraan fungsi varians konstan, linier, kuadrat dll. Sebaliknya pendekatan MCMC dibangun berdasarkan asumsi sampel dari distribusi normal multivariat yang sesuai dengan varian dan kovarian dengan sifat baik dan propogasi kesalahan penuh sehingga ketidakpastian dalam estimasi kovarian diperhitungkan dalam estimasi varians dan sebaliknya.

MLwin menjadi rantai estimasi MCMC dengan estimasi IGLS dan matriks kovariansi varians non-negatif yang pasti mungkin perlu diubah dengan mengubah kovarians menjadi nol pada awalnya sebelum memulai pengambilan sampel.

Untuk contoh yang berhasil lihat

Mengembangkan model bertingkat untuk menganalisis kontekstual, heterogenitas, dan perubahan menggunakan MLwiN 3, Volume 1 (diperbarui September 2017); Volume 2 juga ada di RGate

https://www.researchgate.net/publication/320197425_Vol1Training_manualRevisedSept2017

Lampiran untuk Bab 10