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ą
który może być wyposażony w R z zastosowaniem glm(..,family=poisson)
a glm.nb(...)
, odpowiednio. Jest też quasipoisson
rodzina, która w moim rozumieniu jest dostosowanym Poissonem o tym samym EV i wariancji
,
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:
,
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 P
- Czy ten pomysł ma sens? Czy moja weryfikacja opiera się na okrągłym uzasadnieniu?
- 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
źródło
Odpowiedzi:
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 ilogLik()
aniAIC()
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θja μja
size
parametr 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.μ iJeśli nie ma regresory (tylko osią) parametryzacji NB1 i parametryzacji NB2 zatrudnionych
MASS
„sglm.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ćgamlss
pakietu do zrobieniagamlss(y ~ x, family = NBII)
. Zauważ, że nieco mylącogamlss
używa sięNBI
do parametryzacji NB2 iNBII
NB1. (Ż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.
źródło
gamlss
teraz 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?vignette("countreg", package = "pscl")
.