Bisakah Bobot dan Offset mengarah ke hasil serupa dalam regresi poisson?

Dalam "Panduan Practioner untuk Generalized linear models" dalam paragraf 1.83 dinyatakan bahwa:

"Dalam kasus khusus GLM multiplikatif Poisson dapat ditunjukkan bahwa klaim pemodelan dihitung dengan istilah offset yang sama dengan log paparan yang menghasilkan hasil yang identik dengan pemodelan frekuensi klaim dengan bobot sebelumnya yang ditetapkan sama dengan paparan setiap pengamatan. "

Saya tidak dapat menemukan referensi lebih lanjut dari hasil ini, jadi saya mengambil beberapa pengujian empiris di mana saya tidak dapat menemukan bukti bahwa pernyataan itu benar. Adakah yang bisa memberikan beberapa wawasan mengapa hasil ini mungkin benar / salah.

FYI, saya menggunakan kode R berikut untuk menguji hipotesis, di mana saya tidak bisa mendapatkan hasil yang serupa untuk dua kasus yang disebutkan:

n=1000
m=10

# Generate random data
X = matrix(data = rnorm(n*m)+1, ncol = m, nrow = n)

intercept = 2
coefs = runif(m)
offset = runif(n)
## DGP: exp of Intercept + linear combination X variables + log(offset)
mu = exp(intercept + X%*%coefs + log(offset))
y = rpois(n=n, lambda=mu)

df = data.frame('y'=y, 'X'=X, 'offset' = offset)
formula = paste("y ~",paste(colnames(df)[grepl("X", colnames(df))], collapse = "+"))

#First model using log(offset) as offset
fit1  = glm(formula, family = "poisson", df, offset = log(offset))
#Second model using offset as weights for individual observations
fit2 = glm(formula, family = "poisson", df, weights = offset) 
#Third model using poisson model on y/offset as reference
dfNew = df
dfNew$y = dfNew$y/offset
fit3 = glm(formula, family = "poisson", dfNew)

#Combine coefficients with the true coefficients
rbind(fit1$coefficients, fit2$coefficients, fit3$coefficients, c(intercept,coefs))

Estimasi koefisien yang dihasilkan dari menjalankan kode ini diberikan di bawah:

 >  
    (Intercept)       X.1       X.2       X.3        X.4       X.5       X.6
[1,]    1.998277 0.2923091 0.4586666 0.1802960 0.11688860 0.7997154 0.4786655
[2,]    1.588620 0.2708272 0.4540180 0.1901753 0.07284985 0.7928951 0.5100480
[3,]    1.983903 0.2942196 0.4593369 0.1782187 0.11846876 0.8018315 0.4807802
[4,]    2.000000 0.2909240 0.4576965 0.1807591 0.11658183 0.8005451 0.4780123
              X.7       X.8       X.9      X.10
[1,]  0.005772078 0.9154808 0.9078758 0.3512824
[2,] -0.003705015 0.9117014 0.9063845 0.4155601
[3,]  0.007595660 0.9181014 0.9076908 0.3505173
[4,]  0.005881960 0.9150350 0.9084375 0.3511749
>

dan kita bisa amati koefisiennya tidak identik.

generalized-linear-model modeling poisson-regression weights offset BDP1
sumber

Anda tidak harus memasukkan rm(list=ls() )kode R yang Anda posting di sini! Itu bisa mengejutkan seseorang yang menjalankannya, membuat mereka marah pada Anda. Saya menghapusnya. Saya juga diedit untuk memasukkan hasil dari menjalankan kode. Jika Anda telah melakukan ini pada awalnya, mungkin Anda mendapatkan respons lebih cepat, karena hanya sedikit pembaca yang menjalankan kode itu sendiri.

kjetil b halvorsen

@kjetilbhalvorsen Terima kasih atas umpan baliknya, akan dilakukan di masa mendatang.

BDP1

@ BDP1 Seperti sekarang, kode R Anda tidak menguji klaim yang Anda maksud. Saya sarankan Anda untuk menambahkan istilah berat untuk fit3 dan menambahkan addendum langsung ke pertanyaan.

