Nilai P untuk istilah interaksi dalam model efek campuran menggunakan lme4

Saya menganalisa beberapa data perilaku menggunakan lme4di R, sebagian besar mengikuti tutorial Bodo Winter baik , tapi saya tidak mengerti jika saya menangani interaksi dengan benar. Lebih buruk lagi, tidak ada orang lain yang terlibat dalam penelitian ini menggunakan model campuran, jadi saya sedikit terombang-ambing ketika datang untuk memastikan semuanya benar.

Alih-alih hanya mengirim teriakan minta tolong, saya pikir saya harus melakukan upaya terbaik saya untuk menafsirkan masalah, dan kemudian meminta koreksi kolektif Anda. Beberapa selainnya adalah:

• Saat menulis, saya menemukan pertanyaan ini , menunjukkan bahwa nlmelebih langsung memberikan nilai p untuk istilah interaksi, tapi saya pikir itu masih valid untuk ditanyakan sehubungan dengan lme4.
• Livius'jawaban untuk pertanyaan ini memberikan tautan ke banyak bacaan tambahan, yang akan saya coba selesaikan dalam beberapa hari ke depan, jadi saya akan berkomentar dengan kemajuan apa pun yang membawa.

Dalam data saya, saya memiliki variabel dependen dv, conditionmanipulasi (0 = kontrol, 1 = kondisi eksperimental, yang akan menghasilkan yang lebih tinggi dv), dan juga prasyarat, berlabel appropriate: uji coba berkode 1untuk ini harus menunjukkan efeknya, tetapi uji coba berkode 0mungkin tidak, karena faktor penting hilang.

Saya juga menyertakan dua intersep acak, untuk subject, dan untuk target, yang mencerminkan dvnilai-nilai yang berkorelasi dalam setiap subjek, dan dalam masing-masing dari 14 masalah yang diselesaikan (masing-masing peserta menyelesaikan kontrol dan versi eksperimental dari setiap masalah).

library(lme4)
data = read.csv("data.csv")

null_model        = lmer(dv ~ (1 | subject) + (1 | target), data = data)
mainfx_model      = lmer(dv ~ condition + appropriate + (1 | subject) + (1 | target),
                         data = data)
interaction_model = lmer(dv ~ condition + appropriate + condition*appropriate +
                              (1 | subject) + (1 | target), data = data)
summary(interaction_model)

Keluaran:

## Linear mixed model fit by REML ['lmerMod']
## ...excluded for brevity....
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept) 0.006594 0.0812  
##  target   (Intercept) 0.000557 0.0236  
##  Residual             0.210172 0.4584  
## Number of obs: 690, groups: subject, 38; target, 14
## 
## Fixed effects:
##                                Estimate Std. Error t value
## (Intercept)                    0.2518     0.0501    5.03
## conditioncontrol               0.0579     0.0588    0.98
## appropriate                   -0.0358     0.0595   -0.60
## conditioncontrol:appropriate  -0.1553     0.0740   -2.10
## 
## Correlation of Fixed Effects:
## ...excluded for brevity.

ANOVA kemudian menunjukkan interaction_modelsecara signifikan lebih cocok daripada mainfx_model, dari mana saya menyimpulkan bahwa ada interaksi hadir yang signifikan (p = 0,035).

anova(mainfx_model, interaction_model)

Keluaran:

## ...excluded for brevity....
##                   Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)  
## mainfx_model       6 913 940   -450      901                          
## interaction_model  7 910 942   -448      896  4.44      1      0.035 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dari sana, saya mengisolasi bagian dari data yang memenuhi appropriatepersyaratan (yaitu, appropriate = 1), dan untuk itu cocok dengan model nol, dan model termasuk conditionsebagai efek, membandingkan dua model menggunakan ANOVA lagi, dan lihat, menemukan bahwa conditionadalah prediktor yang signifikan.

good_data = data[data$appropriate == 1, ]
good_null_model   = lmer(dv ~ (1 | subject) + (1 | target), data = good_data)
good_mainfx_model = lmer(dv ~ condition + (1 | subject) + (1 | target), data = good_data)

