Krótka odpowiedź
Nadmierna dyspersja nie ma znaczenia przy szacowaniu wektora współczynników regresji dla średniej warunkowej w modelu quasi / poissona! Będzie Ci dobrze, jeśli zapomnisz o nadmiernej dyspersji tutaj, użyj glmnet z rodziną poissonów i po prostu skup się na tym, czy twój błąd predykcji jest niski.
Kwalifikacja następuje poniżej.
Funkcje Poissona, Quasi-Poissona i funkcje szacowania:
Mówię to powyżej, ponieważ nadmierna dyspersja (OD) w modelu poissona lub quasi-poissona wpływa na wszystko, co ma związek z dyspersją (lub wariancją, skalą lub niejednorodnością lub rozprzestrzenianiem się, czy jakkolwiek to nazwać) i jako taka ma wpływ na standard błędy i przedziały ufności, ale pozostawia szacunki średniej warunkowej (zwanej ) nietknięte. Dotyczy to szczególnie liniowych rozkładów średnich, takich jakyμx⊤β .
Wynika to z faktu, że równania szacunkowe dla współczynników średniej warunkowej są praktycznie takie same dla modeli poissona i quasi-poissona. Quasi-poisson określa funkcję wariancji w kategoriach średniej i dodatkowego parametru (powiedz ) jako (z dla Poissona = 1), ale nie okazuje się mieć znaczenie przy optymalizacji równania szacunkowego. Zatem nie odgrywa żadnej roli w szacowaniu gdy średnia warunkowa i wariancja są proporcjonalne. Dlatego oszacowania punktowe są identyczne dla modeli quasi i Poissona!θVar(y)=θμθθθββ^
Pozwól mi zilustrować przykładem (zauważ, że trzeba przewinąć, aby zobaczyć cały kod i wynik):
> library(MASS)
> data(quine)
> modp <- glm(Days~Age+Sex+Eth+Lrn, data=quine, family="poisson")
> modqp <- glm(Days~Age+Sex+Eth+Lrn, data=quine, family="quasipoisson")
> summary(modp)
Call:
glm(formula = Days ~ Age + Sex + Eth + Lrn, family = "poisson",
data = quine)
Deviance Residuals:
Min 1Q Median 3Q Max
-6.808 -3.065 -1.119 1.819 9.909
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.71538 0.06468 41.980 < 2e-16 ***
AgeF1 -0.33390 0.07009 -4.764 1.90e-06 ***
AgeF2 0.25783 0.06242 4.131 3.62e-05 ***
AgeF3 0.42769 0.06769 6.319 2.64e-10 ***
SexM 0.16160 0.04253 3.799 0.000145 ***
EthN -0.53360 0.04188 -12.740 < 2e-16 ***
LrnSL 0.34894 0.05204 6.705 2.02e-11 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 2073.5 on 145 degrees of freedom
Residual deviance: 1696.7 on 139 degrees of freedom
AIC: 2299.2
Number of Fisher Scoring iterations: 5
> summary(modqp)
Call:
glm(formula = Days ~ Age + Sex + Eth + Lrn, family = "quasipoisson",
data = quine)
Deviance Residuals:
Min 1Q Median 3Q Max
-6.808 -3.065 -1.119 1.819 9.909
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.7154 0.2347 11.569 < 2e-16 ***
AgeF1 -0.3339 0.2543 -1.313 0.191413
AgeF2 0.2578 0.2265 1.138 0.256938
AgeF3 0.4277 0.2456 1.741 0.083831 .
SexM 0.1616 0.1543 1.047 0.296914
EthN -0.5336 0.1520 -3.511 0.000602 ***
LrnSL 0.3489 0.1888 1.848 0.066760 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for quasipoisson family taken to be 13.16691)
Null deviance: 2073.5 on 145 degrees of freedom
Residual deviance: 1696.7 on 139 degrees of freedom
AIC: NA
Number of Fisher Scoring iterations: 5
Jak widać, mimo że w tym zbiorze danych mamy silną naddyspersję 12,21, deviance(modp)/modp$df.residual
współczynniki regresji (szacunki punktowe) w ogóle się nie zmieniają. Ale zauważ, jak zmieniają się standardowe błędy.
Pytanie o efekt nadmiernej dyspersji w karanych modelach Poissona
Modele karane są najczęściej używane do przewidywania i wyboru zmiennych, a nie (jeszcze) do wnioskowania. Tak więc ludzie, którzy używają tych modeli, są zainteresowani parametrami regresji średniej warunkowej, po prostu skurczyli się do zera. Jeśli kara jest taka sama, równania szacunkowe dla środków warunkowych pochodzących z karanego (quasi) prawdopodobieństwa również nie zależą od a zatem nadmierna dyspersja nie ma znaczenia dla oszacowań w modelu typu:θβ
g(μ)=x⊤β+f(β)
ponieważ jest szacowany w ten sam sposób dla każdej funkcji wariancji formy , więc ponownie dla wszystkich modeli, w których średnia warunkowa i wariancja są proporcjonalne. βθμTo jest tak, jak w nie zdecentralizowanych modelach poissona / quasipoissona.
Jeśli nie chcesz brać tego za wartość nominalną i unikać matematyki, możesz znaleźć empiryczne wsparcie w tym, że w glmnet
, jeśli ustawisz parametr regularyzacji na 0 (a zatem ), to skończysz właściwie tam, gdzie lądują modele poissona i quasipoissona (patrz ostatnia kolumna poniżej, gdzie lambda wynosi 0,005).f(β)=0
> library(glmnet)
> y <- quine[,5]
> x <- model.matrix(~Age+Sex+Eth+Lrn,quine)
> modl <- glmnet(y=y,x=x, lambda=c(0.05,0.02,0.01,0.005), family="poisson")
> coefficients(modl)
8 x 4 sparse Matrix of class "dgCMatrix"
s0 s1 s2 s3
(Intercept) 2.7320435 2.7221245 2.7188884 2.7172098
(Intercept) . . . .
AgeF1 -0.3325689 -0.3335226 -0.3339580 -0.3340520
AgeF2 0.2496120 0.2544253 0.2559408 0.2567880
AgeF3 0.4079635 0.4197509 0.4236024 0.4255759
SexM 0.1530040 0.1581563 0.1598595 0.1607162
EthN -0.5275619 -0.5311830 -0.5323936 -0.5329969
LrnSL 0.3336885 0.3428815 0.3459650 0.3474745
Co więc OD robi z modelami regresji karanej? Jak zapewne wiesz, wciąż trwa debata na temat właściwego sposobu obliczania standardowych błędów dla modeli penalizowanych (patrz np. Tutaj ) i glmnet
tak czy inaczej nie jest generowana, prawdopodobnie z tego powodu. Może się zdarzyć, że OD wpłynie na część wnioskowania modelu, tak jak ma to miejsce w przypadku bez kary, ale dopóki nie zostanie osiągnięty konsens w sprawie wnioskowania w tym przypadku, nie będziemy wiedzieć.
Nawiasem mówiąc, można pozostawić cały ten bałagan, jeśli ktoś chce przyjąć pogląd bayesowski, w którym modele karane są jedynie standardowymi modelami z konkretnym przełożeniem.
poisson
aquasipoisson
regresje oceniają współczynniki w ten sam sposób, a różnią się tym, jak oceniają błędy standardowe, a tym samym znaczenie. Jednak w przypadku metody Lasso sposób obliczania błędów standardowych nie osiągnął jeszcze konsensusu, a zatem jego obecne zastosowanie polega głównie na selekcji zmiennych, a nie wnioskowaniu. Jako takie nie ma znaczenia, czy używamyglmnet
poissona, czy quasipoissona, ale co to znaczy, że błąd sprawdzany krzyżowo powinien zostać zminimalizowany.summary(modqp)
sam i zobaczyłem, że ma dokładnie takie same oszacowania współczynników. Wierzę, że twoja odpowiedź przyniesie korzyść większej liczbie osób w tej sprawie, ponieważ nie znalazłem żadnej, więc sugeruję dodanie wyniku podsumowania (modqp) w celu jeszcze lepszego zilustrowania przykładu. Jeszcze raz wielkie dzięki!