aivanov

Jawaban:

(dengan kode R, Anda dapat mengganti "poisson" dengan "quasipoisson" untuk menghindari semua peringatan yang dihasilkan. Tidak ada lagi impor yang akan berubah. Lihat (*) di bawah). Referensi Anda menggunakan istilah "multiplicative glm" yang menurut saya hanya berarti glm dengan tautan log, karena tautan log dapat dianggap sebagai model multiplikatif. Contoh Anda sendiri menunjukkan bahwa klaim itu salah, setidaknya seperti yang kami tafsirkan (Karena parameter yang diestimasi tidak sama). Anda bisa menulis kepada penulis dan bertanya apa artinya. Di bawah ini saya akan berdebat mengapa klaim itu salah.

Membiarkan $\lambda_i$ menjadi parameter poisson dan $\omega_i$ bobot. Membiarkan $\eta_i$ menjadi prediktor linier tanpa offset, dan kemudian $\eta_i+\log(\omega_i)$ menjadi prediktor linier dengan offset. Fungsi probabilitas poisson adalah

f (y_{i}) = e^{- λ_{i}} λ_{i}^{y_{i}} / y_{i}!

$f(y_i) = e^{-\lambda_i} \lambda_i^{y_i}/y_i !$ Kemudian fungsi log likelihood untuk model dengan offset menjadi

ℓ = - \sum_{i} ω_{i} e^{η_{i}} + \sum_{i} y_{i} η_{i} + \sum_{i} y_{i} \log ω_{i} - \sum_{i} \log y_{i}!

$\ell = -\sum_i \omega_i e^{\eta_i} + \sum_i y_i \eta_i +\sum_i y_i\log \omega_i - \sum_i \log y_i!$ sedangkan fungsi log likelihood untuk model dengan bobot menjadi

ℓ^{w} = - \sum_{i} ω_{i} e^{η_{i}} + \sum_{i} y_{i} ω_{i} η_{i} - \sum_{i} ω_{i} \log y_{i}!

$\ell^w = -\sum_i \omega_i e^{\eta_i}+\sum_i y_i \omega_i \eta_i -\sum_i \omega_i \log y_i!$ dan ini jelas tidak sama. Jadi apa yang penulis maksudkan itu tidak jelas bagi saya.

(*) Catatan dari bantuan glmfungsi R :

Bobot 'NULL' 'dapat digunakan untuk menunjukkan bahwa pengamatan yang berbeda memiliki dispersi yang berbeda (dengan nilai dalam' bobot 'berbanding terbalik dengan dispersi); atau setara, ketika elemen 'bobot' adalah bilangan bulat positif w_i, bahwa setiap respons y_i adalah rata-rata dari pengamatan berat unit w_i. Untuk GLM binomial, bobot sebelum digunakan untuk memberikan jumlah percobaan ketika responsnya adalah proporsi keberhasilan: mereka jarang akan digunakan untuk GLM Poisson.

Melihat makna argumen bobot menjelaskan ini, ia memberi sedikit makna dengan fungsi keluarga poisson, yang mengasumsikan parameter skala konstan $\phi=1$ sementara argumen bobot memodifikasi $\phi$ . Ini memang memberi lebih banyak makna dengan fungsi keluarga quasiposson. Lihat jawaban untuk input "bobot" dalam fungsi glm dan lm di R. Jawaban yang diberikan di sana juga membantu dalam melihat mengapa kemungkinan dalam kasus berbobot mengambil bentuk yang diberikan di atas.

Jawaban yang diberikan di sini mungkin relevan: Bagaimana regresi tingkat Poisson sama dengan regresi Poisson dengan istilah offset yang sesuai? dan sangat menarik.

kjetil b halvorsen
sumber

Terima kasih atas jawabannya. Akan melalui kontribusi kemungkinan banyak menjelaskan. Setelah mencari beberapa jawaban untuk jawaban Anda, saya menemukan pertanyaan berikut: stats.stackexchange.com/questions/66791/... Di mana ditunjukkan bahwa dengan membagi respons asli dengan offset, dan mengatur variabel offset sebagai bobot, hasil yang sama dapat diperoleh sebagai model dasar, di mana log (offset) memasuki prediktor linier dengan koefisien konstan 1.