anova(good_null_model, good_mainfx_model)

Keluaran:

## Data: good_data
## models:
## good_null_model: dv ~ (1 | subject) + (1 | target)
## good_mainfx_model: dv ~ condition + (1 | subject) + (1 | target)
##                   Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)  
## good_null_model    4 491 507   -241      483                          
## good_mainfx_model  5 487 507   -238      477  5.55      1      0.018 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Saya tidak melihat terlalu banyak untuk dikatakan di sini. Saya pikir Anda telah melakukan pekerjaan dengan baik.

Ada beberapa cara yang telah dibahas orang untuk menguji efek dan mendapatkan nilai-p untuk model efek campuran yang rumit. Ada ikhtisar yang bagus di sini . Yang terbaik adalah menggunakan metode intensif komputasi (bootstrap atau metode Bayesian), tetapi ini lebih maju untuk kebanyakan orang. Metode terbaik kedua (dan paling nyaman) adalah menggunakan uji rasio kemungkinan. Itulah yang anova()(secara teknis ? Anova.merMod () ) lakukan. Penting untuk hanya menggunakan uji rasio kemungkinan model yang sesuai dengan kemungkinan maksimum penuh , daripada kemungkinan maksimum terbatas.(REML). Di sisi lain, untuk model akhir Anda, dan untuk interpretasi, Anda ingin menggunakan REML. Ini membingungkan banyak orang. Dalam output Anda, kami melihat bahwa Anda cocok dengan model Anda dengan REML (ini karena opsi diatur TRUEsecara default di lmer(). Itu berarti bahwa tes Anda tidak valid, namun, karena ini adalah kesalahan umum, anova.merMod()berisi refitargumen yang oleh default diatur ke TRUE, dan Anda tidak mengubahnya. Jadi pandangan ke depan dari pengembang paket telah menyelamatkan Anda di sana.

Mengenai strategi Anda untuk membongkar interaksi, apa yang Anda lakukan baik-baik saja. Ingatlah bahwa interaksi menggunakan semua data untuk pengujiannya. Dimungkinkan untuk memiliki interaksi yang signifikan tetapi tidak memiliki tes bertingkat yang signifikan, yang membingungkan beberapa orang. (Tapi sepertinya itu tidak terjadi padamu.)

Saya sendiri seorang pemula, dan mengikuti saran dari Zuur dkk. Saya menggunakan lmedari nlmepaket alih-alih lme4ketika saya perlu menambahkan struktur kesalahan hierarkis ke model linier yang sebaliknya. Respons saya mungkin jauh.

Dua komentar:

(1) Saya tidak yakin masuk akal untuk menguji conditiondi subset kapan appropriate==1saja. Jika Anda ingin mendapatkan nilai p untuk efek utama, Anda dapat menggunakan Anovadari paket 'mobil':

require(car)
Anova(M,type="III")# type III Sum of Squares. M was fitted using method=REML

Jika Anda ingin menyelesaikan interaksi, Anda dapat menjalankan perbandingan berpasangan secara langsung (?) Atau melakukan apa yang Anda lakukan tetapi pada kedua himpunan bagian (yaitu juga dengan subset di mana appropriate==0).

(2) Anda mungkin ingin memilih struktur kesalahan Anda terlebih dahulu alih-alih berasumsi (1 | subject) + (1 | target)adalah struktur kesalahan terbaik. Dari apa yang Anda tulis, saya kumpulkan bahwa itu conditionadalah dalam faktor subjek, sementara itu appropriateadalah antara subjek atau antara target. Anda mungkin ingin menambahkan lereng untuk faktor dalam subjek dan / atau dalam target, misalnya: dv ~ condition + appropriate + (1+condition | subject) + (1 | target)menambahkan kemiringan acak untuk faktor dalam subjek condition. Tidak ada kemiringan yang diperlukan untuk faktor antar subjek / target.

Bersulang