Kapan Poisson dan regresi binomial negatif sesuai dengan koefisien yang sama?

Saya perhatikan bahwa dalam regresi R, Poisson dan negatif binomial (NB) selalu cocok dengan koefisien yang sama untuk prediktor kategoris, tetapi tidak kontinu.

Misalnya, berikut ini adalah regresi dengan prediktor kategori:

data(warpbreaks)
library(MASS)

rs1 = glm(breaks ~ tension, data=warpbreaks, family="poisson")
rs2 = glm.nb(breaks ~ tension, data=warpbreaks)

#compare coefficients
cbind("Poisson"=coef(rs1), "NB"=coef(rs2))

Berikut ini adalah contoh dengan prediktor kontinu, di mana Poisson dan NB sesuai dengan koefisien yang berbeda:

data(cars)
rs1 = glm(dist ~ speed, data=cars, family="poisson")
rs2 = glm.nb(dist ~ speed, data=cars)

#compare coefficients
cbind("Poisson"=coef(rs1), "NB"=coef(rs2))

(Tentu saja ini tidak menghitung data, dan modelnya tidak berarti ...)

Lalu saya mengkode ulang prediktor menjadi faktor, dan kedua model tersebut kembali sesuai dengan koefisien yang sama:

library(Hmisc)
speedCat = cut2(cars$speed, g=5) 
#you can change g to get a different number of bins

rs1 = glm(cars$dist ~ speedCat, family="poisson")
rs2 = glm.nb(cars$dist ~ speedCat)

#compare coefficients
cbind("Poisson"=coef(rs1), "NB"=coef(rs2))

Namun, Regresi Binomial Negatif Joseph Hilbe memberikan contoh (6.3, hal 118-119) di mana prediktor kategori, jenis kelamin, cocok dengan koefisien yang sedikit berbeda oleh Poisson ( $b=0.883$ ) dan NB ( $b=0.881$ ). Dia mengatakan: “Rasio tingkat kejadian antara model Poisson dan NB sangat mirip. Ini tidak mengejutkan mengingat kedekatan $\alpha$ [sesuai dengan $1/\theta$ dalam R] ke nol. ”

Saya memahami ini, tetapi dalam contoh di atas, summary(rs2)memberi tahu kita bahwa diperkirakan masing-masing 9,16 dan 7,93.$\theta$

Jadi mengapa koefisiennya persis sama? Dan mengapa hanya untuk prediksi kategori?

---

Edit # 1

Berikut adalah contoh dengan dua prediktor non-ortogonal. Memang, koefisien tidak lagi sama:

data(cars)

#make random categorical predictor
set.seed(1); randomCats1 = sample( c("A","B","C"), length(cars$dist), replace=T)
set.seed(2); randomCats2 = sample( c("C","D","E"), length(cars$dist), replace=T)

rs1 = glm(dist ~ randomCats1 + randomCats2, data=cars, family="poisson")
rs2 = glm.nb(dist ~ randomCats1 + randomCats2, data=cars)

#compare coefficients
cbind("Poisson"=coef(rs1), "NB"=coef(rs2))

Dan, termasuk prediktor lain menyebabkan model sesuai dengan koefisien yang berbeda bahkan ketika prediktor baru kontinu. Jadi, itu ada hubungannya dengan ortogonalitas dari variabel dummy yang saya buat dalam contoh asli saya?

rs1 = glm(dist ~ randomCats1 + speed, data=cars, family="poisson")
rs2 = glm.nb(dist ~ randomCats1 + speed, data=cars)

#compare coefficients
cbind("Poisson"=coef(rs1), "NB"=coef(rs2))

Anda telah menemukan properti GLM yang intim, tetapi umum, cocok dengan kemungkinan maksimum . Hasilnya keluar begitu seseorang mempertimbangkan kasus paling sederhana dari semuanya: Menyesuaikan parameter tunggal ke tunggal pengamatan!

