Dlaczego quasi-Poissona w GLM nie traktuje się jako specjalnego przypadku ujemnego dwumianu?

21

Próbuję dopasować uogólnione modele liniowe do niektórych zestawów danych zliczania, które mogą być rozproszone lub nie. Dwa obowiązujące tutaj rozkłady kanoniczne to Poisson i ujemny dwumianowy (Negbin), z EV i wariancjąμ

VarP=μ

VarNB=μ+μ2θ

który może być wyposażony w R z zastosowaniem glm(..,family=poisson)a glm.nb(...), odpowiednio. Jest też quasipoissonrodzina, która w moim rozumieniu jest dostosowanym Poissonem o tym samym EV i wariancji

VarQP=ϕμ ,

tj. wpadnięcie gdzieś pomiędzy Poissonem i Negbin. Główny problem z rodziną quasipoisson polega na tym, że nie ma dla niego odpowiedniego prawdopodobieństwa, a zatem wiele niezwykle przydatnych testów statystycznych i miar dopasowania (AIC, LR itp.) Jest niedostępnych.

Jeśli porównasz wariancje QP i Negbin, możesz zauważyć, że możesz je zrównać, umieszczając . Kontynuując tę ​​logikę, możesz spróbować wyrazić rozkład quasipoisson jako szczególny przypadek Negbina:ϕ=1+μθ

QP(μ,ϕ)=NB(μ,θ=μϕ1) ,

tzn. Negbin z liniowo zależnym od . Próbowałem zweryfikować ten pomysł, generując losową sekwencję liczb zgodnie z powyższym wzorem i dopasowując go do :μθμglm

#fix parameters

phi = 3
a = 1/50
b = 3
x = 1:100

#generating points according to an exp-linear curve
#this way the default log-link recovers the same parameters for comparison

mu = exp(a*x+b) 
y = rnbinom(n = length(mu), mu = mu, size = mu/(phi-1)) #random negbin generator

#fit a generalized linear model y = f(x)  
glmQP = glm(y~x, family=quasipoisson) #quasipoisson
glmNB = glm.nb(y~x) #negative binomial

> glmQP

Call:  glm(formula = y ~ x, family = quasipoisson)

Coefficients:
(Intercept)            x  
    3.11257      0.01854  
(Dispersion parameter for quasipoisson family taken to be 3.613573)

Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
Null Deviance:      2097 
Residual Deviance: 356.8    AIC: NA

> glmNB

Call:  glm.nb(formula = y ~ x, init.theta = 23.36389741, link = log)

Coefficients:
(Intercept)            x  
    3.10182      0.01873  

Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
Null Deviance:      578.1 
Residual Deviance: 107.8    AIC: 824.7

Oba pasowania odtwarzają parametry, a quasipoisson daje „rozsądne” oszacowanie dla . Możemy teraz również zdefiniować wartość AIC dla quasipoissonu:ϕ

df = 3 # three model parameters: a,b, and phi
phi.fit = 3.613573 #fitted phi value copied from summary(glmQP)
mu.fit = glmQP$fitted.values 

#dnbinom = negbin density, log=T returns log probabilities
AIC = 2*df - 2*sum(dnbinom(y, mu=mu.fit, size = mu.fit/(phi.fit - 1), log=T))
> AIC
[1] 819.329

(Musiałem ręcznie skopiować dopasowaną wartość , ponieważ nie mogłem jej znaleźć w obiekcie)ϕsummary(glmQP)glmQP

Ponieważ oznaczałoby to, że quasipoisson jest, co nie dziwi, lepszym dopasowaniem; więc co najmniej robi to, co powinien, i dlatego może być rozsądną definicją AIC (i, co za tym idzie, prawdopodobieństwa) quasipoissonu. Powstają więc wielkie pytania A I C Q PAICQP<AICNBAICQP

  1. Czy ten pomysł ma sens? Czy moja weryfikacja opiera się na okrągłym uzasadnieniu?
  2. Główne pytanie dla każdego, kto „wymyśla” coś, co wydaje się brakować w dobrze ugruntowanym temacie: jeśli ten pomysł ma sens, dlaczego nie został jeszcze wdrożony glm?

Edycja: rysunek dodany

glm pasuje i pasma sigma + -1

