Bagaimana cara mendapatkan nilai-p (periksa signifikansi) dari suatu efek dalam model campuran lme4?

56

Saya menggunakan lme4 dalam R agar sesuai dengan model campuran

lmer(value~status+(1|experiment)))

di mana nilai kontinu, status dan percobaan adalah faktor, dan saya dapatkan

Linear mixed model fit by REML 
Formula: value ~ status + (1 | experiment) 
  AIC   BIC logLik deviance REMLdev
 29.1 46.98 -9.548    5.911    19.1
Random effects:
 Groups     Name        Variance Std.Dev.
 experiment (Intercept) 0.065526 0.25598 
 Residual               0.053029 0.23028 
Number of obs: 264, groups: experiment, 10

Fixed effects:
            Estimate Std. Error t value
(Intercept)  2.78004    0.08448   32.91
statusD      0.20493    0.03389    6.05
statusR      0.88690    0.03583   24.76

Correlation of Fixed Effects:
        (Intr) statsD
statusD -0.204       
statusR -0.193  0.476

Bagaimana saya bisa tahu bahwa pengaruh status itu signifikan? R hanya melaporkan nilai- dan bukan nilai- .ptp

ECII
sumber
1
Pergi dengan jawaban yang diberikan untuk pertanyaan ini, orang bertanya-tanya dalam apa sebenarnya OP tertarik di sini: menguji koefisien terhadap nol (vanilla -tes satu tidak dalam regresi linier reguler terhadap nol ), atau pengujian untuk meminimalkan varians (Uji kami dapatkan dari banyak jenis ANOVA). Keduanya bertujuan pada hal-hal yang berbeda. Jawaban yang mencerahkan, meskipun bukan tentang model efek campuran, dapat ditemukan di sini . H 0 : β = β null FtH0:β=βnullF
Firebug

Jawaban:

61

Ada banyak informasi tentang topik ini di FAQ GLMM . Namun, dalam kasus khusus Anda, saya sarankan menggunakan

library(nlme)
m1 <- lme(value~status,random=~1|experiment,data=mydata)
anova(m1)

karena Anda tidak memerlukan barang apa pun yang lmermenawarkan (kecepatan lebih tinggi, penanganan efek acak silang, GLMMs ...). lmeharus memberikan persis sama koefisien dan varians estimasi tetapi juga akan df menghitung dan p-nilai untuk Anda (yang melakukan masuk akal dalam "klasik" desain seperti Anda tampaknya memiliki). Anda mungkin juga ingin mempertimbangkan istilah acak ~status|experiment(memungkinkan variasi efek status di seluruh blok, atau setara dengan interaksi status-per-eksperimen). Poster di atas juga benar bahwa tstatistik Anda sangat besar sehingga nilai-p Anda pasti akan <0,05, tetapi saya bisa membayangkan Anda ingin nilai-p "nyata".

Ben Bolker
sumber
3
Saya tidak tahu tentang jawaban ini. lmerbisa dengan mudah melaporkan jenis-nilai p yang sama tetapi tidak untuk alasan yang valid. Saya kira itu adalah komentar bahwa ada nilai-p "nyata" di sini yang mengganggu saya. Anda dapat berargumen bahwa Anda dapat menemukan satu kemungkinan cutoff, dan bahwa cutoff wajar apa pun dilewati. Tapi Anda tidak bisa berdebat ada nilai p nyata.
John
11
Untuk desain klasik (seimbang, bersarang, dll.) Saya pikir saya memang dapat berdebat bahwa ada p-vaue nyata, yaitu probabilitas untuk mendapatkan estimasi beta dari besaran yang diamati atau lebih besar jika hipotesis nol (beta = 0) itu salah ... lme4 tidak menyediakan penyebut df ini, saya percaya, karena lebih sulit untuk mendeteksi secara umum dari struktur model lme4 ketika model yang ditentukan adalah salah satu di mana beberapa heuristik untuk menghitung penyebut klasik df akan bekerja ...
Ben Bolker
mencoba summary(m1)bukan (saya menggunakan ini dengan paket nlme)
jena
36