Satu kalimat jawaban : Jika semua kita peduli sepatutnya sarana terpisah untuk subset menguraikan sampel kami, maka GLMS akan selalu menghasilkan μ j = ˉ y j untuk setiap bagian j , sehingga struktur kesalahan yang sebenarnya dan parametrization kepadatan kedua menjadi tidak relevan ke estimasi (titik)!$\hat\mu_j = \bar y_j$$j$

Sedikit lebih banyak : Menyesuaikan faktor-faktor kortogonal ortogonal dengan kemungkinan maksimum setara dengan cara-cara yang terpisah untuk memisahkan subset dari sampel kami, jadi ini menjelaskan mengapa Poisson dan GLM binomial negatif menghasilkan estimasi parameter yang sama. Memang, hal yang sama terjadi apakah kita menggunakan Poisson, negbin, Gaussian, invers Gaussian atau Regresi Gamma (lihat di bawah). Dalam kasus Poisson dan negbin, fungsi tautan default adalah tautan , tetapi itu adalah herring merah; sementara ini menghasilkan mentah yang sama$\log$ estimasi parameter , kita akan melihat di bawah bahwa properti ini benar-benar tidak ada hubungannya dengan fungsi tautan sama sekali.

Ketika kita tertarik pada parametrization dengan lebih banyak struktur, atau itu tergantung pada prediktor kontinu, maka struktur kesalahan yang diasumsikan menjadi relevan karena hubungan mean-varians dari distribusi karena berkaitan dengan parameter dan fungsi nonlinear yang digunakan untuk memodelkan kondisi. cara.

GLM dan keluarga dispersi eksponensial: Crash course

Sebuah keluarga dispersi eksponensial dalam bentuk alami adalah salah satu seperti bahwa kepadatan log adalah bentuk

$$\log f(y;\,\theta,\nu) = \frac{\theta y - b(\theta)}{\nu} + a(y,\nu) \>.$$

Di sini adalah parameter alami dan ν adalah parameter dispersi . Jika ν diketahui, ini hanya akan menjadi keluarga eksponensial satu parameter standar. Semua GLM yang dipertimbangkan di bawah ini menganggap model kesalahan dari keluarga ini.$\theta$$\nu$$\nu$

Pertimbangkan contoh pengamatan tunggal dari keluarga ini. Jika kita cocok oleh kemungkinan maksimum, kita mendapatkan bahwa y = b ' ( θ ) , terlepas dari nilai ν . Ini mudah meluas ke kasus sampel iid sejak likelihood log menambahkan, menghasilkan ˉ y = b ' ( θ ) .$\theta$$y = b'(\hat\theta)$$\nu$$\bar y = b'(\hat\theta)$

Tetapi, kita juga tahu, karena keteraturan yang bagus dari kepadatan log sebagai fungsi dari , itu $\theta$ Jadi, sebenarnya b ′ ( θ ) = E Y = μ .

$$\frac{\partial}{\partial \theta} \mathbb E \log f(Y;\theta,\nu) = \mathbb E \frac{\partial}{\partial \theta} \log f(Y;\theta,\nu) = 0 \>.$$ $b'(\theta) = \mathbb E Y = \mu$

Karena perkiraan kemungkinan maksimum yang berubah dalam transformasi, ini berarti bahwa untuk keluarga ini kepadatan.$\bar y = \hat\mu$

Sekarang, dalam GLM, kami memodelkan sebagai μ i = g - 1 ( x T i β ) di mana g adalah fungsi tautan. Tetapi jika x i adalah vektor semua nol kecuali untuk 1 tunggal pada posisi j , maka μ i = g ( β j ) . Kemungkinan GLM kemudian memfaktorkan berdasarkan β j 'dan kami melanjutkan seperti di atas. Inilah tepatnya kasus faktor ortogonal.$\mu_i$$\mu_i = g^{-1}(\mathbf x_i^T \beta)$$g$$\mathbf x_i$$j$$\mu_i = g(\beta_j)$$\beta_j$

Apa yang berbeda dari prediktor kontinu?

