Czy „model przeszkodowy” to naprawdę jeden model? Czy tylko dwa oddzielne, sekwencyjne modele?

Rozważ model przeszkodowy przewidujący zliczanie danych yz normalnego predyktora 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 

W tym przypadku mam dane z 69 zerami i 31 dodatnimi. Nie ważne, że jest to z definicji procedura generowania danych proces Poissona, ponieważ moje pytanie dotyczy modeli przeszkód.

Powiedzmy, że chcę poradzić sobie z tymi zerami za pomocą modelu przeszkód. Z mojej lektury na ich temat wydawało się, że modele przeszkód same w sobie nie są rzeczywistymi modelami - po prostu wykonują dwie różne analizy sekwencyjnie. Po pierwsze, regresja logistyczna przewidująca, czy wartość jest dodatnia w stosunku do zera. Po drugie, zerowana regresja Poissona, obejmująca tylko przypadki niezerowe. Ten drugi krok wydawał mi się niewłaściwy, ponieważ (a) wyrzucamy idealnie dobre dane, co (b) może prowadzić do problemów z zasilaniem, ponieważ większość danych to zera, i (c) zasadniczo nie jest „modelem” samym w sobie , ale po prostu uruchamiają kolejno dwa różne modele.

Wypróbowałem więc „model przeszkody” zamiast po prostu osobno uruchomić logistyczną i zerowaną regresję Poissona. Dali mi identyczne odpowiedzi (skrótowo, dla uproszczenia):

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

Wydaje mi się to nie na miejscu, ponieważ wiele różnych matematycznych reprezentacji modelu obejmuje prawdopodobieństwo, że obserwacja jest różna od zera w szacowaniu przypadków zliczania dodatniego, ale modele, które przedstawiłem powyżej, całkowicie się ignorują. Na przykład pochodzi z Rozdziału 5, strona 128 Uogólnionych modeli liniowych Smithsona i Merkle'a dla zmiennych jakościowych i ciągłych ograniczonych zależnych :

... Po drugie, prawdopodobieństwo, że przyjmie dowolną wartość (zero i dodatnie liczby całkowite), musi wynosić jeden. Nie jest to zagwarantowane w równaniu (5.33). Aby poradzić sobie z tym problemem, mnożymy prawdopodobieństwo Poissona przez prawdopodobieństwo sukcesu Bernoulliego π .      Te problemy wymagają od nas wyrażenia powyższego modelu przeszkód jako gdzie , ,$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$są zmiennymi towarzyszącymi modelu Poissona, są zmiennymi logistycznymi modelu regresji, a i są odpowiednimi współczynnikami regresji ... . $\boldsymbol z$$\hat{\boldsymbol{\beta}}$$\hat{\boldsymbol{\gamma}}$

Wykonując dwa modele całkowicie od siebie oddzielone - co wydaje się być tym, co robią modele przeszkodowe - nie widzę, jak jest włączony do przewidywania przypadków zliczania dodatniego. Ale w oparciu o to, jak mogłem zreplikować funkcję, uruchamiając tylko dwa różne modele, nie widzę, jak odgrywa rolę w okrojonym Poissonie regresja w ogóle.$\hat{\pi}$hurdle$\text{logit}^{-1}(z\hat{\gamma})$

Czy rozumiem poprawnie modele przeszkód? Wydaje się, że dwa prowadzą tylko dwa modele sekwencyjne: Po pierwsze, logistyka; Po drugie, Poissona, całkowicie ignorując przypadki, w których . Byłbym wdzięczny, gdyby ktoś mógł wyjaśnić moje zamieszanie związane z biznesem .$y = 0$$\hat{\pi}$

---

Jeśli mam rację, że takie są modele przeszkód, to jaka jest definicja modelu „przeszkód”, bardziej ogólnie? Wyobraź sobie dwa różne scenariusze:

