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
yang dapat dipasang di R menggunakan glm(..,family=poisson)
dan glm.nb(...)
, masing-masing. Ada juga quasipoisson
keluarga, yang dalam pemahaman saya adalah Poisson yang disesuaikan dengan EV dan varians yang sama
,
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:
,
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 itu
- Apakah ide ini masuk akal? Apakah verifikasi saya didasarkan pada alasan sirkuler?
- 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
sumber
Jawaban:
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 danlogLik()
atau diAIC()
sini, dll.size
Jika tidak ada regressors (hanya intercept) yang parametrization NB1 dan parametrization NB2 dipekerjakan oleh
MASS
'sglm.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 menggunakangamlss
paket untuk melakukannyagamlss(y ~ x, family = NBII)
. Perhatikan bahwagamlss
penggunaan membingungkanNBI
untuk parametri NB2 danNBII
untuk 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.
sumber
gamlss
sekarang dan sepertinya itulah yang saya butuhkan. Bisakah Anda menguraikan motivasi untuk menggunakan kuasi-kemungkinan versus ML penuh?vignette("countreg", package = "pscl")
.