Ist ein „Hürdenmodell“ wirklich ein Modell? Oder nur zwei separate, sequentielle Modelle?

Betrachten Sie ein Hürdenmodell, das die Zähldaten yvon einem normalen Prädiktor vorhersagt 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 

In diesem Fall habe ich Zähldaten mit 69 Nullen und 31 positiven Zählwerten. Denken Sie im Moment nicht daran, dass dies per Definition des Datenerzeugungsverfahrens ein Poisson-Prozess ist, da es sich bei meiner Frage um Hürdenmodelle handelt.

Nehmen wir an, ich möchte diese überschüssigen Nullen mit einem Hürdenmodell behandeln. Aus meiner Lektüre über sie ging hervor, dass Hürdenmodelle an sich keine tatsächlichen Modelle sind - sie führen nur zwei verschiedene Analysen nacheinander durch. Erstens eine logistische Regression, die vorhersagt, ob der Wert positiv gegenüber Null ist oder nicht. Zweitens eine nullverengte Poisson-Regression, bei der nur die Nicht-Null-Fälle berücksichtigt werden. Dieser zweite Schritt fühlte sich für mich falsch an, weil er (a) vollkommen gute Daten wegwirft, was (b) zu Stromproblemen führen kann, da viele der Daten Nullen sind, und (c) im Grunde genommen kein "Modell" an und für sich ist , aber nur nacheinander laufen zwei verschiedene Modelle.

Also habe ich ein "Hürdenmodell" ausprobiert, anstatt nur die logistische und die nullverengte Poisson-Regression getrennt auszuführen. Sie gaben mir identische Antworten (ich kürze die Ausgabe der Kürze halber ab):

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

Dies scheint mir unverständlich zu sein, da viele verschiedene mathematische Darstellungen des Modells die Wahrscheinlichkeit beinhalten, dass eine Beobachtung bei der Schätzung von positiven Zählfällen ungleich Null ist, aber die Modelle, die ich oben ausgeführt habe, ignorieren sich vollständig. Dies ist zum Beispiel aus Kapitel 5, Seite 128 von Smithson & Merkles verallgemeinerten linearen Modellen für kategoriale und stetig begrenzte abhängige Variablen :