• Wyobraź sobie modelowanie konkurencyjności ras wyborczych, patrząc na wyniki konkurencji (1 - (odsetek głosów zwycięzcy - odsetek głosów drugich). Jest to [0, 1), ponieważ nie ma żadnych powiązań (np. 1). Model przeszkód ma tutaj sens, ponieważ istnieje jeden proces (a) czy wybory były bezsporne? oraz (b) jeśli nie, co przewidywało konkurencyjność? Najpierw wykonujemy regresję logistyczną, aby przeanalizować 0 vs. (0, 1). Następnie wykonujemy regresję beta, aby przeanalizować przypadki (0, 1).

• Wyobraź sobie typowe badanie psychologiczne. Odpowiedzi wynoszą [1, 7], podobnie jak tradycyjna skala Likerta, z ogromnym efektem pułapu na poziomie 7. Można zrobić model przeszkody, w którym regresja logistyczna wynosi [1, 7) vs. 7, a następnie regresja Tobit dla wszystkich przypadków, w których zaobserwowane odpowiedzi wynoszą <7.

Czy bezpiecznie byłoby nazwać obie te modele „przeszkodą” , nawet jeśli oszacuję je za pomocą dwóch modeli sekwencyjnych (w pierwszym przypadku logistyka, a następnie beta, w drugim logistyka, a potem Tobit)?

Rozdzielenie prawdopodobieństwa dziennika

Prawdą jest, że większość modeli przeszkód można oszacować osobno (powiedziałbym, zamiast sekwencyjnie ). Powodem jest to, że prawdopodobieństwo dziennika można rozłożyć na dwie części, które można osobno zmaksymalizować. To dlatego, że π jest tylko współczynnik skalowania w (5,34), który staje się dodatek termin w dzienniku prawdopodobieństwa.$\hat \pi$

W zapisie Smithsona i Merkle:

$$\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)$

$\ell_1(\gamma)$ (binary logit model) and $\ell_2(\beta)$ (zero-truncated Poisson model) can be maximized separately, leading to the same parameter estimates, covariances, etc. as in the case where they are maximized jointly.

The same logic also works if the zero hurdle probability $\pi$ is not parametrized through a logit model but any other binary regression model, e.g., a count distribution right-censored at 1. And, of course, $f(\cdot)$ could also be another count distribution, e.g., negative binomial. The whole separation only breaks down if there are shared parameters between the zero hurdle and the truncated count part.

A prominent example would be if negative binomial distributions with separate $\mu$ but common $\theta$ parameters are employed in the two components of the model. (This is available in hurdle(..., separate = FALSE, dist = "negbin", zero.dist = "negbin") in the countreg package from R-Forge, the successor to the pscl implementation.)

Concrete questions

(a) Wyrzucanie idealnie dobrych danych: w twoim przypadku tak, ogólnie nie. Masz dane z jednego modelu Poissona bez zer zerowych (choć wiele zer ). Dlatego nie jest konieczne szacowanie osobnych modeli dla zer i zer. Jednak jeśli dwie części są naprawdę napędzane różnymi parametrami, należy to uwzględnić.

(b) Może to prowadzić do problemów z zasilaniem, ponieważ większość danych to zera: niekoniecznie. Tutaj masz jedną trzecią obserwacji, które są „sukcesami” (przejazdy przez przeszkody). Nie byłoby to uważane za ekstremalne w modelu regresji binarnej. (Oczywiście, jeśli oszacowanie osobnych modeli nie jest konieczne, zyskałbyś moc.)

(c) Basically not a 'model' in and of itself, but just sequentially running two different models: This is more philosophical and I won't try to give "one" answer. Instead, I will point out pragmatic points of view. For model estimation, it can be convenient to emphasize that the models are separate because - as you show - you might not need a dedicated function for the estimation. For model application, e.g., for predictions or residuals etc., it can be more convenient to see this as a single model.

(d) Would it be safe to call both of these situations 'hurdle' models: In principle yes. However, jargon may vary across communities. For example, the zero-hurdle beta regression is more commonly (and very confusingly) called zero-inflated beta regression. Personally, I find the latter very misleading because the beta distribution has no zeros that could be inflated - but it's the standard term in the literature anyway. Moreover, the tobit model is a censored model and hence not a hurdle model. It could be extended, though, by a probit (or logit) model plus a truncated normal model. In the econometrics literature this is known as the Cragg two-part model.

