Plot diagnostik untuk regresi jumlah

88

Plot diagnostik apa (dan mungkin tes formal) yang menurut Anda paling informatif untuk regresi di mana hasilnya adalah variabel hitungan?

Saya terutama tertarik pada Poisson dan model binomial negatif, serta rekan-rekan nol-inflated dan rintangan masing-masing. Sebagian besar sumber yang saya temukan hanya memplot nilai residual vs nilai tanpa diskusi tentang seperti apa plot ini "seharusnya".

Kebijaksanaan dan referensi sangat dihargai. Kisah belakang tentang mengapa saya menanyakan hal ini, jika relevan, adalah pertanyaan saya yang lain .

Diskusi terkait:

setengah lulus
sumber

Jawaban:

101

Inilah yang biasanya saya sukai (untuk ilustrasi saya menggunakan data quine yang terlalu banyak disebarkan dan tidak mudah dimodelkan dari hari-hari murid absen dari sekolah MASS):

  1. Tes dan grafik data penghitungan asli dengan memplot frekuensi yang diamati dan frekuensi yang sesuai (lihat bab 2 dalam Friendly ) yang didukung oleh vcdpaket dalam Rsebagian besar. Misalnya, dengan goodfitdan rootogram:

    library(MASS)
    library(vcd)
    data(quine) 
    fit <- goodfit(quine$Days) 
    summary(fit) 
    rootogram(fit)

    atau dengan plot Ord yang membantu mengidentifikasi model data hitungan mana yang mendasarinya (misalnya, di sini kemiringan positif dan intersepnya positif yang berbicara untuk distribusi binomial negatif):

    Ord_plot(quine$Days)

    atau dengan plot "XXXXXXness" di mana XXXXX adalah distribusi pilihan, katakan plot Poissoness (yang berbicara menentang Poisson, coba juga type="nbinom"):

    distplot(quine$Days, type="poisson")
  2. Periksa pengukuran good-of-fit yang biasa (seperti statistik rasio kemungkinan vs model nol atau serupa):

    mod1 <- glm(Days~Age+Sex, data=quine, family="poisson")
    summary(mod1)
    anova(mod1, test="Chisq")
  3. Periksa over / underdispersion dengan melihat residual deviance/dfatau pada statistik uji formal (misalnya, lihat jawaban ini ). Di sini kami memiliki penayangan berlebih:

    library(AER)
    deviance(mod1)/mod1$df.residual
    dispersiontest(mod1)
  4. Periksa poin pengaruh dan pengaruh , misalnya, dengan yang ada influencePlotdalam carpaket. Tentu saja di sini banyak poin yang sangat berpengaruh karena Poisson adalah model yang buruk:

    library(car)
    influencePlot(mod1)
  5. Periksa nol inflasi dengan memasang model data hitungan dan mitra zeroinflated / hurdle dan membandingkannya (biasanya dengan AIC). Di sini model nol yang digelembungkan akan lebih cocok daripada Poisson sederhana (sekali lagi mungkin karena penyebaran berlebihan):

    library(pscl)
    mod2 <- zeroinfl(Days~Age+Sex, data=quine, dist="poisson")
    AIC(mod1, mod2)
  6. Plot residual (mentah, penyimpangan atau diskalakan) pada sumbu y vs nilai prediksi (log) (atau prediktor linier) pada sumbu x. Di sini kita melihat beberapa residu yang sangat besar dan penyimpangan substansial residu penyimpangan dari normal (berbicara melawan Poisson; Edit: @ FlorianHartig menjawab bahwa normalitas residu ini tidak diharapkan sehingga ini bukan petunjuk konklusif):

    res <- residuals(mod1, type="deviance")
    plot(log(predict(mod1)), res)
    abline(h=0, lty=2)
    qqnorm(res)
    qqline(res)
  7. Jika berminat, plot setengah plot probabilitas normal residual dengan memplot residu absolut terurut vs. nilai normal yang diharapkan Atkinson (1981) . Fitur khusus akan mensimulasikan 'garis' referensi dan amplop dengan interval kepercayaan yang disimulasikan / dibooting (tidak ditampilkan):

    library(faraway)
    halfnorm(residuals(mod1))
  8. ±

    plot(Days~Age, data=quine) 
    prs  <- predict(mod1, type="response", se.fit=TRUE)
    pris <- data.frame("pest"=prs[[1]], "lwr"=prs[[1]]-prs[[2]], "upr"=prs[[1]]+prs[[2]])
    points(pris$pest ~ quine$Age, col="red")
    points(pris$lwr  ~ quine$Age, col="pink", pch=19)
    points(pris$upr  ~ quine$Age, col="pink", pch=19)

Ini akan memberi Anda banyak informasi berguna tentang analisis Anda dan sebagian besar langkah berfungsi untuk semua distribusi data jumlah standar (misalnya, Poisson, Binomial Negatif, COM Poisson, Power Laws).

