Dlaczego reszty Pearsona z ujemnej regresji dwumianowej są mniejsze niż z regresji Poissona?

Mam te dane:

set.seed(1)
predictor  <- rnorm(20)
set.seed(1)
counts <- c(sample(1:1000, 20))
df <- data.frame(counts, predictor)

Przeprowadziłem regresję Poissona

poisson_counts <- glm(counts ~ predictor, data = df, family = "poisson")

I ujemna regresja dwumianowa:

require(MASS)
nb_counts <- glm.nb(counts ~ predictor, data = df)

Następnie obliczyłem statystyki dyspersji dla regresji Poissona:

sum(residuals(poisson_counts, type="pearson")^2)/df.residual(poisson_counts)

# [1] 145.4905

I ujemna regresja dwumianowa:

sum(residuals(nb_counts, type="pearson")^2)/df.residual(nb_counts)

# [1] 0.7650289

Czy ktokolwiek jest w stanie wyjaśnić, BEZ UŻYCIA RÓWNIKÓW, dlaczego statystyka dyspersji dla ujemnej regresji dwumianowej jest znacznie mniejsza niż statystyka dyspersji dla regresji Poissona?

— luciano
źródło

Odpowiedzi:

Jest to dość proste, ale „bez użycia równań” jest znaczną przeszkodą. Mogę to wyjaśnić słowami, ale te słowa koniecznie odzwierciedlają równania. Mam nadzieję, że będzie to dla ciebie akceptowalne / nadal miało jakąś wartość. (Odpowiednie równania nie są trudne.)

Istnieje kilka rodzajów resztek. Surowe reszty to po prostu różnica między zaobserwowanymi wartościami odpowiedzi (w twoim przypadku counts) a przewidywanymi wartościami odpowiedzi modelu. Reszty Pearsona dzielą je przez odchylenie standardowe (pierwiastek kwadratowy z funkcji wariancji dla konkretnej wersji uogólnionego modelu liniowego, którego używasz).

Odchylenie standardowe związane z rozkładem Poissona jest mniejsze niż odchylenie dwumianowe . Zatem, gdy dzielisz przez większy mianownik, iloraz jest mniejszy.

Ponadto ujemny dwumian jest bardziej odpowiedni dla twojego przypadku, ponieważ twój countsbędzie rozłożony jako jednolity w populacji. Oznacza to, że ich wariancja nie będzie równa ich średniej.

— gung - Przywróć Monikę
źródło

Chociaż PO prosi o wyjaśnienie niematematyczne, nadal miło byłoby zobaczyć matematyczne (lub równie rygorystyczne i jasne) uzasadnienie tej odpowiedzi. Po przeczytaniu pytania moja intuicja brzmiała: „Ponieważ Poisson jest (ograniczającym) specjalnym przypadkiem NB, a NB ma więcej parametrów, istnieje większa elastyczność dopasowania, więc oczywiście jakakolwiek rozsądna miara reszt nie powinna się zwiększać podczas wymiany Poisson GLM przez NB GLM. ” Zastanawiam się, czy taka intuicja była naprawdę poprawna.

— whuber

Jeśli , . Jeśli , i . Zatem wariancja Poissona jest równa średniej, wariancja NegBin jest większa niż średnia ( ). Dlatego „odchylenie standardowe związane z rozkładem Poissona jest mniejsze niż odchylenie dwumianowe”.

X \sim Poisson (λ)

$X\sim\text{Poisson}(\lambda)$

E [X] = V [X] = λ

$E[X]=V[X]=\lambda$

X \sim NegBin (r, p)

$X\sim\text{NegBin}(r,p)$

E [X] = p r / (1 - p)

$E[X]=pr/(1-p)$

V [X] = p r / (1 - p)^{2}

$V[X]=pr/(1-p)^2$

p < 1 \Rightarrow (1 - p)^{2} < (1 - p)

$p<1\Rightarrow (1-p)^2<(1-p)$

— Sergio

@Sergio Sedno sprawy polega na tym, że w modelu Poissona pracujemy z oszacowaniem zamiast z samym , aw modelu NB podobnie pracujemy z dwoma oszacowaniami i . . Twoje porównanie nie ma zatem bezpośredniego zastosowania. Bez faktycznego spisania wzorów dla MLE w obu modelach wcale nie jest oczywiste, jakie muszą być relacje między tymi zestawami oszacowań. Co więcej, reszta Pearsona jest współczynnikiem, a spór o wariancje odnosi się tylko do mianowników, co stanowi zaledwie połowę historii.

\hat{λ}

$\hat\lambda$

λ

$\lambda$

\hat{r}

$\hat{r}$

\hat{p}

$\hat{p}$

— whuber

Szacunki MLE są spójne. Problem polega na tym, że kiedy, jak mówi Gung, „liczba będzie rozłożona w populacji jako jednolita. Oznacza to, że ich wariancja nie będzie równa swojej średniej”, nigdy nie będziesz w stanie uzyskać oszacowanej wariancji Poissona większej niż szacowana Poissona, nawet jeśli twoje szacunki są obiektywne i spójne. To problem błędnej specyfikacji.

— Sergio

