Bagaimana menyesuaikan model campuran dengan variabel respons antara 0 dan 1?

15

Saya mencoba menggunakan lme4::glmer()agar sesuai dengan model campuran binomial umum (GLMM) dengan variabel dependen yang bukan biner, tetapi variabel kontinu antara nol dan satu. Orang dapat menganggap variabel ini sebagai probabilitas; sebenarnya itu adalah probabilitas sebagaimana dilaporkan oleh subyek manusia (dalam percobaan yang saya bantu analisis). Yaitu itu bukan fraksi "diskrit", tetapi variabel kontinu.

glmer()Panggilan saya tidak berfungsi seperti yang diharapkan (lihat di bawah). Mengapa? Apa yang dapat saya?

Sunting nanti: jawaban saya di bawah ini lebih umum daripada versi asli pertanyaan ini, jadi saya memodifikasi pertanyaannya menjadi lebih umum juga.


Keterangan lebih lanjut

Tampaknya dimungkinkan untuk menggunakan regresi logistik tidak hanya untuk DV biner tetapi juga untuk DV kontinu antara nol dan satu. Memang saat saya lari

glm(reportedProbability ~ a + b + c, myData, family="binomial")

Saya mendapat pesan peringatan

Warning message:
In eval(expr, envir, enclos) : non-integer #successes in a binomial glm!

tetapi kecocokan yang sangat masuk akal (semua faktor bersifat kategoris, jadi saya dapat dengan mudah memeriksa apakah prediksi model dekat dengan lintas-subjek-rata-rata, dan semuanya).

Namun, apa yang sebenarnya ingin saya gunakan adalah

glmer(reportedProbability ~ a + b + c + (1 | subject), myData, family="binomial")

Ini memberi saya peringatan identik, mengembalikan model, tetapi model ini jelas sangat tidak aktif; perkiraan efek tetap sangat jauh dari glm()yang dan dari lintas subjek berarti. (Dan saya perlu memasukkan glmerControl(optimizer="bobyqa")ke dalam glmerpanggilan, kalau tidak itu tidak menyatu sama sekali.)

amuba kata Reinstate Monica
sumber
1
Bagaimana dengan mentransformasikan probabilitas terlebih dahulu? Bisakah Anda mendapatkan sesuatu yang lebih dekat dengan distribusi normal dengan mengatakan, transformasi logit? Atau arcsin-sqrt? Itu akan menjadi preferensi saya daripada menggunakan glmer. Atau dalam solusi retas Anda, Anda juga dapat mencoba menambahkan efek acak untuk setiap pengamatan ke akun penyimpangan karena pilihan bobot Anda.
Aaron - Pasang kembali Monica
Terima kasih. Ya, saya bisa masuk ke DV dan kemudian menggunakan model campuran Gaussian (lmer), tetapi ini juga semacam peretasan, dan saya sudah membaca bahwa itu tidak disarankan. Saya akan mencoba efek acak untuk setiap pengamatan! Saat ini, saya mencoba model beta campuran; lme4 tidak bisa mengatasinya, tetapi glmmadmb bisa. Ketika saya menjalankan glmmadmb(reportedProbability ~ a + b + c + (1 | subject), myData, family="beta"), saya mendapatkan interval kepercayaan dan kewajaran yang benar, tetapi konvergensi gagal peringatan: - / Mencoba mencari cara untuk meningkatkan jumlah iterasi. Beta mungkin cocok untuk saya karena saya tidak memiliki DV = 0 atau DV = 1 case.
Amoeba mengatakan Reinstate Monica
Saya tidak tahu untuk glmer tetapi untuk glm ini dapat membantu: stats.stackexchange.com/questions/164120/… :
1
@ Harun: Saya mencoba menambahkan + (1 | rowid)panggilan telepon saya dan ini menghasilkan perkiraan stabil dan interval kepercayaan stabil, terlepas dari pilihan berat badan saya (saya mencoba 100 dan 500). Saya juga mencoba menjalankan lmer pada logit (melaporkanProbabilitas) dan saya mendapatkan hal yang hampir persis sama. Jadi kedua solusi tampaknya bekerja dengan baik! Beta MM dengan glmmadmb juga memberikan hasil yang sangat dekat, tetapi karena beberapa alasan gagal untuk menyatu sepenuhnya dan membutuhkan waktu lama untuk berjalan. Pertimbangkan memposting jawaban yang mencantumkan opsi ini dan menjelaskan sedikit perbedaan dan pro / kontra! (Interval kepercayaan yang saya sebutkan semuanya Wald.)
Amuba mengatakan Reinstate Monica
1
Dan mereka benar-benar yakin tentang nilai mereka seperti 0,9, atau apakah mereka juga memiliki "margin kesalahan"? Bisakah Anda berasumsi bahwa kepercayaan diri yang dilaporkan oleh subjek berbeda sama persis?

