Apakah "model rintangan" benar-benar satu model? Atau hanya dua model berurutan yang terpisah?

Pertimbangkan model rintangan yang memprediksi data jumlah ydari prediktor normal x:

set.seed(1839)
# simulate poisson with many zeros
x <- rnorm(100)
e <- rnorm(100)
y <- rpois(100, exp(-1.5 + x + e))

# how many zeroes?
table(y == 0)

FALSE  TRUE 
   31    69 

Dalam hal ini, saya memiliki data jumlah dengan 69 angka nol dan 31 jumlah positif. Nevermind untuk saat ini, menurut definisi prosedur pembuatan data, proses Poisson, karena pertanyaan saya adalah tentang model rintangan.

Katakanlah saya ingin menangani kelebihan nol ini dengan model rintangan. Dari bacaan saya tentang mereka, sepertinya model rintangan bukanlah model aktual — mereka hanya melakukan dua analisis berbeda secara berurutan. Pertama, regresi logistik yang memprediksi apakah nilainya positif versus nol. Kedua, regresi Poisson nol terpotong dengan hanya memasukkan kasus-kasus yang tidak nol. Langkah kedua ini terasa salah bagi saya karena (a) membuang data yang sangat baik, yang (b) dapat menyebabkan masalah daya karena banyak data nol, dan (c) pada dasarnya bukan "model" dalam dan dari dirinya sendiri , tetapi hanya menjalankan dua model secara berurutan.

Jadi saya mencoba "model rintangan" versus hanya menjalankan regresi Poisson logistik dan nol terpotong secara terpisah. Mereka memberi saya jawaban yang sama (saya menyingkat hasilnya, demi singkatnya):

> # hurdle output
> summary(pscl::hurdle(y ~ x))

Count model coefficients (truncated poisson with log link):
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)  -0.5182     0.3597  -1.441   0.1497  
x             0.7180     0.2834   2.533   0.0113 *

Zero hurdle model coefficients (binomial with logit link):
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -0.7772     0.2400  -3.238 0.001204 ** 
x             1.1173     0.2945   3.794 0.000148 ***

> # separate models output
> summary(VGAM::vglm(y[y > 0] ~ x[y > 0], family = pospoisson()))

Coefficients: 
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)  -0.5182     0.3597  -1.441   0.1497  
x[y > 0]      0.7180     0.2834   2.533   0.0113 *

> summary(glm(I(y == 0) ~ x, family = binomial))

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   0.7772     0.2400   3.238 0.001204 ** 
x            -1.1173     0.2945  -3.794 0.000148 ***
---

Ini nampak bagi saya karena banyak representasi matematis yang berbeda dari model termasuk probabilitas bahwa pengamatan adalah tidak nol dalam estimasi jumlah kasus positif, tetapi model yang saya jalankan di atas sepenuhnya mengabaikan satu sama lain. Sebagai contoh, ini dari Bab 5, halaman 128 dari Generalized Linear Models dari Smithson & Merkle untuk Variabel Ketergantungan Terbatas Kategori dan Kontinu :

... Kedua, probabilitas bahwa mengasumsikan nilai apa pun (nol dan bilangan bulat positif) harus sama dengan satu. Ini tidak dijamin dalam Persamaan (5.33). Untuk mengatasi masalah ini, kami mengalikan probabilitas Poisson dengan probabilitas keberhasilan Bernoulli .      Masalah ini mengharuskan kita untuk mengekspresikan model rintangan di atas sebagai mana , ,$y$$\pi$