Ketika prediktor kontinu atau kategorik, tetapi tidak dapat direduksi menjadi bentuk ortogonal, maka kemungkinan tidak lagi menjadi faktor individu dengan rata-rata terpisah tergantung pada parameter terpisah. Pada titik ini, struktur kesalahan dan fungsi tautan memang ikut berperan.

Jika salah satu engkol melalui (membosankan) aljabar, persamaan kemungkinan menjadi untuk semua j = 1 , ... , p di mana λ i = x T i β . Di sini,parameter β dan ν masuk secara implisit melalui hubungan tautan μ i = g ( λ i ) = g ( x T i β ) dan varians σ 2 i .

$$\sum_{i=1}^n \frac{(y_i - \mu_i)x_{ij}}{\sigma_i^2}\frac{\partial \mu_i}{\partial \lambda_i} = 0\>,$$ $j = 1,\ldots,p$$\lambda_i = \mathbf x_i^T \beta$$\beta$$\nu$$\mu_i = g(\lambda_i) = g(\mathbf x_i^T \beta)$$\sigma_i^2$

Dengan cara ini, fungsi tautan dan model kesalahan yang diasumsikan menjadi relevan dengan estimasi.

Contoh: Model kesalahan (hampir) tidak masalah

Dalam contoh di bawah ini, kami menghasilkan data acak binomial negatif tergantung pada tiga faktor kategori. Setiap pengamatan berasal dari satu kategori dan parameter dispersi yang sama ( ) digunakan.$k = 6$

Kami kemudian cocok dengan data ini menggunakan lima GLM yang berbeda, masing-masing dengan tautan : ( a ) binomial negatif, ( b ) Poisson, ( c ) Gaussian, ( d ) Inverse Gaussian dan ( e ) Gamma GLMs. Semua ini adalah contoh keluarga dispersi eksponensial.$\log$

Dari tabel, kita dapat melihat bahwa estimasi parameter identik , meskipun beberapa GLM ini adalah untuk data diskrit dan yang lain untuk kontinu, dan beberapa untuk data non-negatif sementara yang lain tidak.

      negbin  poisson gaussian invgauss    gamma
XX1 4.234107 4.234107 4.234107 4.234107 4.234107
XX2 4.790820 4.790820 4.790820 4.790820 4.790820
XX3 4.841033 4.841033 4.841033 4.841033 4.841033

Peringatan dalam pos berasal dari fakta bahwa prosedur pemasangan akan gagal jika pengamatan tidak termasuk dalam domain kepadatan tertentu. Misalnya, jika kami memiliki hitungan yang dihasilkan secara acak dalam data di atas, maka Gamma GLM akan gagal untuk menyatu karena GLM Gamma membutuhkan data yang benar-benar positif.$0$

Contoh: Fungsi tautan (hampir) tidak masalah

$\log$$\log(\hat \beta)$$\log(\hat \beta^2)$$\hat\beta$

> coefs.po
         log       id     sqrt
XX1 4.234107 4.234107 4.234107
XX2 4.790820 4.790820 4.790820
XX3 4.841033 4.841033 4.841033

Peringatan dalam pos hanya merujuk pada fakta bahwa taksiran mentah akan berbeda dengan fungsi tautan, tetapi taksiran taksiran parameter rata-rata tidak akan.

Kode r

# Warning! This code is a bit simplified for compactness.
library(MASS)
n <- 5
m <- 3
set.seed(17)
b <- exp(5+rnorm(m))
k <- 6

# Random negbin data; orthogonal factors
y <- rnbinom(m*n, size=k, mu=rep(b,each=n))
X <- factor(paste("X",rep(1:m,each=n),sep=""))

# Fit a bunch of GLMs with a log link
con <- glm.control(maxit=100)
mnb <- glm(y~X+0, family=negative.binomial(theta=2))
mpo <- glm(y~X+0, family="poisson")
mga <- glm(y~X+0, family=gaussian(link=log), start=rep(1,m), control=con)
miv <- glm(y~X+0, family=inverse.gaussian(link=log), start=rep(2,m), control=con)
mgm <- glm(y~X+0, family=Gamma(link=log), start=rep(1,m), control=con)    
coefs <- cbind(negbin=mnb$coef, poisson=mpo$coef, gaussian=mga$coef
                   invgauss=miv$coef, gamma=mgm$coef)

