Mengatasi heteroskedastisitas dalam Poisson GLMM

Saya memiliki data pengumpulan jangka panjang, dan saya ingin menguji, apakah jumlah hewan yang dikumpulkan dipengaruhi oleh efek cuaca. Model saya terlihat seperti di bawah ini:

glmer(SumOfCatch ~ I(pc.act.1^2) +I(pc.act.2^2) + I(pc.may.1^2) + I(pc.may.2^2) + 
                   SampSize + as.factor(samp.prog) + (1|year/month), 
      control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=1e9,npt=5)), 
      family="poisson", data=a2)

Penjelasan dari variabel yang digunakan:

• SumOfCatch: jumlah hewan yang dikumpulkan
• pc.act.1, pc.act.2: sumbu komponen utama yang mewakili kondisi cuaca selama pengambilan sampel
• pc.may.1, pc.may.2: sumbu PC yang mewakili kondisi cuaca di bulan Mei
• SampSize: jumlah perangkap lubang, atau mengumpulkan transek dari panjang standar
• samp.prog: metode pengambilan sampel
• tahun: tahun pengambilan sampel (dari 1993 hingga 2002)
• bulan: bulan pengambilan sampel (dari Agustus hingga November)

Residual model pas menunjukkan ketidakhomogenan yang cukup (heteroskedastisitas?) Ketika diplotkan dengan nilai pas (lihat Gbr.1):

Pertanyaan utama saya adalah: apakah ini masalah yang membuat keandalan model saya dipertanyakan? Jika demikian, apa yang bisa saya lakukan untuk menyelesaikannya?

Sejauh ini saya sudah mencoba yang berikut ini:

• kontrol untuk penyebaran berlebihan dengan mendefinisikan efek acak tingkat observasi, yaitu menggunakan ID unik untuk setiap pengamatan, dan menerapkan variabel ID ini sebagai efek acak; meskipun data saya memang menunjukkan overdispersi yang cukup besar, ini tidak membantu karena residu menjadi lebih buruk (lihat Gambar. 2)

• Saya memasang model tanpa efek acak, dengan quasi-Poisson glm dan glm.nb; juga menghasilkan plot residu vs dipasang serupa dengan model asli

Sejauh yang saya tahu, mungkin ada cara untuk estimasi kesalahan standar yang konsisten heteroskedastisitas, tapi saya gagal menemukan metode seperti itu untuk Poisson (atau jenis lain dari) GLMM di R.

Menanggapi @FlorianHartig: jumlah pengamatan dalam dataset saya adalah N = 554, saya pikir ini adalah hal yang wajar. nomor untuk model seperti itu, tetapi tentu saja, semakin banyak lebih meriah. Saya memposting dua angka, yang pertama adalah plot residual skala DHARMa (disarankan oleh Florian) dari model utama.

Angka kedua adalah dari model kedua, di mana satu-satunya perbedaan adalah bahwa ia berisi efek acak tingkat observasi (yang pertama tidak).

MEMPERBARUI

Gambar hubungan antara variabel cuaca (sebagai prediktor, yaitu sumbu x) dan keberhasilan pengambilan sampel (respons):

PEMBARUAN II.

Angka yang menunjukkan nilai prediktor vs residu:

Sulit untuk menilai kecocokan Poisson (atau GLM bernilai integer lainnya dalam hal ini) dengan Pearson atau residu penyimpangan, karena juga Poisson GLMM yang sangat cocok akan menunjukkan residu penyimpangan tidak homogen.

Ini terutama terjadi jika Anda melakukan GLMM dengan RE tingkat observasi, karena dispersi yang dibuat oleh OL-RE tidak dianggap oleh residu Pearson.

Untuk menunjukkan masalah ini, kode berikut membuat data Poisson overdispersed, yang kemudian dilengkapi dengan model yang sempurna. Residu Pearson sangat mirip dengan plot Anda - karenanya, mungkin tidak ada masalah sama sekali.

Masalah ini diselesaikan oleh paket DHARMa R , yang mensimulasikan dari model yang dipasang untuk mengubah residu dari setiap GL (M) M menjadi ruang standar. Setelah ini dilakukan, Anda dapat menilai / menguji masalah residual secara visual, seperti penyimpangan dari distribusi, ketergantungan residual pada prediktor, heteroskedastisitas, atau autokorelasi dengan cara normal. Lihat sketsa paket untuk contoh yang dikerjakan. Anda dapat melihat di plot bawah bahwa model yang sama sekarang terlihat baik-baik saja, sebagaimana mestinya.

Jika Anda masih melihat heteroskedastisitas setelah berkomplot dengan DHARMa, Anda harus memodelkan dispersi sebagai fungsi dari sesuatu, yang bukan masalah besar, tetapi kemungkinan akan mengharuskan Anda pindah ke JAG atau perangkat lunak Bayesia lainnya.

library(DHARMa)
library(lme4)

testData = createData(sampleSize = 200, overdispersion = 1, randomEffectVariance = 1, family = poisson())

fittedModel <- glmer(observedResponse ~ Environment1 + (1|group) + (1|ID), family = "poisson", data = testData, control=glmerControl(optCtrl=list(maxfun=20000) ))

# standard Pearson residuals
plot(fittedModel, resid(., type = "pearson") ~ fitted(.) , abline = 0)

# DHARMa residuals
plot(simulateResiduals(fittedModel))