Mengapa kuasi-Poisson di GLM tidak diperlakukan sebagai kasus khusus binomial negatif?

21

Saya mencoba menyesuaikan model linier umum untuk beberapa set data jumlah yang mungkin atau mungkin tidak disebarkan secara berlebihan. Dua distribusi kanonik yang berlaku di sini adalah Poisson dan Negative Binomial (Negbin), dengan EV dan variansμ

VSebuahrP=μ

VSebuahrNB=μ+μ2θ

yang dapat dipasang di R menggunakan glm(..,family=poisson)dan glm.nb(...), masing-masing. Ada juga quasipoissonkeluarga, yang dalam pemahaman saya adalah Poisson yang disesuaikan dengan EV dan varians yang sama

VSebuahrQP=ϕμ ,

yaitu jatuh di suatu tempat di antara Poisson dan Negbin. Masalah utama dengan keluarga quasipoisson adalah bahwa tidak ada kemungkinan yang sesuai untuk itu, dan karenanya banyak tes statistik yang sangat berguna dan langkah-langkah kecocokan (AIC, LR dan sebagainya) tidak tersedia.

Jika Anda membandingkan varian QP dan Negbin, Anda mungkin memperhatikan bahwa Anda dapat menyamakannya dengan meletakkan . Melanjutkan logika ini, Anda dapat mencoba untuk mengekspresikan distribusi quasipoisson sebagai kasus khusus dari Negbin:ϕ=1+μθ

QP(μ,ϕ)=NB(μ,θ=μϕ-1) ,

yaitu Negbin dengan secara linear tergantung pada . Saya mencoba memverifikasi ide ini dengan membuat urutan angka acak sesuai dengan rumus di atas dan menyesuaikannya dengan :θμglm

#fix parameters

phi = 3
a = 1/50
b = 3
x = 1:100

#generating points according to an exp-linear curve
#this way the default log-link recovers the same parameters for comparison

mu = exp(a*x+b) 
y = rnbinom(n = length(mu), mu = mu, size = mu/(phi-1)) #random negbin generator

#fit a generalized linear model y = f(x)  
glmQP = glm(y~x, family=quasipoisson) #quasipoisson
glmNB = glm.nb(y~x) #negative binomial

> glmQP

Call:  glm(formula = y ~ x, family = quasipoisson)

Coefficients:
(Intercept)            x  
    3.11257      0.01854  
(Dispersion parameter for quasipoisson family taken to be 3.613573)

Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
Null Deviance:      2097 
Residual Deviance: 356.8    AIC: NA

> glmNB

Call:  glm.nb(formula = y ~ x, init.theta = 23.36389741, link = log)

Coefficients:
(Intercept)            x  
    3.10182      0.01873  

Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
Null Deviance:      578.1 
Residual Deviance: 107.8    AIC: 824.7

Keduanya cocok mereproduksi parameter, dan quasipoisson memberikan estimasi 'masuk akal' untuk . Kita sekarang juga dapat mendefinisikan nilai AIC untuk quasipoisson:ϕ

df = 3 # three model parameters: a,b, and phi
phi.fit = 3.613573 #fitted phi value copied from summary(glmQP)
mu.fit = glmQP$fitted.values 

#dnbinom = negbin density, log=T returns log probabilities
AIC = 2*df - 2*sum(dnbinom(y, mu=mu.fit, size = mu.fit/(phi.fit - 1), log=T))
> AIC
[1] 819.329

(Aku harus secara manual menyalin dilengkapi nilai dari , karena saya tidak bisa menemukannya di objek)ϕsummary(glmQP)glmQP

Karena , ini akan menunjukkan bahwa quasipoisson, secara mengejutkan, lebih cocok; jadi setidaknya melakukan apa yang harus dilakukan, dan karenanya mungkin definisi yang masuk akal untuk AIC (dan dengan ekstensi, kemungkinan) dari quasipoisson. Pertanyaan besar yang tersisa bagi saya adalah saat ituSEBUAHsayaCQP<SEBUAHsayaCNBSEBUAHsayaCQP

  1. Apakah ide ini masuk akal? Apakah verifikasi saya didasarkan pada alasan sirkuler?
  2. Pertanyaan utama bagi siapa pun yang 'menciptakan' sesuatu yang tampaknya hilang dari topik mapan: jika ide ini masuk akal, mengapa itu belum diterapkan glm?

Edit: gambar ditambahkan

glm cocok dan + -1 sigma band