Jawaban:

20

Masuk akal untuk memulai dengan kasus sederhana tanpa efek acak.

Ada empat cara untuk berurusan dengan variabel respons nol-ke-satu terus-menerus yang berperilaku seperti pecahan atau probabilitas ( ini adalah utas kami yang paling kanonik / terunggul / dilihat pada topik ini, tetapi sayangnya tidak semua keempat opsi dibahas di sana):

  1. hal=m/nnnN

    glm(p ~ a+b+c, myData, family="binomial", weights=n)
  2. halhal01

    betareg(p ~ a+b+c, myData)
  3. Logit mengubah respons dan menggunakan regresi linier. Ini biasanya tidak disarankan.

    lm(log(p/(1-p)) ~ a+b+c, myData)
  4. Cocok dengan model binomial tapi kemudian hitung kesalahan standar dengan memperhitungkan over-dispersi. Kesalahan standar dapat dihitung dengan berbagai cara:

    • (a) meningkatkan kesalahan standar melalui estimasi penayangan berlebih ( satu , dua ). Ini disebut "quasi-binomial" GLM.

    • (B) kesalahan standar yang kuat melalui estimator sandwich ( satu , dua , tiga , empat ). Ini disebut "log fraksional" dalam ekonometrik.


    (A) dan (b) tidak identik (lihat komentar ini , dan bagian 3.4.1 dan 3.4.2 dalam buku ini , dan pos SO ini dan juga yang ini dan yang ini ), tetapi cenderung memberikan hasil yang serupa. Opsi (a) diimplementasikan glmsebagai berikut:

    glm(p ~ a+b+c, myData, family="quasibinomial")