Momo
sumber
6
Jawaban menyeluruh yang bagus! Sangat membantu untuk juga menjalankan diagnostik ini dengan data yang disimulasikan Poisson untuk melatih mata saya dengan seperti apa plotnya.
setengah jalan
Haruskah saya memberikan penjelasan lebih lanjut tentang apa yang dilakukan plot atau apakah tidak apa-apa seperti ini?
Momo
2
Catatan samping yang menarik: Saya menemukan bahwa distribusi NB tampaknya tidak sesuai dengan data NB yang disimulasikan berdasarkan uji GOF, rootogram, plot Ord, dan plot NB-ness. Pengecualian tampaknya sangat "jinak" data NB yang hampir simetris - mu tinggi, theta tinggi.
setengah jalan
1
Hm, apakah Anda yakin Anda menggunakan type = "nbinomial" sebagai argumen? Misalnya fm <- glm.nb (Hari ~., Data = quine); dummy <- rnegbin (fitting (fm), theta = 4.5) berfungsi dengan baik.
Momo
@Momo, terima kasih - Saya melakukan sesuatu seperti x = rnegbin (n = 1000, mu = 10, theta = 1); fit = goodfit (x, type = "nbinomial"); ringkasan (fit). Pengaturan theta = 4,5 memang meningkatkan kecocokan, tetapi masih sering p <0,05 dan rootogram bisa terlihat sangat buruk. Hanya supaya saya mengerti perbedaan antara simulasi kami: pada Anda, setiap nilai dummy disimulasikan dari parameter rata-rata yang berbeda (nilai terpasang (fm)), kan? Di saya, mereka semua memiliki rata-rata 10.
setengah lulus
14

Untuk pendekatan menggunakan plot diagnostik standar tetapi ingin tahu seperti apa bentuknya, saya suka makalahnya:

 Buja, A., Cook, D. Hofmann, H., Lawrence, M. Lee, E.-K., Swayne,
 D.F and Wickham, H. (2009) Statistical Inference for exploratory
 data analysis and model diagnostics Phil. Trans. R. Soc. A 2009
 367, 4361-4383 doi: 10.1098/rsta.2009.0120

Salah satu pendekatan yang disebutkan di sana adalah untuk membuat beberapa set data simulasi di mana asumsi-asumsi yang menarik adalah benar dan membuat plot diagnostik untuk set data-data yang disimulasikan ini dan juga membuat plot diagnostik untuk data nyata. letakkan semua plot ini di layar pada saat yang sama (secara acak menempatkan plot berdasarkan data nyata). Sekarang Anda memiliki referensi visual tentang seperti apa plot seharusnya dan jika asumsi berlaku untuk data nyata maka plot tersebut akan terlihat seperti yang lain (jika Anda tidak tahu mana data sebenarnya, maka asumsi yang diuji kemungkinan dekat) cukup untuk benar), tetapi jika plot data nyata terlihat jelas berbeda dari yang lain, maka itu berarti bahwa setidaknya satu asumsi tidak berlaku. The vis.testfungsi dalam paket TeachingDemos untuk R membantu menerapkan ini sebagai ujian.

Greg Snow
sumber
6
Contoh dengan data di atas, untuk catatan: mod1 <- glm (Hari ~ Usia + Jenis Kelamin, data = quine, family = "poisson"); if (interactive ()) {vis.test (residual (mod1, type = "response"), vt.qqnorm, nrow = 5, ncol = 5, npage = 3)}
setengah jalan
13

Ini adalah pertanyaan lama, tetapi saya pikir akan bermanfaat untuk menambahkan bahwa paket DHARMa R saya (tersedia dari CRAN, lihat di sini ) sekarang menyediakan residu standar untuk GLM dan GLMM, berdasarkan pada pendekatan simulasi yang serupa dengan apa yang disarankan oleh @GregSnow .

Dari deskripsi paket:

Paket DHARMa menggunakan pendekatan berbasis simulasi untuk membuat sisa skala yang dapat ditafsirkan dengan mudah dari model campuran linier umum yang sesuai. Yang didukung saat ini adalah semua kelas 'merMod' dari 'lme4' ('lmerMod', 'glmerMod'), 'glm' (termasuk 'negbin' dari 'MASS', tetapi tidak termasuk distribusi kuasi) dan kelas model 'lm'. Atau, simulasi yang dibuat secara eksternal, misalnya simulasi prediksi posterior dari perangkat lunak Bayesian seperti 'JAGS', 'STAN', atau 'BUGS' dapat diproses juga. Residu yang dihasilkan distandarisasi ke nilai antara 0 dan 1 dan dapat diinterpretasikan secara intuitif sebagai residu dari regresi linier. Paket ini juga menyediakan sejumlah fungsi plot dan uji untuk masalah kesalahan spesifikasi model,

@Momo - Anda mungkin ingin memperbarui rekomendasi Anda 6, itu menyesatkan. Normalitas residual penyimpangan pada umumnya tidak diharapkan di bawah Poisson , seperti yang dijelaskan dalam sketsa DHARMa atau di sini ; dan melihat residu penyimpangan (atau residual standar lainnya) yang berbeda dari garis lurus dalam plot qqnorm karena itu pada umumnya tidak ada perhatian sama sekali . Paket DHARMa menyediakan plot qq yang dapat diandalkan untuk mendiagnosis penyimpangan dari Poisson atau keluarga GLM lainnya. Saya telah membuat contoh yang menunjukkan masalah dengan residu penyimpangan di sini .

Florian Hartig
sumber
4

Ada fungsi yang disebut glm.diag.plotsdalam paket boot, untuk menghasilkan plot diagnostik untuk GLM. Apa fungsinya:

Membuat plot residu penyimpangan jackknife terhadap prediktor linier, plot skor normal dari residual deviansi standar, plot perkiraan statistik Cook terhadap leverage / (1-leverage), dan plot kasus statistik Cook.

kdarras
sumber