R

“跨欄模型”真的是一種模型嗎?或者只是兩個獨立的順序模型?

  • December 30, 2017

y考慮一個從正常預測器預測計數數據的障礙模型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 

在這種情況下,我有 69 個零和 31 個正計數的計數數據。暫時不要介意,根據數據生成過程的定義,這是一個泊松過程,因為我的問題是關於障礙模型的。

假設我想通過障礙模型處理這些多餘的零。從我對它們的閱讀來看,障礙模型本身似乎並不是真正的模型——它們只是按順序進行兩種不同的分析。首先,邏輯回歸預測該值是否為正與零。其次,零截斷泊松回歸*僅包括非零情況。*我覺得這第二步是錯誤的,因為它 (a) 丟棄了非常好的數據,這 (b) 可能導致電源問題,因為大部分數據都是零,並且 (c) 本身基本上不是“模型” ,但只是順序運行兩個不同的模型。

所以我嘗試了一個“障礙模型”,而不是單獨運行邏輯和零截斷泊松回歸。他們給了我相同的答案(為簡潔起見,我將輸出縮寫):

> # 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 ***
---

這對我來說似乎是因為模型的許多不同的數學表示包括在正計數情況的估計中觀察值非零的概率,但我上面運行的模型*完全忽略了另一個。*例如,這是來自 Smithson & Merkle’s Generalized Linear Models for Categorical and Continuous Limited Dependent Variables的第 128 頁第 5 章:

…第二,概率假定任何值(零和正整數)必須等於一。這在公式 (5.33) 中無法保證。為了解決這個問題,我們將泊松概率乘以伯努利成功概率.

     這些問題需要我們將上述障礙模型表示為

在哪裡,,是泊松模型的協變量,是邏輯回歸模型的協變量,並且和是各自的回歸係數….

通過將兩個模型彼此完全分開——這似乎是障礙模型所做的——我不明白如何被納入陽性計數病例的預測中。但是基於我如何hurdle通過運行兩個不同的模型來複製該功能,我不知道如何完全在截斷泊松回歸中起作用。

我是否正確理解障礙模型?他們似乎有兩個只是運行兩個順序模型:首先,一個邏輯;第二,泊松,完全忽略的情況. 如果有人能解決我對商業。


如果我是正確的,那就是障礙模型,更一般地說,“障礙”模型的定義是什麼?想像兩種不同的場景:

  • 想像一下通過查看競爭力分數(1 - (獲勝者的投票比例 - 亞軍的投票比例))來模擬選舉競賽的競爭力。這是 [0, 1),因為沒有關係(例如,1)。障礙模型在這裡是有意義的,因為有一個過程 (a) 選舉沒有競爭?(b) 如果不是,那麼預測的競爭力是什麼?所以我們首先做一個邏輯回歸來分析 0 vs. (0, 1)。然後我們做 beta 回歸來分析 (0, 1) 的情況。
  • 想像一個典型的心理學研究。響應是 [1, 7],就像傳統的李克特量表一樣,在 7 處具有巨大的天花板效應。人們可以做一個障礙模型,即 [1, 7) 與 7 的邏輯回歸,然後對所有情況下的 Tobit 回歸觀察到的反應 < 7。

將這兩種情況稱為“障礙”模型是否安全,即使我使用兩個順序模型(邏輯模型,然後是第一種情況下的 beta,後勤情況,然後是第二種情況下的 Tobit)來估計它們?

分離對數似然

大多數障礙模型可以單獨估計是正確的(我會說,而不是順序)。原因是對數似然可以分解為可以分別最大化的兩部分。這是因為只是 (5.34) 中的一個比例因子,它成為對數似然中的一個附加項。

在史密森和默克爾的符號中:

在哪裡是(未截斷的)泊松分佈的密度,並且是從零截斷的因子。 那麼很明顯(二進制 logit 模型)和(零截斷泊松模型)可以單獨最大化,導致與聯合最大化的情況下相同的參數估計、協方差等。

如果障礙概率為零,同樣的邏輯也適用不是通過 logit 模型進行參數化,而是通過任何其他二元回歸模型進行參數化,例如,在 1 處右刪失的計數分佈。當然,也可以是另一個計數分佈,例如負二項分佈。只有在零障礙和截斷計數部分之間存在共享參數時,整個分離才會失效。