W przypadku modelu Poissona, jeśli oczekiwanie dla tej obserwacji wynosi jego wariancja wynosi , a zatem resztkowa wartość Pearson $i$ $Y_i$ $\mu_i$ $\mu_i$

\frac{y_{i} - {\hat{μ}}_{i}}{\sqrt{{\hat{μ}}_{i}}}

$\frac{y_i-\hat\mu_i}{\sqrt{\hat\mu_i}}$

gdzie jest oszacowaniem średniej. Wyjaśniono tutaj parametryzację ujemnego modelu dwumianowego stosowanego w MASS . Jeśli oczekiwaniem dla tej obserwacji jest jej wariancja wynosi , a zatem pozostałość Pearsona $\hat\mu$ $i$ $Y_i$ $\mu_i$ $\mu_i + \frac{\mu^2}{\theta}$

\frac{y_{i} - {\tilde{μ}}_{i}}{\sqrt{{\tilde{μ}}_{i} + \frac{{\tilde{μ}}^{' 2}}{θ}}}

$\frac{y_i-\tilde\mu_i}{\sqrt{\tilde\mu_i+\frac{\tilde\mu'^2}{\theta}}}$

gdzie jest oszacowaniem średniej. Im mniejsza wartość - tj. Bardziej dodatkowa wariancja Poissona - tym mniejsza wartość resztkowa w porównaniu do jej ekwiwalentu Poissona. [Ale jak zauważył @whuber, oszacowania średnich nie są takie same, , ponieważ procedura estymacji waży obserwacje zgodnie z ich zakładaną wariancją. Jeśli mielibyście wykonać pomiary dla tego wzoru predykcyjnego, zbliżyliby się i ogólnie dodanie parametru powinno dać lepsze dopasowanie do wszystkich obserwacji, chociaż nie wiem, jak to dokładnie wykazać. Mimo to szacowane ilości populacji są większe, jeśli model Poissona utrzymuje, więc nie powinno to być zaskoczeniem.] $\tilde\mu$ $\theta$ $\hat\mu\neq\tilde\mu$ $i$

— Scortchi - Przywróć Monikę
źródło

Dziękujemy za wprowadzenie niektórych równań. Ale czy w dwóch modelach będzie mieć takie same wartości? (Nie wydaje mi się.) Jeśli nie, to jak można porównać dwie pozostałości Pearsona?

μ_{i}

$\mu_i$

— whuber

@whuber W tym przypadku okazuje się, że dopasowane wartości dla obu modeli są prawie identyczne. W końcu „prawdziwy” model ma po prostu przecięcie i zasadniczo modeluje średnią, ponieważ w symulacji nie ma związku między x i Y.

— jsk

@jsk Tak, przejrzałem dane i uruchomiłem kod. (BTW, możliwa jest zmiana danych i uzyskanie zasadniczo tej samej statystyki dyspersji dla dwóch modeli). Niestety, wasz punkt, który jest prawidłowy, nadal nie rozwiązuje konkretnego pytania ani nie odnosi się do (domyślnego) pytania ogólnego na temat porównanie reszt Poissona z resztami NB, ponieważ oszacowane wariancje również mogą być prawie identyczne. Jednym potencjalnie mylącym aspektem obecnej odpowiedzi jest użycie symbolu „ ” w odniesieniu do tego, co (w zasadzie) mogłoby być różnymi szacunkami w dwóch modelach tych samych danych.

μ_{i}

$\mu_i$

— whuber

@ whuber Rzeczywiście, masz ważne punkty na temat korzystania z . Co ciekawe, nie mogę znaleźć sposobu na symulację danych, które spowodowałyby niższą statystykę dyspersji dla Poissona niż NB. Być może nie jest to możliwe? Zgadzam się, że ma to sens intuicyjnie. Nie jest to łatwe do udowodnienia, ponieważ nie ma zamkniętego rozwiązania dla mle, gdy masz glm z funkcją link inną niż tożsamość. Ale tak, łatwo jest uczynić te dwie statystyki dyspersji bardzo podobnymi.

μ_{i}

$\mu_i$

— jsk

@jsk - jednym teoretycznym argumentem, aby podejrzewać, że model NB zawsze będzie pasował lepiej niż Poisson, jest to, że można napisać NB jako rozkład związku poissona-gamma. Masz więc a następnie daje ujemny model dwumianowy . Teraz dodanie tych parametrów pozwala modelowi zbliżyć przewidywaną średnią do obserwowanej wartości (gdy zobaczysz , zmniejszając resztkową).

(y_{i} | λ, v_{i}, r) \sim P o i s s o n (λ v_{i})

$( y_i|\lambda, v_i, r)\sim Poisson (\lambda v_i)$

(v_{i} | λ, r) \sim G a m m a (r, r)

$(v_i|\lambda, r)\sim Gamma(r, r)$

(y_{i} | λ, r) \sim N B (r, \frac{λ}{r + λ})

$( y_i|\lambda, r)\sim NB (r,\frac {\lambda}{r+\lambda} )$

v_{i}

$v_i$

y_{i} > λ

$y_i>\lambda$

v_{i} > 1

$v_i> 1$

— prawdopodobieństwo