$$P(Y=y|\boldsymbol{x,z,\beta,\gamma}) = \begin{cases} 1-\hat\pi &\text{for } y=0 \\ \hat\pi\times\frac{\exp(-\hat\lambda)\hat\lambda^y/y!}{1-\exp(-\hat\lambda)} &\text{for } y=1,2,\ldots \end{cases} \tag{5.34}$$ $\hat\lambda=\exp(\boldsymbol{x\beta})$$\hat\pi = {\rm logit}^{-1}(\boldsymbol{z\gamma})$$\boldsymbol x$adalah kovariat untuk model Poisson, adalah kovariat untuk model regresi logistik, dan dan adalah koefisien regresi masing-masing ... . $\boldsymbol z$$\hat{\boldsymbol{\beta}}$$\hat{\boldsymbol{\gamma}}$

Dengan melakukan dua model yang benar-benar terpisah satu sama lain — yang tampaknya merupakan model rintangan — saya tidak melihat bagaimana dimasukkan ke dalam prediksi kasus jumlah positif. Tetapi berdasarkan bagaimana saya dapat mereplikasi fungsi hanya dengan menjalankan dua model yang berbeda, saya tidak melihat bagaimana memainkan peran dalam Poisson terpotong regresi sama sekali.$\hat{\pi}$hurdle$\text{logit}^{-1}(z\hat{\gamma})$

Apakah saya memahami model rintangan dengan benar? Mereka tampaknya hanya menjalankan dua model berurutan: Pertama, logistik; Kedua, Poisson, benar-benar mengabaikan kasus di mana . Saya akan sangat menghargai jika seseorang dapat menjernihkan kebingungan saya dengan bisnis .$y = 0$$\hat{\pi}$

---

Jika saya benar bahwa itulah model rintangan, apa definisi dari model "rintangan", secara umum? Bayangkan dua skenario berbeda:

• Bayangkan memodelkan daya saing ras pemilu dengan melihat skor daya saing (1 - (proporsi suara pemenang - proporsi suara pemenang runner up)). Ini [0, 1), karena tidak ada ikatan (misalnya, 1). Model rintangan masuk akal di sini, karena ada satu proses (a) apakah pemilihan tidak dipertanyakan? dan (b) jika tidak, apa yang meramalkan daya saing? Jadi pertama-tama kami melakukan regresi logistik untuk menganalisis 0 vs (0, 1). Kemudian kami melakukan regresi beta untuk menganalisis (0, 1) kasus.

• Bayangkan sebuah studi psikologis yang khas. Responsnya adalah [1, 7], seperti skala Likert tradisional, dengan efek langit-langit yang sangat besar pada 7. Seseorang dapat melakukan model rintangan yang merupakan regresi logistik [1, 7) vs 7, dan kemudian regresi Tobit untuk semua kasus di mana tanggapan yang diamati adalah <7.

Apakah aman untuk menyebut kedua situasi ini model "rintangan" , bahkan jika saya memperkirakannya dengan dua model berurutan (logistik dan kemudian beta dalam kasus pertama, logistik dan kemudian Tobit dalam yang kedua)?

Memisahkan kemungkinan log

Memang benar bahwa sebagian besar model rintangan dapat diperkirakan secara terpisah (saya akan mengatakan, bukannya berurutan ). Alasannya adalah bahwa log-kemungkinan dapat didekomposisi menjadi dua bagian yang dapat dimaksimalkan secara terpisah. Ini karena hanyalah faktor penskalaan pada (5.34) yang menjadi istilah tambahan dalam log-likelihood.$\hat \pi$

Dalam notasi Smithson & Merkle: mana adalah kepadatan distribusi Poisson (yang tidak dikunci) dan adalah faktor dari pemotongan nol.