Anda bisa menggunakan paket lmerTest . Anda cukup menginstal / memuatnya dan model lmer akan diperpanjang. Jadi misalnya

library(lmerTest)
lmm <- lmer(value~status+(1|experiment)))
summary(lmm)
anova(lmm)

akan memberi Anda hasil dengan nilai-p. Jika nilai-p adalah indikasi yang tepat agak sedikit diperdebatkan, tetapi jika Anda ingin memilikinya, ini adalah cara untuk mendapatkannya.

pbx101
sumber
28

Jika Anda dapat menangani meninggalkan nilai-p ( dan Anda harus ), Anda dapat menghitung rasio kemungkinan yang akan mewakili bobot bukti untuk efek status melalui:

#compute a model where the effect of status is estimated
unrestricted_fit = lmer(
    formula = value ~ (1|experiment) + status
    , REML = F #because we want to compare models on likelihood
)
#next, compute a model where the effect of status is not estimated
restricted_fit = lmer(
    formula = value ~ (1|experiment)
    , REML = F #because we want to compare models on likelihood
)
#compute the AIC-corrected log-base-2 likelihood ratio (a.k.a. "bits" of evidence)
(AIC(restricted_fit)-AIC(unrestricted_fit))*log2(exp(1))
Mike Lawrence
sumber
16
Perhatikan bahwa rasio kemungkinan asimptotik, yaitu jangan memperhitungkan ketidakpastian dalam estimasi varians residual ...
Ben Bolker
5
Saya tertarik pada baris terakhir Anda. Apa interpretasi hasilnya? Apakah ada sumber yang bisa saya lihat?
mguzmann
13

Masalahnya adalah bahwa perhitungan nilai-p untuk model-model ini tidak sepele, lihat diskusi di sini sehingga penulis lme4paket ini sengaja memilih untuk tidak memasukkan nilai-p dalam output. Anda mungkin menemukan metode penghitungan ini, tetapi mereka tidak selalu benar.

Michelle
sumber
9

Pertimbangkan apa yang Anda tanyakan. Jika Anda hanya ingin tahu apakah nilai p keseluruhan untuk efek status melewati semacam nilai cutoff sewenang-wenang, seperti 0,05, maka itu mudah. Pertama, Anda ingin mengetahui efek keseluruhannya. Anda bisa mendapatkannya dari anova.

m <- lmer(...) #just run your lmer command but save the model
anova(m)

Sekarang Anda memiliki nilai F. Anda bisa mengambilnya dan mencarinya di beberapa tabel F. Pilih saja denom serendah mungkin. derajat kebebasan. Batasnya akan ada sekitar 20. F Anda mungkin lebih besar dari itu tetapi saya bisa salah. Bahkan jika tidak, lihat jumlah derajat kebebasan dari perhitungan ANOVA konvensional di sini menggunakan jumlah percobaan yang Anda miliki. Menempel nilai itu di Anda turun menjadi sekitar 5 untuk cutoff. Sekarang Anda dengan mudah lulus dalam studi Anda. Df 'benar' untuk model Anda akan menjadi sesuatu yang lebih tinggi dari itu karena Anda memodelkan setiap titik data yang bertentangan dengan nilai agregat yang akan dimodelkan oleh ANOVA.

Jika Anda benar-benar menginginkan nilai-p yang tepat, tidak ada yang namanya kecuali Anda bersedia membuat pernyataan teoretis tentang hal itu. Jika Anda membaca Pinheiro & Bates (2001, dan mungkin beberapa buku lagi tentang masalah ini ... lihat tautan lain dalam jawaban ini) dan Anda keluar dengan argumen untuk df tertentu maka Anda dapat menggunakannya. Tapi Anda sebenarnya tidak mencari nilai p yang tepat. Saya menyebutkan ini karena Anda seharusnya tidak melaporkan nilai-p yang tepat, hanya saja cutoff Anda dilewati.