użytkownik 28400
źródło
1
(+1) Witamy w Cross Validated! I dziękuję za doskonałe pytanie (choć kilka komentarzy w kodzie może być miłych dla osób, które nie używają R). Myślę, że mogłeś wymyślić na nowo model NB1 (choć nie śledziłem go jeszcze szczegółowo). Zauważ też, że nie istnieje rozkład quasi-Poissona - dlatego nie ma prawdopodobieństwa ani AIC - odnosi się tylko do sposobu dopasowania środków i wariancji.
Scortchi - Przywróć Monikę
2
Dzięki! W międzyczasie dodałem kilka komentarzy, mam nadzieję, że to wyjaśni. Rozumiem, że rozkład quasi-Poissona nie istnieje per se - co ja naprawdę próbuje dowiedzieć się, dlaczego QP jest jeszcze coś w ogóle, biorąc pod uwagę, że rozkład NB1 istnieje i nie ma żadnej z jednostek typu problemów QP za (patrz odpowiedź Achimsa na pozorne rozwiązanie).
user28400,
1
@Scortchi --- rzeczywiście nie jest taki rozkład ... Jeżeli , a , a jest wykładniczy rodziny o średniej i wariancji . Jeśli . To niekoniecznie nadaje się do zliczania danych (z wyjątkiem jako przybliżenie), ponieważ jest zdefiniowana na . Y = K X Y μ = k λ k μ k 1 0 , k , 2 k , . . .XPois(λ)Y=kXYμ=kλkμk10,k,2k,...
Glen_b
1
@Glen_b: Czy ludzie naprawdę nazywają to quasi-Poisson? W każdym razie jest to ładna ilustracja - gdy używasz modelu „quasiPoisson”, tak naprawdę nie zakładasz, że ten rozkład, NB1 lub jakikolwiek inny, to tylko związek między średnią i wariancją, który sprawia, że ​​twoje oszacowania współczynników i ich standardowych błędów lepiej, gdy próbka staje się większa.
Scortchi - Przywróć Monikę
1
@Scortchi Jest to jedyny wykładniczy rozkład rodziny, który spełnia założenia quasi-Poissona, więc coś w rodzaju - czasami widziałem, jak ludzie zwracają uwagę, że to rozkład sugeruje to założenie. Oczywiście, kiedy ludzie go używają, prawie * nigdy nie zamierzają, aby ich dane pochodziły z tego konkretnego rozkładu - jest to po prostu przybliżony opis relacji między ich średnią a wariancją. (Może to mieć sens przy bardzo prostych założeniach w niektórych wnioskach ubezpieczeniowych - całkowity koszt roszczenia, gdzie liczba roszczeń wynosi Poissona, a koszt roszczenia jest faktycznie stały.)
Glen_b

Odpowiedzi:

24

Model quasi-Poissona nie jest modelem pełnego prawdopodobieństwa (ML), ale modelem quasi-ML. Wystarczy użyć funkcji estymacji (lub funkcji score) z modelu Poissona, aby oszacować współczynniki, a następnie zastosować pewną funkcję wariancji, aby uzyskać odpowiednie błędy standardowe (lub raczej pełną macierz kowariancji) do przeprowadzenia wnioskowania. Dlatego glm()nie dostarcza i logLik()ani AIC()tutaj itp.

Jak słusznie zaznaczasz, model z tą samą funkcją oczekiwania i wariancji może zostać osadzony w strukturze ujemnej dwumianowej (NB), jeśli sizeparametr zmienia się wraz z oczekiwaniem . W literaturze nazywa się to zwykle parametryzacją NB1. Zobacz na przykład książkę Cameron i Trivedi (Analiza regresji danych zliczania) lub „Analiza mikrodanych” Winkelmann & Boes.μ iθiμi

Jeśli nie ma regresory (tylko osią) parametryzacji NB1 i parametryzacji NB2 zatrudnionych MASS„s glm.nb()zbieżne. Z regresorami różnią się. W literaturze statystycznej częściej używana jest parametryzacja NB2, ale niektóre pakiety oprogramowania oferują również wersję NB1. Na przykład w R możesz użyć gamlsspakietu do zrobienia gamlss(y ~ x, family = NBII). Zauważ, że nieco myląco gamlssużywa się NBIdo parametryzacji NB2 i NBIINB1. (Żargon i terminologia nie są jednakowe we wszystkich społecznościach).

W takim razie możesz oczywiście zapytać, dlaczego używać quasi-Poissona, jeśli jest dostępna NB1? Nadal istnieje subtelna różnica: poprzednia używa quasi-ML i uzyskuje oszacowanie z dyspersji od kwadratowych odchyleń (lub Pearsona) reszt. Ten ostatni używa pełnego ML. W praktyce różnica często nie jest duża, ale motywacje do korzystania z obu modeli są nieco inne.

Achim Zeileis
źródło
1
Dzięki! Bardzo pomocna odpowiedź, eksperymentuję gamlssteraz i wygląda na to, że właśnie tego potrzebowałem. Czy mógłbyś rozwinąć motywację użycia quasi-prawdopodobieństwa w porównaniu do pełnego ML?
user28400,
2
Zakładasz mniej: po prostu zakładasz (1) logarytmiczno-liniowy związek między oczekiwaniem a regresorami (2) liniowy związek między wariancją a oczekiwaniem. Reszta prawdopodobieństwa pozostaje całkowicie nieokreślona. Jako alternatywa dla (2), praktykujący czasami stosują tak zwane „solidne” standardowe błędy kanapkowe, które pozwoliłyby na bardziej ogólne wzorce heteroskedastyczności. Oczywiście można również zastosować NB1 z błędami standardowymi kanapek ... Kilka innych komentarzy znajduje się w naszym vignette("countreg", package = "pscl").
Achim Zeileis,