# Fit a bunch of Poisson GLMs with different links.
mpo.log  <- glm(y~X+0, family=poisson(link="log"))
mpo.id   <- glm(y~X+0, family=poisson(link="identity"))
mpo.sqrt <- glm(y~X+0, family=poisson(link="sqrt"))   
coefs.po <- cbind(log=mpo$coef, id=log(mpo.id$coef), sqrt=log(mpo.sqrt$coef^2))

Untuk melihat apa yang terjadi di sini, pertama-tama berguna untuk melakukan regresi tanpa intersep, karena intersep dalam regresi kategori dengan hanya satu prediktor tidak ada artinya:

> rs1 = glm(breaks ~ tension-1, data=warpbreaks, family="poisson")
> rs2 = glm.nb(breaks ~ tension-1, data=warpbreaks)

Karena Poisson dan regresi binomial negatif menentukan log dari parameter rata-rata, maka untuk regresi kategoris, dengan mengekspresiasikan koefisien akan memberikan Anda parameter rata-rata aktual untuk setiap kategori:

>  exp(cbind("Poisson"=coef(rs1), "NB"=coef(rs2)))
          Poisson       NB
tensionL 36.38889 36.38889
tensionM 26.38889 26.38889
tensionH 21.66667 21.66667

Parameter ini sesuai dengan rata-rata aktual pada nilai kategori yang berbeda:

> with(warpbreaks,tapply(breaks,tension,mean))
       L        M        H 
36.38889 26.38889 21.66667 

$\lambda$

$\lambda$

$\lambda$$\theta$$\theta$$\lambda$

$$\begin{align} L(X,\lambda,\theta) &= \prod \left(\frac{\theta}{\lambda+\theta}\right)^\theta \frac{\Gamma(\theta + x_i)}{x_i!\Gamma(\theta)}\left(\frac{\lambda}{\lambda+\theta}\right)^{x_i}\\ \log L(X,\lambda,\theta) &= \sum\theta\left(\text{log}\theta-\text{log}(\lambda+\theta)\right) +x_i\left(\text{log}\lambda-\text{log}(\lambda+\theta)\right) +\log\left(\frac{\Gamma(\theta + x_i)}{x_i!\Gamma(\theta)}\right)\\ \frac{d}{d\lambda}\log L(X,\lambda,\theta) &= \sum \frac{x_i}{\lambda}-\frac{\theta+x_i}{\lambda+\theta} =n\left(\frac{\bar{x}}{\lambda}-\frac{\bar{x}+\theta}{\lambda+\theta}\right), \end{align}$$ so the maximum is attained when $\lambda=\bar{x}$.

The reason that you don't get the same coefficients for continuous data is because in a continuous regression, $\text{log}(\lambda)$ is no longer going to be a piecewise constant function of the predictor variables, but a linear one. Maximizing the likelihood function in this case will not reduce to independently fitting a value $\lambda$ for disjoint subsets of the data, but will rather be a nontrivial problem that is solved numerically, and is likely to produce different results for different likelihood functions.

Similarly, if you have multiple categorical predictors, despite the fact that the fitted model will ultimately specify $\lambda$ as a piecewise constant function, in general there will not be enough degrees of freedom to allow $\lambda$ to be determined independently for each constant segment. For example, suppose that you have $2$ predictors with $5$ categories each. In this case, your model has $10$ degrees of freedom, whereas there are $5*5=25$ unique different combinations of the categories, each of which will have its own fitted value of $\lambda$. So, assuming that the intersections of these categories are non-empty (or at least that $11$ of them are nonempty), the likelihood maximization problem again becomes nontrivial and will generally produce different outcomes for Poisson versus negative binomial or any other distribution.