Software comments

countregPakiet R-Forge w https://R-Forge.R-project.org/R/?group_id=522 jest następcą wdrożenie do hurdle()/ zeroinfl()z pscl. Głównym powodem, dla którego nie jest (nadal) w CRAN, jest to, że chcemy zmienić predict()interfejs, być może w sposób, który nie jest w pełni kompatybilny wstecz. W przeciwnym razie implementacja jest dość stabilna. W porównaniu z pscltym ma kilka fajnych funkcji, np .:

• zerotrunc()Funkcja, która używa dokładnie tego samego kodu, jak hurdle()dla zerowego ściętej części modelu. W ten sposób oferuje alternatywę dla VGAM.

• Co więcej, działa jako d / p / q / r dla rozkładów zliczania zerowanego, przeszkodowego i zerowanego. Ułatwia to postrzeganie ich jako „jednego” modelu, a nie oddzielnych modeli.

• For assessing the goodness of fit, graphical displays like rootograms and randomized quantile residual plots are available. (See Kleiber & Zeileis, 2016, The American Statistician, 70(3), 296–303. doi:10.1080/00031305.2016.1173590.)

Simulated data

Your simulated data comes from a single Poisson process. If e is treated as a known regressor then it would be a standard Poisson GLM. If e is an unknown noise component, then there is some unobserved heterogeneity causing a little bit of overdispersion which could be captured by a negative binomial model or some other kind of continuous mixture or random effect etc. However, as the effect of e is rather small here, none of this makes a big difference. Below, I'm treating e as a regressor (i.e., with true coefficient of 1) but you could also omit this and use negative binomial or Poisson models. Qualitatively, all of these lead to similar insights.

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

Odzwierciedla to, że wszystkie trzy modele mogą konsekwentnie oszacować prawdziwe parametry. Patrząc na odpowiednie błędy standardowe pokazuje, że w tym scenariuszu (bez potrzeby stosowania części przeszkody) Poisson GLM jest bardziej wydajny:

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

Standardowe kryteria informacyjne wybrałyby prawdziwy Poisson GLM jako najlepszy model:

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

Test Walda prawidłowo wykryłby, że dwa elementy modelu przeszkody nie różnią się znacząco:

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

Wreszcie zarówno rootogram(p)i qqrplot(p)pokazują, że Poissona GLM wpisuje dane i bardzo dobrze, że nie ma nadmiaru zer lub wskazówki dotyczące dalszych misspecifications.

I agree the difference between zero-inflated and hurdle models is hard to understand. Both are a kind of mixture model. From what I can tell, the important difference is, in a zero-inflated model, you mix a mass at zero with a distribution \textit{that can also take the value zero}. For a hurdle model, you mix a mass at zero with a distribution that only takes values greater than 0. Thus, in the zero-inflated model you can distinguish between 'structural zeros' (corresponding to the mass at zero) and 'sampling zeros' corresponding to the chance occurrence of 0 from the model you are mixing in. Of course, this identification depends strongly on making the right choice of distribution! But, if you have a zero-inflated Poisson, for example, you can distinguish between zeros that come from the Poisson component (sampling zeros) and zeros that come from the mass at zero (structural zeros). If you have a zero-inflated model and the distribution you are mixing in has no mass at zero, it could be interpreted as a hurdle model.

Regarding the philosophical aspect, "when should we consider something a single model and when two separate models", it may be interesting to note that the sample estimates of the model-parameters are correlated.

In the plot below with a simulation you mostly see the correlation between the slope and the intercept of the counts part. But there is also some slight relation between the counts part and the hurdle part. If you change the parameters, e.g. make the lambda in the Poisson distribution smaller or the sample size smaller, then the correlation becomes stronger.

So I would say that you should not consider it as two separate models. Or at least there is some relation even when in practice you can compute the two estimates independent from each other.

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

---

It doesn't make much sense that there is a correlation between the two parts. But it might probably be due to discrete levels of the estimates for the parameters in the Poisson model, and how these relate to the number of zero's.