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

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ą$\mu$

$Var_P = \mu$

$Var_{NB} = \mu + \frac{\mu^2}{\theta}$

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

$Var_{QP} = \phi\mu$ ,

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:$\phi = 1 + \frac{\mu}{\theta}$

$QP\,(\mu,\phi) = NB\,(\mu,\theta = \frac{\mu}{\phi-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 :$\theta$$\mu$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:$\phi$

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)$\phi$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 $AIC_{QP} < AIC_{NB}$$AIC_{QP}$

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

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.μ $\theta_i$$\mu_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.