$$\begin{eqnarray*} \ell(\beta, \gamma; y, x, z) & = & \ell_1(\gamma; y, z) + \ell_2(\beta; y, x) \\ & = & \sum_{i: y_i = 0} \log\left\{1 - \mathrm{logit}^{-1}(z_i^\top \gamma)\right\} + \sum_{i: y_i > 0} \log\left\{\mathrm{logit}^{-1}(z_i^\top \gamma)\right\} + \\ & & \sum_{i: y_i > 0} \left[ \log \left\{f(y_i; \exp(x_i^\top \beta)\right\} - \log\left\{ 1 - f(0; \exp(x_i^\top \beta)\right\}\right] \end{eqnarray*}$$ $f(y; \lambda) = \exp(-\lambda) \lambda^y/y!$$1 - f(0; \lambda) = 1 - \exp(-\lambda)$

Maka menjadi jelas bahwa (model biner) dan (model Poisson nol-terpotong) dapat dimaksimalkan secara terpisah, yang mengarah ke estimasi parameter yang sama, kovarian, dll. Seperti dalam kasus ini di mana mereka dimaksimalkan bersama.$\ell_1(\gamma)$$\ell_2(\beta)$

Logika yang sama juga berfungsi jika probabilitas rintangan nol tidak ditentukan melalui model logit tetapi model regresi biner lainnya, misalnya, distribusi jumlah yang disensor dengan benar pada 1. Dan, tentu saja, juga dapat distribusi hitungan lain, misalnya binomial negatif. Seluruh pemisahan hanya rusak jika ada parameter bersama antara rintangan nol dan bagian hitung terpotong.$\pi$$f(\cdot)$

Contoh yang menonjol adalah jika distribusi binomial negatif dengan parameter terpisah tetapi umum digunakan dalam dua komponen model. (Ini tersedia di dalam paket dari R-Forge, penerus implementasi.)$\mu$$\theta$hurdle(..., separate = FALSE, dist = "negbin", zero.dist = "negbin")countregpscl

Pertanyaan konkret

(a) Membuang data yang sangat bagus: Dalam kasus Anda ya, secara umum tidak. Anda memiliki data dari model Poisson tunggal tanpa kelebihan nol (meskipun banyak nol ). Oleh karena itu, tidak perlu memperkirakan model terpisah untuk nol dan bukan nol. Namun, jika kedua bagian tersebut benar-benar didorong oleh parameter yang berbeda maka perlu untuk memperhitungkan ini.

(B) Dapat menyebabkan masalah listrik karena banyak data nol: Tidak harus. Di sini, Anda memiliki sepertiga dari pengamatan yang merupakan "keberhasilan" (rintangan rintangan). Ini tidak akan dianggap sangat ekstrim dalam model regresi biner. (Tentu saja, jika tidak perlu memperkirakan model yang terpisah Anda akan mendapatkan kekuatan.)

(c) Pada dasarnya bukan 'model' di dalam dan dari dirinya sendiri, tetapi hanya menjalankan dua model secara berurutan: Ini lebih filosofis dan saya tidak akan mencoba memberikan "satu" jawaban. Sebagai gantinya, saya akan menunjukkan sudut pandang pragmatis. Untuk estimasi model , akan lebih mudah untuk menekankan bahwa model terpisah karena - seperti yang Anda tunjukkan - Anda mungkin tidak memerlukan fungsi khusus untuk estimasi tersebut. Untuk aplikasi model , misalnya, untuk prediksi atau residu, dll., Akan lebih mudah untuk melihatnya sebagai model tunggal.

(D) Apakah aman untuk menyebut kedua model situasi ini 'rintangan': Pada prinsipnya ya. Namun, jargon dapat bervariasi di seluruh komunitas. Misalnya, regresi beta zero-hurdle lebih umum (dan sangat membingungkan) disebut regresi beta zero-inflated. Secara pribadi, saya menemukan yang terakhir sangat menyesatkan karena distribusi beta tidak memiliki nol yang dapat meningkat - tapi itu adalah istilah standar dalam literatur. Selain itu, model tobit adalah model yang disensor dan karenanya bukan model rintangan. Ini bisa diperluas, dengan model probit (atau logit) plus model normal terpotong . Dalam literatur ekonometrik ini dikenal sebagai model dua bagian Cragg.

Komentar perangkat lunak

The countregpaket pada R-Forge di https://R-Forge.R-project.org/R/?group_id=522 adalah implementasi penerus hurdle()/ zeroinfl()dari pscl. Alasan utama bahwa itu (masih) tidak pada CRAN adalah bahwa kami ingin merevisi predict()antarmuka, mungkin dengan cara yang tidak sepenuhnya kompatibel ke belakang. Kalau tidak, implementasinya cukup stabil. Dibandingkan dengan psclitu hadir dengan beberapa fitur bagus, misalnya:

• Sebuah zerotrunc()fungsi yang menggunakan kode yang sama persis seperti hurdle()untuk bagian nol-dipotong model. Dengan demikian, ia menawarkan alternatif VGAM.

• Selain itu, ini berfungsi sebagai d / p / q / r untuk distribusi penghitungan nol terpotong, rintangan, dan nol. Ini memfasilitasi melihat ini sebagai model "satu" daripada model terpisah.

• Untuk menilai goodness of fit, tampilan grafis seperti rootogram dan plot residu kuantil acak tersedia. (Lihat Kleiber & Zeileis, 2016, The American Statistician , 70 (3), 296–303. Doi: 10.1080 / 00031305.2016.1173590 .)

Data simulasi

Data simulasi Anda berasal dari proses Poisson tunggal. Jika ediperlakukan sebagai regresi yang dikenal maka itu akan menjadi Poisson GLM standar. Jika ekomponen noise yang tidak diketahui, maka ada beberapa heterogenitas yang tidak teramati yang menyebabkan sedikit penyebaran berlebihan yang dapat ditangkap oleh model binomial negatif atau beberapa jenis campuran kontinu atau efek acak dll. Namun, karena efeknya eagak kecil di sini , semua ini tidak membuat perbedaan besar. Di bawah, saya memperlakukan esebagai regressor (yaitu, dengan koefisien true 1) tetapi Anda juga bisa menghilangkan ini dan menggunakan model binomial atau Poisson negatif. Secara kualitatif, semua ini mengarah pada wawasan serupa.

## Poisson GLM
p <- glm(y ~ x + e, family = poisson)
## Hurdle Poisson (zero-truncated Poisson + right-censored Poisson)
library("countreg")
hp <- hurdle(y ~ x + e, dist = "poisson", zero.dist = "poisson")
## all coefficients very similar and close to true -1.5, 1, 1
cbind(coef(p), coef(hp, model = "zero"), coef(hp, model = "count"))
##                   [,1]       [,2]      [,3]
## (Intercept) -1.3371364 -1.2691271 -1.741320
## x            0.9118365  0.9791725  1.020992
## e            0.9598940  1.0192031  1.100175

Ini mencerminkan bahwa ketiga model dapat secara konsisten memperkirakan parameter sebenarnya. Melihat kesalahan standar yang sesuai menunjukkan bahwa dalam skenario ini (tanpa perlu bagian rintangan) Poisson GLM lebih efisien:

serr <- function(object, ...) sqrt(diag(vcov(object, ...)))
cbind(serr(p), serr(hp, model = "zero"), serr(hp, model = "count"))
##                  [,1]      [,2]      [,3]
## (Intercept) 0.2226027 0.2487211 0.5702826
## x           0.1594961 0.2340700 0.2853921
## e           0.1640422 0.2698122 0.2852902

Kriteria informasi standar akan memilih Poisson GLM sejati sebagai model terbaik:

AIC(p, hp)
##    df      AIC
## p   3 141.0473
## hp  6 145.9287

Dan tes Wald akan mendeteksi dengan benar bahwa dua komponen model rintangan tidak berbeda secara signifikan:

hurdletest(hp)
## Wald test for hurdle models
## 
## Restrictions:
## count_((Intercept) - zero_(Intercept) = 0
## count_x - zero_x = 0
## count_e - zero_e = 0
## 
## Model 1: restricted model
## Model 2: y ~ x + e
## 
##   Res.Df Df  Chisq Pr(>Chisq)
## 1     97                     
## 2     94  3 1.0562     0.7877

Akhirnya keduanya rootogram(p)dan qqrplot(p)menunjukkan bahwa Poisson GLM sangat cocok dengan data dan bahwa tidak ada kelebihan nol atau petunjuk tentang kesalahan spesifikasi lebih lanjut.

Saya setuju perbedaan antara model zero-inflated dan hurdle sulit dipahami. Keduanya adalah sejenis model campuran. Dari apa yang saya tahu, perbedaan penting adalah, dalam model nol-inflasi, Anda mencampur massa di nol dengan distribusi \ textit {yang juga bisa mengambil nilai nol}. Untuk model rintangan, Anda mencampur massa nol dengan distribusi yang hanya mengambil nilai lebih besar dari 0. Dengan demikian, dalam model nol Anda dapat membedakan antara 'nol struktural' (sesuai dengan massa nol) dan 'sampel nol 'sesuai dengan kemungkinan terjadinya 0 dari model yang Anda campur. Tentu saja, identifikasi ini sangat bergantung pada pilihan distribusi yang tepat! Tapi, jika Anda memiliki Poisson nol, misalnya, Anda dapat membedakan antara nol yang berasal dari komponen Poisson (nol pengambilan sampel) dan nol yang berasal dari massa nol (nol struktural). Jika Anda memiliki model inflasi nol dan distribusi yang Anda campur tidak memiliki massa nol, itu bisa ditafsirkan sebagai model rintangan.

Mengenai aspek filosofis, "kapan kita harus mempertimbangkan sesuatu model tunggal dan ketika dua model terpisah" , mungkin menarik untuk dicatat bahwa perkiraan sampel parameter-model berkorelasi.

Dalam plot di bawah ini dengan simulasi Anda sebagian besar melihat korelasi antara kemiringan dan intersep dari bagian hitungan. Tetapi ada juga sedikit hubungan antara bagian penghitungan dan bagian rintangan. Jika Anda mengubah parameter, misalnya membuat lambda dalam distribusi Poisson lebih kecil atau ukuran sampel lebih kecil, maka korelasinya menjadi lebih kuat.

Jadi saya akan mengatakan bahwa Anda tidak harus menganggapnya sebagai dua model terpisah . Atau setidaknya ada beberapa hubungan bahkan ketika dalam praktiknya Anda dapat menghitung dua perkiraan yang independen satu sama lain.

set.seed(1839)

Nrep <- 3000
Ns <- 100
pars <- matrix(rep(0,3*Nrep),Nrep)
colnames(pars) <- c("count_intercept","count_slope","hurdle_intercept")

# simulation-loop
# Note that a truncated poisson is used to generate data
# this will make the parameters from the hurdle function easier to interpret and compare
for (i in 1:Nrep) {
  x <- rnorm(Ns,0,1)
  e <- rbinom(Ns,1,exp(-0.7))
  y <- e*truncdist::rtrunc(n=Ns,spec='pois',a=0,b=Inf,lambda=exp(-1.5 + x))
  mod <- pscl::hurdle(y ~ 1+x|1, link="log")
  pars[i,1]<-mod$coefficients$count[1]
  pars[i,2]<-mod$coefficients$count[2]
  pars[i,3]<-mod$coefficients$zero[1]
}  

# viewing data
plotpars <- pars[pars[,1]>-7,] #clipping
pairs(plotpars,cex=0.7,pch=21,
      col= rgb(0,0,0,0.03),
      bg = rgb(0,0,0,0.03))

# demonstrating linear relation / significant correlation
summary(lm(pars[,1] ~ pars[,3]))

---

Tidak masuk akal jika ada korelasi antara dua bagian. Tapi itu mungkin karena tingkat diskrit estimasi untuk parameter dalam model Poisson, dan bagaimana ini berhubungan dengan jumlah nol.