pengguna28400
sumber
1
(+1) Selamat Datang di Cross Divalidasi! Dan terima kasih atas pertanyaan yang sangat bagus (meskipun beberapa komentar dalam kode mungkin bagus untuk orang yang tidak menggunakan R). Saya pikir Anda mungkin telah menemukan kembali model NB1 (meskipun saya belum mengikutinya secara detail). Perhatikan juga bahwa tidak ada distribusi kuasi-Poisson - itulah sebabnya tidak ada kemungkinan atau AIC - itu hanya merujuk pada cara pemasangan sarana & varian.
Scortchi
2
Terima kasih! Saya telah menambahkan beberapa komentar untuk sementara waktu, saya harap semuanya beres. Saya mengerti bahwa distribusi quasi-Poisson tidak ada per se - apa yang saya benar-benar coba cari tahu adalah mengapa QP bahkan merupakan hal yang sama sekali, mengingat distribusi NB1 ada dan tidak memiliki masalah kuasi-QP apa pun (lihat jawaban Achims untuk resolusi yang jelas).
user28400
1
XPois(λ)Y=kXYμ=kλkμk10,k,2k,...
1
@ Glen_b: Apakah orang benar-benar menyebutnya quasi-Poisson? Bagaimanapun itu adalah ilustrasi yang bagus - ketika Anda menggunakan model "quasiPoisson" Anda tidak benar-benar mengasumsikan distribusi itu, atau NB1, atau lainnya, hanya hubungan antara mean & varians yang membuat perkiraan koefisien & kesalahan standar mereka lebih baik karena sampel semakin besar.
Scortchi
1
@Scortchi Ini adalah satu-satunya distribusi keluarga eksponensial yang memenuhi asumsi quasi-Poisson, jadi semacam - kadang-kadang saya melihat orang-orang menunjukkan bahwa itu adalah distribusi yang disiratkan oleh asumsi tersebut. Tentu saja ketika orang menggunakannya, mereka hampir * tidak pernah berniat bahwa data mereka berasal dari distribusi tertentu - itu hanya dimaksudkan sebagai deskripsi kasar tentang bagaimana makna dan varians mereka berhubungan. (Mungkin masuk akal di bawah asumsi yang sangat sederhana dalam beberapa aplikasi asuransi - total biaya klaim, di mana jumlah klaim adalah Poisson dan biaya per klaim secara efektif konstan.)
Glen_b -Reinstate Monica

Jawaban:

24

The quasi-Poisson bukan model kemungkinan maksimum penuh (ML) tetapi model kuasi-ML. Anda cukup menggunakan fungsi estimasi (atau fungsi skor) dari model Poisson untuk memperkirakan koefisien, dan kemudian menggunakan fungsi varians tertentu untuk mendapatkan kesalahan standar yang sesuai (atau lebih tepatnya matriks kovarians penuh) untuk melakukan inferensi. Oleh karena itu, glm()tidak menyediakan dan logLik()atau di AIC()sini, dll.

sizeθsayaμsaya

Jika tidak ada regressors (hanya intercept) yang parametrization NB1 dan parametrization NB2 dipekerjakan oleh MASS's glm.nb()bersamaan. Dengan regresi mereka berbeda. Dalam literatur statistik parametri NB2 lebih sering digunakan tetapi beberapa paket perangkat lunak juga menawarkan versi NB1. Misalnya dalam R, Anda dapat menggunakan gamlsspaket untuk melakukannya gamlss(y ~ x, family = NBII). Perhatikan bahwa gamlsspenggunaan membingungkan NBIuntuk parametri NB2 dan NBIIuntuk NB1. (Tetapi jargon dan terminologi tidak disatukan di semua komunitas.)

Maka Anda bisa bertanya, tentu saja, mengapa menggunakan quasi-Poisson jika ada NB1 tersedia? Masih ada perbedaan kecil: Yang pertama menggunakan kuasi-ML dan memperoleh estimasi dari dispersi dari residu kuadrat (atau Pearson). Yang terakhir menggunakan ML penuh. Dalam praktiknya, perbedaannya sering tidak besar tetapi motivasi untuk menggunakan kedua model sedikit berbeda.

Achim Zeileis
sumber
1
Terima kasih! Jawaban yang sangat membantu, saya bereksperimen dengan gamlsssekarang dan sepertinya itulah yang saya butuhkan. Bisakah Anda menguraikan motivasi untuk menggunakan kuasi-kemungkinan versus ML penuh?
user28400
2
Anda menganggap lebih sedikit: Anda hanya menganggap (1) hubungan log-linier antara harapan dan regresi (2) hubungan linear antara varians dan harapan. Kemungkinan lainnya tidak sepenuhnya ditentukan. Sebagai alternatif untuk (2), praktisi kadang-kadang menggunakan apa yang disebut "standar kesalahan sandwich" yang memungkinkan pola heteroskedastisitas yang lebih umum. Tentu saja, orang juga dapat menggunakan NB1 dengan kesalahan standar sandwich ... Beberapa komentar lagi ada di kami vignette("countreg", package = "pscl").
Achim Zeileis