BDP1

Saya juga berusaha untuk mendapatkan bahwa kontribusi kemungkinan dari model baru ini sama dengan kontribusi dari model yang baru disebutkan, tetapi tidak dapat melakukannya, mengingat

(y_{i} / w_{i})!

$(y_i/w_i)!$ istilah yang muncul untuk yang terakhir.

BDP1

Setidaknya dalam percobaan di R-script yang disediakan, klaim ini tampaknya benar. ini dapat dengan mudah dilihat dengan menambahkan bobot = argumen offset pada baris fit3 = glm (...)

BDP1

Saya tidak mengerti apa yang Anda katakan di sini, jika menurut Anda eksperimen Anda menguatkan, mengapa tidak mengedit perubahan itu pada percobaan ke OP dan menjelaskan mengapa menurut Anda eksperimen itu menguatkan klaim. Saya mencoba menjelaskan mengapa itu tidak benar, apa yang salah dengan argumen saya?

kjetil b halvorsen

Kutipan dari Panduan Practioner's sebenarnya benar, hanya saja belum diterapkan dengan benar oleh OP. Ini telah ditunjukkan oleh Alan Chalk dalam jawaban lain di sini, oleh saya stats.stackexchange.com/questions/430001 dan juga oleh Cokes stats.stackexchange.com/questions/264071 (meskipun kesimpulan yang benar dikubur di baris terakhir dari jawaban itu).

Gordon Smyth

Maaf tidak hanya menambahkan komentar di atas, tapi saya tidak punya cukup perwakilan.

Klaim asli - tetapi dimodifikasi sedikit - ternyata benar.

Dua model berikut memberikan jawaban yang persis sama dalam R menggunakan poisson glm dengan log-link:

model y, gunakan offset yang sedang log (offset)
model y / offset, gunakan bobot sama dengan offset

menyesuaikan kode asli Anda sedikit menunjukkan jawaban yang identik:

n=1000
m=10

# Generate random data
X = matrix(data = rnorm(n*m)+1, ncol = m, nrow = n)

intercept = 2
coefs = runif(m)
offset = runif(n)
## DGP: exp of Intercept + linear combination X variables + log(offset)
mu = exp(intercept + X%*%coefs + log(offset))
y = rpois(n=n, lambda=mu)

df = data.frame('y' = y,
                'y_over_offset' = y/offset,
                'X' = X,
                'offset' = offset)

formula_offset = paste("y ~",paste(colnames(df)[grepl("X", colnames(df))], collapse = "+"))
formula_weights = paste("y_over_offset ~",paste(colnames(df)[grepl("X", colnames(df))], collapse = "+"))

#First model using log(offset) as offset
fit1  = glm(formula_offset, family = "poisson", df, offset = log(offset))
#Second model using offset as weights for individual observations
fit2 = glm(formula_weights, family = "poisson", df, weights = offset) 


#Combine coefficients with the true coefficients
rbind(fit1$coefficients, fit2$coefficients, c(intercept,coefs))

Semoga itu memberikan jawaban yang identik.

Dimungkinkan untuk menunjukkan bahwa kedua model secara statistik setara (ada kertas CAS di suatu tempat yang menunjukkan ini - saya akan memposting tautan jika saya punya waktu).

Kebetulan, jika Anda melakukan regresi penalti maka cara paket yang berbeda seperti glmnet dan H2o mengukur penyimpangan untuk dua cara yang berbeda dalam mendefinisikan model, dapat menyebabkan hasil yang berbeda.

Alan Chalk
sumber

Pertanyaan singkat Anda dapat melihat bahwa dengan fit2 Anda tidak memiliki AIC dan plot antara 2 cocok adalah plot yang sedikit berbeda (fit1); plot (fit2) - juga tahukah Anda mengapa ini terjadi?]

Charl Francois Marais