Empat cara yang sama tersedia dengan efek acak.

  1. Menggunakan weightsargumen ( satu , dua ):

    glmer(p ~ a+b+c + (1|subject), myData, family="binomial", weights=n)

    Menurut tautan kedua di atas, mungkin ide yang baik untuk memodelkan penayangan berlebih, lihat di sana (dan juga # 4 di bawah).

  2. Menggunakan model beta campuran:

    glmmadmb(p ~ a+b+c + (1|subject), myData, family="beta")

    atau

    glmmTMB(p ~ a+b+c + (1|subject), myData, 
            family=list(family="beta",link="logit"))

    Jika ada nol atau satu dalam data respons, maka seseorang dapat menggunakan model beta nol / satu meningkat glmmTMB .

  3. Menggunakan transformasi logit dari respons:

    lmer(log(p/(1-p)) ~ a+b+c + (1|subject), myData)
  4. Akuntansi untuk penyebaran berlebihan dalam model binomial. Ini menggunakan trik yang berbeda: menambahkan efek acak untuk setiap titik data:

    myData$rowid = as.factor(1:nrow(myData))
    glmer(p ~ a+b+c + (1|subject) + (1|rowid), myData, family="binomial",
          glmerControl(optimizer="bobyqa"))

    Untuk beberapa alasan ini tidak berfungsi sebagaimana mestinya karena glmer()mengeluh tentang non-integer pdan menghasilkan perkiraan yang tidak masuk akal. Solusi yang saya buat adalah dengan menggunakan konstanta palsu weights=kdan memastikan bahwa p*kselalu bilangan bulat. Ini membutuhkan pembulatan ptetapi dengan memilih kyang cukup besar itu seharusnya tidak terlalu menjadi masalah. Hasilnya tampaknya tidak tergantung pada nilai k.

    k = 100
    glmer(round(p*k)/k ~ a+b+c + (1|subject) + (1|rowid), myData, 
          family="binomial", weights=rowid*0+k, glmerControl(optimizer="bobyqa"))

    Pembaruan kemudian (Jan 2018): Ini mungkin pendekatan yang tidak valid. Lihat diskusi di sini . Saya harus menyelidiki ini lebih lanjut.


Dalam kasus khusus saya, opsi # 1 tidak tersedia.

Opsi # 2 sangat lambat dan memiliki masalah dengan konvergensi: glmmadmbmembutuhkan lima-sepuluh menit untuk berjalan (dan masih mengeluh bahwa itu tidak konvergen!), Sedangkan lmerbekerja dalam sepersekian detik dan glmermembutuhkan beberapa detik. Pembaruan: Saya mencoba glmmTMBseperti yang disarankan dalam komentar oleh @BenBolker dan berfungsi hampir secepat glmer, tanpa masalah konvergensi. Jadi ini yang akan saya gunakan.

Opsi # 3 dan # 4 menghasilkan perkiraan yang sangat mirip dan interval kepercayaan Wald yang sangat mirip (diperoleh dengan confint). Saya bukan penggemar # 3 karena itu agak curang. Dan # 4 terasa agak berantakan.

Terima kasih banyak kepada @ Harun yang mengarahkan saya ke # 3 dan # 4 dalam komentarnya.

amuba kata Reinstate Monica
sumber
1
Jawaban yang bagus, dijelaskan dengan baik dan terhubung dengan model tanpa efek acak. Saya tidak akan menyebut kecurangan # 3 (transformasi), jenis-jenis transformasi itu sangat umum dalam analisis seperti ini. Sebagai gantinya saya akan mengatakan bahwa # 3 dan # 4 membuat asumsi tentang hubungan tentang distribusi data, dan begitu juga tentang hubungan antara mean dan varians, dan hanya karena # 4 memodelkan skala bahwa data dikumpulkan bukan berarti asumsi-asumsi itu akan menjadi lebih baik.
Aaron - Reinstate Monica
1
# 3 mengasumsikan logit dari probabilitas adalah normal dengan varians konstan, sedangkan # 4 mengasumsikan varians sebanding dengan p (1-p). Dari deskripsi Anda tentang kecocokan, ini tampaknya cukup mirip sehingga tidak terlalu menjadi masalah. Dan # 3 hampir pasti lebih standar (tergantung pada audiens Anda) jadi jika diagnostiknya masuk akal, itu yang saya inginkan.
Aaron - Reinstate Monica
1
kemungkinan lain adalah menggunakan glmmTMB ; setelah menginstal dengan devtools::install_github("glmmTMB/glmmTMB",sub="glmmTMB"), menggunakan glmmTMB(p ~ a+b+c + (1|subject), myData, family=list(family="beta",link="logit"))harus bekerja ...
Ben Bolker
@BenBolker Terima kasih! Apakah ada alasan untuk memilih glmmTMB daripada glmmADMB (untuk model beta) atau sebaliknya? Apakah salah satu dari paket ini lebih baru atau lebih aktif dikembangkan? Selain itu, bolehkah saya bertanya pendekatan apa di antara yang tercantum dalam jawaban ini - gaussian glmm setelah transformasi logit, beta glmm, atau binomial glmm dengan istilah (1 | rowid) - menurut Anda umumnya lebih disukai?
Amuba kata Reinstate Monica
1
Saya lebih suka beta GLMM jika memungkinkan - ini adalah model statistik yang dimaksudkan untuk mengukur perubahan proporsi lintas kovariat / kelompok. glmmTMBlebih cepat dan lebih stabil daripada glmmADMB, dan di bawah (sedikit) pengembangan lebih aktif, meskipun tidak matang.
Ben Bolker