Rozważ model przeszkodowy przewidujący zliczanie danych y
z 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 , ,
są zmiennymi towarzyszącymi modelu Poissona, są zmiennymi logistycznymi modelu regresji, a i są odpowiednimi współczynnikami regresji ... .
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.hurdle
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 .
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)?
źródło
pscl::hurdle
, ale wygląda to tak samo w równaniu 5 tutaj: cran.r-project.org/web/packages/pscl/vignettes/countreg.pdf A może ja wciąż brakuje mi czegoś podstawowego, co sprawiłoby, że kliknęło by mnie?hurdle()
. W naszej sparowanej / winiecie staramy się jednak podkreślić bardziej ogólne elementy składowe.Odpowiedzi:
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.π^
W zapisie Smithsona i Merkle:
The same logic also works if the zero hurdle probabilityπ 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(⋅) 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μ but common θ parameters are employed in the two components of the model. (This is available in
hurdle(..., separate = FALSE, dist = "negbin", zero.dist = "negbin")
in thecountreg
package from R-Forge, the successor to thepscl
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
countreg
Pakiet R-Forge w https://R-Forge.R-project.org/R/?group_id=522 jest następcą wdrożenie dohurdle()
/zeroinfl()
zpscl
. 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 zpscl
tym ma kilka fajnych funkcji, np .:zerotrunc()
Funkcja, która używa dokładnie tego samego kodu, jakhurdle()
dla zerowego ściętej części modelu. W ten sposób oferuje alternatywę dlaVGAM
.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. Ife
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 ofe
is rather small here, none of this makes a big difference. Below, I'm treatinge
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.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:
Standardowe kryteria informacyjne wybrałyby prawdziwy Poisson GLM jako najlepszy model:
Test Walda prawidłowo wykryłby, że dwa elementy modelu przeszkody nie różnią się znacząco:
Wreszcie zarówno
rootogram(p)
iqqrplot(p)
pokazują, że Poissona GLM wpisuje dane i bardzo dobrze, że nie ma nadmiaru zer lub wskazówki dotyczące dalszych misspecifications.źródło
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.
źródło
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.
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.
źródło
truncdist::rtrunc(n = 100, spec = 'pois', a = 0, b = Inf, lambda = exp(-1.5 + rnorm(100)))
yields an error (using version 1.0.2):Error in if (G.a == G.b) { : the condition has length > 1
. In any case, usingrhpois()
from packagecountreg
on R-Forge is easier for simulating from a hurdle Poisson model with given hurdle crossing probabilitypi
and underlying (untruncated) Poisson expectationlambda
. If I use these I get only very small empirical correlations between the zero hurdle and truncated count parts.dgp <- function(n = 100, b = c(-0.5, 2), g = c(0.5, -2)) { x <- runif(n, -1, 1) ; y <- rhpois(n, lambda = exp(b[1] + b[2] * x), pi = plogis(g[1] + g[2] * x)); data.frame(x = x, y = y) }
Simulation:set.seed(1); cf <- t(replicate(3000, coef(hurdle(y ~ x, data = dgp()))))
. Evaluation:pairs(cf)
andcor(cf)
. CheckingcolMeans(cf)
also shows that the estimation has worked reasonably well.