一個突出的例子是,如果負二項分佈具有單獨的但很常見參數用於模型的兩個組件。(這可以在R-Forgehurdle(..., separate = FALSE, dist = "negbin", zero.dist = "negbin")countreg包中找到,它是實現的繼承者pscl。)

具體問題

**(a) 丟棄非常好的數據:**在你的情況下是的,一般來說不是。您有來自單個泊松模型的數據,沒有多餘的零(儘管有很多零)。因此,沒有必要為零和非零估計單獨的模型。但是,如果這兩個部分確實是由不同的參數驅動的,那麼就需要考慮這一點。

**(b) 由於大部分數據為零,可能會導致電源問題:**不一定。在這裡,您有三分之一的觀察結果是“成功”(跨欄)。這在二元回歸模型中不會被認為是非常極端的。(當然,如果不需要估計單獨的模型,您將獲得力量。)

**(c) 本身基本上不是一個“模型”,而只是順序運行兩個不同的模型:**這更具哲學性,我不會嘗試給出“一個”答案。相反,我會指出務實的觀點。對於模型估計,強調模型是獨立的會很方便,因為正如您所展示的那樣,您可能不需要專門的函數來進行估計。對於模型應用,例如預測或殘差等,將其視為單個模型會更方便。

**(d) 將這兩種情況都稱為“障礙”模型是否安全:**原則上是的。但是,行話可能因社區而異。例如,零障礙貝塔回歸更常見(並且非常容易混淆)稱為零膨脹貝塔回歸。就個人而言,我發現後者非常具有誤導性,因為 beta 分佈沒有可能被誇大的零點——但無論如何它是文獻中的標準術語。此外,tobit 模型是審查模型,因此不是障礙模型。但是,它可以通過一個概率(或 logit)模型加上一個截斷的正態模型來擴展。在計量經濟學文獻中,這被稱為 Cragg 兩部分模型。

軟件評論

位於https://R-Forge.R-project.org/R/?group_id=522countreg的R-Forge 上的包是 to / from的後續實現。它(仍然)不在 CRAN 上的主要原因是我們想要修改接口,可能以一種不完全向後兼容的方式。否則實現是相當穩定的。與之相比,它具有一些不錯的功能,例如:hurdle()``zeroinfl()``pscl``predict()``pscl

  • 使用與模型的零截斷部分zerotrunc()完全相同的代碼的函數。hurdle()因此,它提供了VGAM.
  • 此外,它作為 d/p/q/r 函數用於零截斷、跨欄和零膨脹計數分佈。這有助於將這些視為“一個”模型而不是單獨的模型。
  • 為了評估擬合優度,可以使用根圖和隨機分位數殘差圖等圖形顯示。(參見 Kleiber & Zeileis, 2016, The American Statistician , 70 (3), 296–303.doi:10.1080/00031305.2016.1173590。)

模擬數據

您的模擬數據來自單個泊松過程。如果e被視為已知的回歸量,那麼它將是標準的 Poisson GLM。如果e是未知噪聲分量,則存在一些未觀察到的異質性,導致一點點過度分散,這可以通過負二項式模型或其他類型的連續混合或隨機效應等捕獲。但是,由於e這裡的影響很小,這些都沒有太大的區別。下面,我將e其視為回歸量(即真實係數為 1),但您也可以省略它並使用負二項式或泊松模型。定性地說,所有這些都導致​​了類似的見解。

## 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

這反映了所有三個模型都可以一致地估計真實參數。查看相應的標準錯誤表明,在這種情況下(不需要障礙部分)泊松 GLM 更有效:

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

標准信息標準將選擇真正的 Poisson GLM 作為最佳模型:

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

Wald 檢驗可以正確檢測到障礙模型的兩個組成部分沒有顯著差異:

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

最後兩者都rootogram(p)表明qqrplot(p)Poisson GLM 非常適合數據,並且沒有過多的零或進一步錯誤指定的提示。

根圖+qqrplot

引用自:https://stats.stackexchange.com/questions/320924

comments powered by Disqus