... Zweitens muss die Wahrscheinlichkeit, dass irgendeinen Wert annimmt (null und die positiven ganzen Zahlen), gleich eins sein. Dies ist in Gleichung (5.33) nicht garantiert. Um dieses Problem zu lösen, multiplizieren wir die Poisson-Wahrscheinlichkeit mit der Bernoulli-Erfolgswahrscheinlichkeit .      Diese Probleme erfordern, dass wir das obige Hürdenmodell als ausdrücken wobei , ,$y$P ( Y = y | x , z , β , γ ) = { 1 - π für  y = 0 π$\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$sind die Kovariaten für das Poisson-Modell, sind die Kovariaten für das logistische Regressionsmodell, und und sind die jeweiligen Regressionskoeffizienten ... . $\boldsymbol z$$\hat{\boldsymbol{\beta}}$$\hat{\boldsymbol{\gamma}}$

Indem ich die beiden Modelle vollständig voneinander trenne - was Hürdenmodelle zu sein scheinen -, sehe ich nicht, wie in die Vorhersage von Fällen mit positiver Zählung einbezogen wird. Aufgrund der Tatsache, dass ich die Funktion nur mit zwei verschiedenen Modellen replizieren konnte , kann ich nicht erkennen, wie eine Rolle im verkürzten Poisson spielt Regression überhaupt.$\hat{\pi}$hurdle$\text{logit}^{-1}(z\hat{\gamma})$

Verstehe ich Hürdenmodelle richtig? Es scheint, als würden zwei nur zwei aufeinanderfolgende Modelle ausführen: Erstens eine Logistik; Zweitens ein Poisson, bei dem Fälle mit völlig ignoriert werden . Ich würde mich freuen, wenn jemand meine Verwechslung mit dem .$y = 0$$\hat{\pi}$

Wenn ich richtig liege, was sind das für Hürdenmodelle, was ist die Definition eines "Hürdenmodells" im Allgemeinen? Stellen Sie sich zwei verschiedene Szenarien vor:

• Stellen Sie sich vor, Sie modellieren die Wettbewerbsfähigkeit von Wahlkämpfen anhand der Wettbewerbsfähigkeitswerte (1 - (Stimmenanteil der Gewinner - Stimmenanteil der Zweitplatzierten)). Dies ist [0, 1), weil es keine Bindungen gibt (zB 1). Ein Hürdenmodell ist hier sinnvoll, weil es einen Prozess gibt (a) War die Wahl unbestritten? und (b) wenn nicht, welche vorausgesagte Wettbewerbsfähigkeit? Also führen wir zuerst eine logistische Regression durch, um 0 vs. (0, 1) zu analysieren. Dann führen wir eine Beta-Regression durch, um die (0, 1) -Fälle zu analysieren.

• Stellen Sie sich eine typische psychologische Studie vor. Die Antworten lauten [1, 7] wie bei einer herkömmlichen Likert-Skala mit einem riesigen Deckeneffekt von 7. Man könnte ein Hürdenmodell erstellen, das eine logistische Regression von [1, 7) gegenüber 7 und dann eine Tobit-Regression für alle Fälle, in denen beobachtete Antworten sind <7.

Wäre es sicher, beide Situationen als "Hürden" -Modelle zu bezeichnen , selbst wenn ich sie mit zwei sequentiellen Modellen einschätze (logistisch und dann Beta im ersten Fall, logistisch und dann Tobit im zweiten)?

Trennen der Log-Wahrscheinlichkeit

$\hat \pi$

$\ell_1(\gamma)$$\ell_2(\beta)$

$\pi$$f(\cdot)$

$\mu$$\theta$hurdle(..., separate = FALSE, dist = "negbin", zero.dist = "negbin")countregpscl

Konkrete Fragen

(a) Perfekt gute Daten wegwerfen: In Ihrem Fall ja, im Allgemeinen nein. Sie haben Daten aus einem einzelnen Poisson-Modell ohne überzählige Nullen (wenn auch viele Nullen ). Daher ist es nicht erforderlich, separate Modelle für die Nullen und Nichtnullen zu schätzen. Wenn die beiden Teile jedoch wirklich von unterschiedlichen Parametern gesteuert werden, muss dies berücksichtigt werden.

(b) Kann zu Stromversorgungsproblemen führen, da viele Daten Nullen sind: Nicht unbedingt. Hier haben Sie ein Drittel der Beobachtungen, die "Erfolge" (Hürdenkreuzungen) sind. Dies würde in einem binären Regressionsmodell nicht als extrem angesehen. (Wenn es nicht notwendig ist, separate Modelle zu schätzen, würden Sie natürlich an Leistung gewinnen.)

(c) Im Grunde genommen kein "Modell" an sich, sondern nur zwei verschiedene Modelle nacheinander: Dies ist philosophischer und ich werde nicht versuchen, "eine" Antwort zu geben. Stattdessen werde ich auf pragmatische Gesichtspunkte hinweisen. Für das Modellschätzung kann es bequem zu betonen, dass die Modelle getrennt sind , weil - wie Sie zeigen - Sie nicht eine spezielle Funktion für die Schätzung benötigen. Für das Modell Anwendung , zum Beispiel für Vorhersagen oder Rückstände usw., kann es günstiger sein , dies als ein einziges Modell zu sehen.

(d) Wäre es sicher, diese beiden Situationen als Hürdenmodelle zu bezeichnen: Im Prinzip ja. Die Fachsprache kann jedoch in den verschiedenen Communities unterschiedlich sein. Beispielsweise wird die Beta-Regression ohne Hürde häufiger (und sehr verwirrend) als Beta-Regression ohne Inflation bezeichnet. Persönlich finde ich letzteres sehr irreführend, da die Beta-Verteilung keine Nullen hat, die aufgeblasen werden könnten - aber es ist der Standardbegriff in der Literatur trotzdem. Darüber hinaus ist das Tobit-Modell ein zensiertes Modell und somit kein Hürdenmodell . Es könnte jedoch um ein Probit- (oder Logit-) Modell und ein abgeschnittenes normales Modell erweitert werden. In der ökonometrischen Literatur ist dies als zweiteiliges Cragg-Modell bekannt.

Software-Kommentare

Das countregPaket auf R-Forge bei https://R-Forge.R-project.org/R/?group_id=522 ist der Nachfolger Umsetzung zu hurdle()/ zeroinfl()ab pscl. Der Hauptgrund dafür, dass es (noch) nicht bei CRAN ist, ist, dass wir die predict()Schnittstelle möglicherweise auf eine Weise überarbeiten möchten, die nicht vollständig abwärtskompatibel ist. Ansonsten ist die Implementierung ziemlich stabil. Im Vergleich dazu gibt pscles ein paar nette Features, zB:

• Eine zerotrunc()Funktion, die genau den gleichen Code verwendet wie hurdle()für den Teil des Modells, bei dem keine Kürzungen vorgenommen wurden. Somit bietet es eine Alternative zu VGAM.

• Darüber hinaus fungiert es als d / p / q / r für die Zählerverteilungen mit Null-Kürzung, Hürde und Null-Inflation. Dies erleichtert das Betrachten dieser Modelle als "ein" Modell und nicht als separate Modelle.

• Zur Beurteilung der Anpassungsgüte stehen grafische Darstellungen wie Rootogramme und randomisierte Quantilrestkurven zur Verfügung. (Siehe Kleiber & Zeileis, 2016, The American Statistician , 70 (3), 296–303. Doi: 10.1080 / 00031305.2016.1173590 .)

Simulierte Daten

Ihre simulierten Daten stammen aus einem einzelnen Poisson-Prozess. Wenn dies eals bekannter Regressor behandelt wird, handelt es sich um einen Standard-Poisson-GLM. Wenn ees sich um eine unbekannte Rauschkomponente handelt, gibt es eine unbeobachtete Heterogenität, die ein wenig Überdispersion verursacht und durch ein negatives Binomialmodell oder eine andere Art von kontinuierlichem Gemisch oder zufälligen Effekt usw. erfasst werden kann. Der Effekt von eist hier jedoch eher gering , nichts davon macht einen großen Unterschied. Im Folgenden werde ich eals Regressor behandelt (dh mit dem wahren Koeffizienten 1), aber Sie können dies auch weglassen und negative Binomial- oder Poisson-Modelle verwenden. Qualitativ führen all diese Erkenntnisse zu ähnlichen Erkenntnissen.

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

Dies zeigt, dass alle drei Modelle die wahren Parameter konsistent schätzen können. Ein Blick auf die entsprechenden Standardfehler zeigt, dass das Poisson-GLM in diesem Szenario (ohne die Notwendigkeit eines Hürdenabschnitts) effizienter ist:

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

Standardinformationskriterien würden das wahre Poisson GLM als bestes Modell auswählen:

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

Und ein Wald-Test würde richtig erkennen, dass sich die beiden Komponenten des Hürdenmodells nicht signifikant unterscheiden:

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

Schließlich beides rootogram(p)und qqrplot(p)zeigen, dass der Poisson GLM sehr gut zu den Daten passt und dass es keine überzähligen Nullen oder Hinweise auf weitere Fehlspezifikationen gibt.

Ich stimme zu, dass der Unterschied zwischen Null-Inflations- und Hürden-Modellen schwer zu verstehen ist. Beide sind eine Art Mischmodell. Nach allem, was ich sagen kann, besteht der wichtige Unterschied darin, dass Sie in einem Modell mit Null-Inflation eine Masse bei Null mit einer Verteilung \ textit {mischen, die auch den Wert Null annehmen kann. Für ein Hürdenmodell mischen Sie eine Masse bei Null mit einer Verteilung, die nur Werte größer als 0 annimmt. Daher können Sie im Modell mit Null-Inflation zwischen "strukturellen Nullen" (die der Masse bei Null entsprechen) und "Stichproben-Nullen" unterscheiden Dies entspricht dem zufälligen Auftreten von 0 in dem Modell, in das Sie sich einmischen. Natürlich hängt diese Identifikation stark von der Auswahl der richtigen Verteilung ab! Aber wenn Sie zum Beispiel einen Poisson ohne Luftdruck haben, Sie können zwischen Nullen, die von der Poisson-Komponente stammen (Stichproben-Nullen) und Nullen, die von der Masse bei Null stammen (strukturelle Nullen), unterscheiden. Wenn Sie ein Null-Inflationsmodell haben und die Verteilung, die Sie einmischen, keine Masse bei Null hat, kann dies als Hürdenmodell interpretiert werden.

In Bezug auf den philosophischen Aspekt, "wann sollten wir ein einzelnes Modell betrachten und wann zwei separate Modelle" , kann es interessant sein zu bemerken, dass die Stichprobenschätzungen der Modellparameter korreliert sind.

In der folgenden Darstellung mit einer Simulation sehen Sie hauptsächlich die Korrelation zwischen der Steigung und dem Achsenabschnitt des Zählteils. Es gibt aber auch eine leichte Beziehung zwischen dem Zählteil und dem Hürdenanteil. Wenn Sie die Parameter ändern, z. B. das Lambda in der Poisson-Verteilung oder die Stichprobengröße verkleinern, wird die Korrelation stärker.

Daher würde ich sagen, dass Sie es nicht als zwei separate Modelle betrachten sollten. Oder zumindest besteht eine gewisse Beziehung, selbst wenn Sie in der Praxis die beiden Schätzungen unabhängig voneinander berechnen können.

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]))

Es macht nicht viel Sinn, dass es eine Korrelation zwischen den beiden Teilen gibt. Möglicherweise liegt dies jedoch daran, dass die Schätzwerte für die Parameter im Poisson-Modell diskret sind und in welchem ​​Verhältnis sie zur Anzahl der Nullen stehen.