Anda harus benar-benar mempertimbangkan jawaban Mike Lawrence karena seluruh gagasan hanya bertahan dengan titik lulus untuk nilai-p sebagai informasi akhir dan paling penting untuk diekstrak dari data Anda pada umumnya salah arah (tetapi mungkin tidak ada dalam kasus Anda karena kami tidak Saya tidak punya cukup informasi untuk diketahui). Mike menggunakan versi hewan peliharaan dari perhitungan LR yang menarik, tetapi mungkin sulit untuk menemukan banyak dokumentasi di dalamnya. Jika Anda melihat pemilihan model dan interpretasi menggunakan AIC Anda mungkin menyukainya.

John
sumber
9

Sunting: Metode ini tidak lagi didukung di versi lme4 yang lebih baru. Gunakan paket lmerTest seperti yang disarankan dalam jawaban ini oleh pbx101 .

Ada posting di daftar R oleh penulis lme4 tentang mengapa nilai-p tidak ditampilkan. Dia menyarankan menggunakan sampel MCMC sebagai gantinya, yang Anda lakukan menggunakan pvals.fnc dari paket languageR:

library("lme4")
library("languageR")
model=lmer(...)
pvals.fnc(model)

Lihat http://www2.hawaii.edu/~kdrager/MixedEffectsModels.pdf untuk contoh dan detail.

Jeff
sumber
3
lme4 tidak lagi mendukung ini. Posting ini dapat diperbarui untuk orang-orang cadangan yang harus mengetahui hal ini seperti yang baru saja saya lakukan.
timothy.s.lau
5

Apakah Anda tertarik mengetahui apakah efek gabungan dari statusberpengaruh signifikan value? Jika demikian, Anda dapat menggunakan Anovafungsi dalam carpaket (jangan bingung dengan anovafungsi di pangkalan R).

dat <- data.frame(
  experiment = sample(c("A","B","C","D"), 264, replace=TRUE), 
  status = sample(c("D","R","A"), 264, replace=TRUE), 
  value = runif(264)   
)
require(lme4)
(fm <- lmer(value~status+(1|experiment), data=dat))

require(car)
Anova(fm)

Lihatlah ?Anovasetelah memuat carpaket.

smillig
sumber
Adakah yang tahu bagaimana cara car::Anova()menghindari masalah lengket seputar perhitungan nilai-p yang ditautkan Michelle?
Mike Lawrence
Aku tidak, tapi tebakanku adalah menghindari masalah lengket dengan mengabaikannya! Setelah membaca ulang posting asli, saya merasa bahwa saya mungkin telah salah mengerti pertanyaan itu. Jika OP menginginkan nilai p yang tepat untuk parameter efek tetap, ia sedang dalam masalah. Tetapi jika OP hanya ingin tahu apakah itu signifikan, saya pikir nilai-t lebih besar daripada ketidakpastian dalam bagaimana nilai-p tepatnya akan dihitung. (Dengan kata lain, itu penting.)
smillig
1
Saya pikir itu pasti ide yang bagus untuk mengarahkan kembali ke perhitungan ANOVA untuk mengetahui efek keseluruhan dari statistik tetapi saya tidak yakin bahwa nilai-p itu baik. anovaPerintah reguler akan memberi Anda nilai F.
John
Saya pikir ini sedikit lebih lengket daripada yang terlihat. Menjalankan ANOVA's valid ketika Anda ingin meminimalkan varians, tetapi dari pertanyaan kata saya pikir OP ingin menetapkan efek marginal dari variabel, yaitu uji koefisien terhadap nol.
Firebug
0

Fungsi pvals.fncini tidak lagi didukung oleh lme4. Menggunakan paket paket lmerTest, dimungkinkan untuk menggunakan metode lain untuk menghitung nilai-p, seperti perkiraan Kenward-Roger

model=lmer(value~status+1|experiment)
anova(model, ddf="Kenward-Roger")
Stefano Cacciatore
sumber
0

Hanya dengan memuat paket afex akan mencetak nilai-p dalam output fungsi lmer dari paket lme4 (Anda tidak perlu menggunakan afex; muat saja):

library(lme4)  #for mixed model
library(afex)  #for p-values

Ini akan secara otomatis menambahkan kolom nilai p ke output lmer (model Anda) untuk efek tetap.

Maryam Nasseri
sumber