Bagaimana saya bisa mempercepat perhitungan efek tetap dalam GLMM?

9

Saya sedang melakukan studi simulasi yang memerlukan perkiraan bootstrap yang diperoleh dari model campuran linier umum (sebenarnya, produk dari dua perkiraan untuk efek tetap, satu dari GLMM dan satu dari LMM). Untuk melakukan penelitian dengan baik akan membutuhkan sekitar 1000 simulasi dengan 1000 atau 1500 replikasi bootstrap setiap kali. Ini membutuhkan banyak waktu di komputer saya (berhari-hari).

How can I speed up the computation of these fixed effects?

Untuk lebih spesifik, saya memiliki subjek yang diukur berulang kali dalam tiga cara, sehingga memunculkan variabel X, M, dan Y, di mana X dan M adalah kontinu dan Y adalah biner. Kami memiliki dua persamaan regresi mana Y adalah variabel laten yang mendasari berkelanjutan untuk dan kesalahan tidak iid. Statistik yang ingin kita bootstrap adalah . Jadi, setiap replikasi bootstrap membutuhkan pemasangan LMM dan GLMM. Kode R saya adalah (menggunakan lme4)

M=α0+α1X+ϵ1
Y α 1 β 2
Y=β0+β1X+β2M+ϵ2
Y
α1β2
    stat=function(dat){
        a=fixef(lmer(M~X+(X|person),data=dat))["X"]
        b=fixef(glmer(Y~X+M+(X+M|person),data=dat,family="binomial"))["M"]
        return(a*b)
    }

Saya menyadari bahwa saya mendapatkan taksiran yang sama untuk jika saya pas sebagai model linier, sehingga menghemat waktu, tetapi trik yang sama tidak berhasil untuk .β 2α1β2

Apakah saya hanya perlu membeli komputer yang lebih cepat? :)

BR
sumber
1
@BR apa leher botolnya di sini? Pada dasarnya apa yang menyita waktu Rprof.
suncoolsu
1
Salah satu caranya adalah dengan mengabaikan bagian "campuran" dari GLMM. Sesuai dengan GLM biasa, perkiraan tidak akan banyak berubah, tetapi kesalahan standar mereka mungkin akan
probabilityislogic
@probabilityislogic. Selain komentar Anda, saya juga berpikir, apakah jawaban akan berbeda banyak tergantung pada ukuran kelompok dan perilaku individu dalam kelompok. Seperti yang dikatakan Gelman dan Hill: hasil model efek campuran akan antara pooling dan no pooling. (Obv. Ini untuk model hierarki Bayesian, tetapi model campuran adalah cara klasik untuk melakukan hal yang sama.)
suncoolsu
@probabilityislogic: Itu bekerja untuk LMM, tetapi tampaknya gagal untuk GLMM (artinya saya menjalankan model dengan dan tanpa M ekstra pada data yang sama dan berakhir dengan hasil yang sangat berbeda). Kecuali, tentu saja, ada kesalahan dalam implementasi glmer.
BR
@suncoolsu: bottleneck adalah estimasi GLMM, yang mungkin beberapa detik (terutama dengan beberapa efek acak). Tapi lakukan itu 1000 * 1000 kali, dan itu 280 jam perhitungan. Pemasangan GLM membutuhkan waktu sekitar 1/100 waktu.
BR

Jawaban:

4

Seharusnya membantu menentukan nilai awal, meskipun sulit untuk mengetahui berapa banyak. Saat Anda melakukan simulasi dan bootstrap, Anda harus mengetahui nilai 'true' atau perkiraan un-bootstrap atau keduanya. Coba gunakan yang ada di start =opsi glmer.

Anda juga dapat mempertimbangkan untuk melihat apakah toleransi untuk menyatakan konvergensi lebih ketat dari yang Anda butuhkan. Saya tidak jelas bagaimana cara mengubahnya dari lme4dokumentasi.

onestop
sumber
4

Dua kemungkinan lain juga dipertimbangkan, sebelum membeli komputer baru.

  1. Komputasi paralel - bootstrap mudah dijalankan secara paralel. Jika komputer Anda cukup baru, Anda mungkin memiliki empat inti. Lihatlah perpustakaan multicore di R.
  2. Komputasi awan juga kemungkinan dan cukup murah. Saya memiliki kolega yang telah menggunakan cloud amazon untuk menjalankan skrip R. Mereka mendapati bahwa itu cukup hemat biaya.
csgillespie
sumber
1
Terima kasih atas jawabannya. Entah bagaimana, saya telah mengabaikan fakta bahwa saya memiliki dua core (komputer saya tidak terlalu baru). Saya seharusnya sudah melihat multicore sejak lama.
BR
2

Ini mungkin komputer yang lebih cepat. Tapi di sini ada satu trik yang mungkin berhasil.

Menghasilkan simultation dari , tetapi hanya tergantung pada Y , kemudian hanya melakukan OLS atau LMM pada simulasi Y * nilai.YYY

Andaikan fungsi tautan Anda adalah . ini mengatakan bagaimana Anda dapatkan dari probabilitas Y = 1 ke Y * nilai, dan kemungkinan besar logistik fungsi g ( z ) = l o g ( zg(.)Y=1Y.g(z)=log(z1z)

Jadi jika Anda mengasumsikan distribusi sampel bernouli untuk , dan kemudian menggunakan jeffrey sebelum probabilitas, Anda mendapatkan beta posterior untuk p B e t a ( Y o b s + 1YYBernoulli(p). Simulasi dari ini harus seperti pencahayaan, dan jika tidak, maka Anda memerlukan komputer yang lebih cepat. Selanjutnya, sampel independen, jadi tidak perlu memeriksa diagnostik "konvergensi" seperti di MCMC, dan Anda mungkin tidak perlu sampel sebanyak - 100 dapat bekerja dengan baik untuk kasus Anda. Jika Anda memiliki binomialYs, maka ganti saja angka1di posterior di atas denganni, jumlah percobaan binomial untuk setiapYi.pBeta(Yobs+12,1Yobs+12)Ys1niYi

psimYsim=g(psim)YssayamN×SNS

gmler()lmer()Ylmer()b

Sebuah=...
b=0
dHai s=1,...,S
best=lmer(Ys...)
b=b+1s(best-b)
end
retkamurn(Sebuahb)

Beri tahu saya jika saya perlu menjelaskan sesuatu yang sedikit lebih jelas

probabilityislogic
sumber
Terima kasih atas jawabannya, saya butuh sedikit waktu untuk mencernanya (dan saya sudah punya rencana untuk Sabtu malam). Cukup berbeda sehingga tidak jelas bagi saya jika memberikan jawaban yang sama dengan pendekatan GLMM, tetapi saya harus lebih memikirkannya.
BR