Strona GR

4.A Models for Over-Dispersed Count Data

Użyjemy danych z Longa (1990) dotyczących liczby publikacji tworzonych przez biochemików z tytułem doktora, aby zilustrować zastosowanie modeli Poissona, nadmiernie rozproszonego Poissona, ujemnego dwumianu i Poissona z zerową inflacją.

Zmiennymi w zbiorze danych są

.

zmienna	opis
art	artykuły w ostatnich trzech latach doktoratu.
fem	kodowany jeden dla kobiet
mar	kodowany jeden, jeśli zamężna
kid5	liczba dzieci poniżej szóstego roku życia
phd	prestiż doktoratu.D. program
ment	artykuły przez mentora w ostatnich trzech latach

Dane te zostały również przeanalizowane przez Longa i Freese’a (2001) i są dostępne na stronie internetowej Stata:

> library(foreign)> ab names(ab) "art" "fem" "mar" "kid5" "phd" "ment"> r c(mean=r, var=r, ratio=r/r) mean var ratio 1.692896 3.709742 2.191358 

Średnia liczba artykułów wynosi 1,69, a wariancja 3,71, czyli nieco więcej niż dwukrotność średniej. Dane są nadmiernie rozproszone, ale oczywiście nie uwzględniliśmy jeszcze żadnych zmiennych.

Model Poissona

Zastosujmy model używany przez Longa i Freese’a(2001), prosty model addytywny wykorzystujący wszystkie pięć predyktorów.

> mp summary(mp)Call:glm(formula = art ~ fem + mar + kid5 + phd + ment, family = poisson, data = ab)Deviance Residuals: Min 1Q Median 3Q Max -3.5672 -1.5398 -0.3660 0.5722 5.4467 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.304617 0.102981 2.958 0.0031 ** femWomen -0.224594 0.054613 -4.112 3.92e-05 ***marMarried 0.155243 0.061374 2.529 0.0114 * kid5 -0.184883 0.040127 -4.607 4.08e-06 ***phd 0.012823 0.026397 0.486 0.6271 ment 0.025543 0.002006 12.733 
We see that the model obviously doesn't fit the data. The five-percent critical value for a chi-squared with 909 d.f. is
>qchisq(0.95, df.residual(mp)) 980.2518> deviance(mp) 1634.371> pr sum(pr^2) 1662.547
i dewiancja i chi kwadrat Pearsona są w 1600s.
 4.A.1. Extra-Poisson Variation
Założymy teraz, że wariancja jest proporcjonalna, a nie równa średniej, i oszacujemy parametr skali φ dzieląc chi kwadrat Pearsona przez jego d.f.:
> phi round(c(phi,sqrt(phi)),4) 1.8290 1.3524
Widzimy, że wariancja jest o około 83% większa od średniej. Oznacza to, że powinniśmy skorygować błędy standardowe, mnożąc je przez 1,35, pierwiastek kwadratowy z 1,83.
R może wykonać to obliczenie za nas, jeśli użyjemy rodziny quasipoisson:
> mq summary(mq)Call:glm(formula = art ~ fem + mar + kid5 + phd + ment, family = quasipoisson, data = ab)Deviance Residuals: Min 1Q Median 3Q Max -3.5672 -1.5398 -0.3660 0.5722 5.4467 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.304617 0.139273 2.187 0.028983 * femWomen -0.224594 0.073860 -3.041 0.002427 ** marMarried 0.155243 0.083003 1.870 0.061759 . kid5 -0.184883 0.054268 -3.407 0.000686 ***phd 0.012823 0.035700 0.359 0.719544 ment 0.025543 0.002713 9.415 
The estimates are exactly the same as before, but the standard errorsare about 35% larger. We can verify this fact easily. First we writea useful function to extract standard errors and then use it onour fits:
> se round(data.frame(p=coef(mp), q=coef(mq),+ se.p=se(mp), se.q=se(mq), ratio=se(mq)/se(mp)), 4) p q se.p se.q ratio(Intercept) 0.3046 0.3046 0.1030 0.1393 1.3524femWomen -0.2246 -0.2246 0.0546 0.0739 1.3524marMarried 0.1552 0.1552 0.0614 0.0830 1.3524kid5 -0.1849 -0.1849 0.0401 0.0543 1.3524phd 0.0128 0.0128 0.0264 0.0357 1.3524ment 0.0255 0.0255 0.0020 0.0027 1.3524
Alternatywnym podejściem jest dopasowanie modelu Poissona i użycie estymatora błędów standardowych therobust lub sandwich. Zwykle daje to wyniki bardzo zbliżone do modelu Poissona z nadmiernym rozproszeniem.
 4.A.2 Negative Binomial Regression
Zastosujemy teraz model dwumianowy ujemny z tymi samymi predyktorami. W tym celu potrzebujemy funkcji glm.nb() z pakietuMASS.
> library(MASS)> mnb summary(mnb)Call:glm.nb(formula = art ~ fem + mar + kid5 + phd + ment, data = ab, init.theta = 2.264387695, link = log)Deviance Residuals: Min 1Q Median 3Q Max -2.1678 -1.3617 -0.2806 0.4476 3.4524 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.256144 0.137348 1.865 0.062191 . femWomen -0.216418 0.072636 -2.979 0.002887 ** marMarried 0.150489 0.082097 1.833 0.066791 . kid5 -0.176415 0.052813 -3.340 0.000837 ***phd 0.015271 0.035873 0.426 0.670326 ment 0.029082 0.003214 9.048 1/mnb$theta 0.4416205> -2*(logLik(mp)-logLik(mnb))'log Lik.' 180.196 (df=6)
R's thetajest precyzją efektu multiplikatywnego losowego i odpowiada 1/σ2 w węzłach. Oszacowanie odpowiada szacowanej wariancji 0,44 i jest wysoce istotne.
Aby sprawdzić istotność tego parametru można pomyśleć o obliczeniu dwukrotnej różnicy log-likelihoods pomiędzy tym modelem a modelem Poissona, 180,2, i potraktowaniu jej jako chi-squared z jednym d.f. Zwykła asymptotyka nie ma jednak zastosowania, ponieważ hipoteza zerowa znajduje się na granicy przestrzeni parametrów.
Istnieją pewne prace pokazujące, że lepszym przybliżeniem jest traktowanie statystyki jako mieszanki 50:50 zera i chi kwadrat z jednym d.f. Alternatywnie, traktowanie statystyki jako chi kwadrat z jednym d.f. daje konserwatywny test. Tak czy inaczej, mamy przytłaczający dowód na nadmierną dyspersję.
Do testowania hipotez o współczynnikach regresji możemy użyć testów Walda lub testów współczynnika prawdopodobieństwa, które są możliwe, ponieważ przyjęliśmy pełne założenia dystrybucyjne.
Istnieje również rodzina negative.binomial dla glm i można jej użyć pod warunkiem, że podane są parametrytheta. (Jest to oparte na wyniku, że ujemny dwumian jest w rodzinie glm dla stałej wariancji). Oto szybkie sprawdzenie:
> mnbg all(abs(coef(mnb)==coef(mnbg)) 
The standard errors would differ, however, because  allows for the fact that was estimated, whereas  does not.
Unobserved Heterogeneity
R has a function  to compute the density of a gamma distribution with given shape and scale (or its reciprocal the rate). In particular, when the random effect hasvariance  the density is. This can beused to draw the density.
> v = 1/mnb$theta> gd = data.frame(x = seq(0.0, 3, 0.1), g = dgamma(x, shape = 1/v, scale = v))> ggplot(gd, aes(x, g)) + geom_line()> ggsave("gamdenr.png", width=500/72, height=400/72, dpi=72)

 Możemy również obliczyć kwantyle używającqgamma(p, shape, rate = 1, scale = 1/rate, lower.tail = TRUE). W szczególności, gdy efekt losowy ma wariancję v kwartyle są qgamma((1:3)/4, shap e= 1/v, scale = v).
> qgamma((1:3)/4, shape = 1/v, scale = v) 0.5114167 0.8572697 1.3347651
Biochemicy z Q1 rozkładu nieobserwowanej heterogeniczności publikują o 49% mniej prac niż oczekiwano na podstawie ich obserwowanych cech, podczas gdy ci z mediany publikują o 14% mniej, a ci z Q3 publikują o 33% więcej niż oczekiwano.
Porównanie oszacowań i błędów standardowych
Porównajmy oszacowania parametrów i błędy standardowe w modelach Poissona, nadmiernie rozproszonego Poissona i ujemnego dwumianu:
> round(data.frame(+ p=coef(mp),q=coef(mq),nb=coef(mnb),+ se.p=se(mp),se.q=se(mq),se.nb=se(mnb)),4) p q nb se.p se.q se.nb(Intercept) 0.3046 0.3046 0.2561 0.1030 0.1393 0.1373femWomen -0.2246 -0.2246 -0.2164 0.0546 0.0739 0.0726marMarried 0.1552 0.1552 0.1505 0.0614 0.0830 0.0821kid5 -0.1849 -0.1849 -0.1764 0.0401 0.0543 0.0528phd 0.0128 0.0128 0.0153 0.0264 0.0357 0.0359ment 0.0255 0.0255 0.0291 0.0020 0.0027 0.0032
Oszacowania ujemnego dwumianu nie różnią się zbytnio od oszacowań opartych na modelu Poissona, a oba zestawy doprowadziłyby do tych samych wniosków. 
Patrząc na błędy standardowe, widzimy, że oba podejścia do nadmiernego rozproszenia prowadzą do bardzo podobnych szacunkowych błędów standardowych, a zwykła regresja Poissona zaniża błędy standardowe.
Dobroć dopasowania
Możemy ocenić dobroć dopasowania ujemnego modelu dwumianowego za pomocą odchylenia
> deviance(mnbg) 1004.281
Ujemny model dwumianowy pasuje lepiej niż model Poissona, ale nadal ma odchylenie powyżej pięcioprocentowej wartości krytycznej980.25.
Funkcja wariancji
Nadmiernie rozproszony model Poissona i ujemny model dwumianowy mają różne funkcje wariancji. Jednym ze sposobów sprawdzenia, który z nich może być bardziej odpowiedni, jest utworzenie grup opartych na predyktorze liniowym, obliczenie średniej i wariancji dla każdej grupy, a następnie wykreślenie zależności średnia-wariancja.
Tutaj znajdują się grupy oparte na liniowym predyktorze dwumianowym ujemnym, utworzone przy użyciu cutz przerwami na 5(5)95 percentylu, aby uzyskać dwadzieścia grup o w przybliżeniu równej wielkości. Zauważ, że predict() oblicza predyktor liniowy, chyba że określono inaczej. Aby przewidywać w skali odpowiedzi, dodaj opcję type="response". 
> xb g m v png("c4afig1.png",width=500,height=400)> plot(m, v, xlab="Mean", ylab="Variance", + main="Mean-Variance Relationship")> mtext("Articles Published by Ph.D. Biochemists",padj=-0.5)> x lines(x, x*phi, lty="dashed")> lines(x, x*(1+x/mnb$theta))> legend("topleft", lty=c("dashed","solid"), + legend=c("Q. Poisson","Neg. Binom."), inset=0.05)> dev.off();null device 1 

Wykres przedstawia średnią względem wariancji i nakłada krzywe odpowiadające nadmiernie rozproszonemu modelowi Poissona, w którym wariancja wynosi φμ, oraz ujemnemu modelowi dwumianowemu, w którym wariancja wynosi μ(1+μσ2).
Funkcja wariancji Poissona wykonuje całkiem dobrą robotę dla większości danych, ale nie jest w stanie uchwycić wysokich wariancji najbardziej produktywnych uczonych. The negative binomial wariancja funkcja być zbyt różny ale, być kwadratowy, móc wzrastać szybki i robić lepszy praca przy the wysoki końcówka. Wnioskujemy, że ujemny model dwumianowy zapewnia lepszy opis danych niż nadmiernie rozproszony model Poissona.
 4.A.3 Modele Poissona z zerem
 Częstym zjawiskiem w przypadku danych zliczeniowych jest nadmiar zer w porównaniu z tym, co jest oczekiwane w modelu Poissona. Jest to faktycznie problem z naszymi danymi: 
> zobs zpoi c(obs=mean(zobs), poi=mean(zpoi)) obs poi 0.3005464 0.2092071 
 Widzimy, że 30,0% naukowców w próbie nie opublikowało żadnego artykułu w ciągu ostatnich trzech lat od doktoratu, ale model Poissona przewiduje, że tylko 20,9% nie miałoby żadnych publikacji. Najwyraźniej model ten niedoszacowuje prawdopodobieństwa zerowych zliczeń.
Jednym ze sposobów modelowania tego typu sytuacji jest założenie, że dane pochodzą z mieszaniny dwóch populacji, jednej, w której zliczenia są zawsze zerowe, i drugiej, w której zliczenia mają rozkład Poissona o średniej μ. W tym modelu zero zliczeń może pochodzić z każdej z tych populacji, podczas gdy dodatnie zliczenia pochodzą tylko z drugiej.
W kontekście publikacji doktorów biochemii, możemy sobie wyobrazić, że niektórzy z nich mieli na myśli pracę, w której publikacje nie byłyby ważne, podczas gdy inni celowali w pracę akademicką, gdzie oczekiwano rekordu publikacji. Członkowie pierwszej grupy opublikowaliby zero artykułów, podczas gdy członkowie drugiej grupy opublikowaliby 0,1,2,..., liczbę, która może być przyjęta jako rozkład Poissona.
Rozkład wyniku może być modelowany w kategoriach dwóch parametrów, π prawdopodobieństwo "zawsze zero", i μ, średnia liczba publikacji dla tych, którzy nie są w grupie "zawsze zero". Naturalnym sposobem wprowadzenia kowariancji jest modelowanie logitu prawdopodobieństwa π zawsze zerowego i logu średniej μ dla tych, którzy nie są w klasie zawsze zerowej.
Ten typ modelu może być dopasowany w R używając funkcji zeroinfl() w pakiecie pscl. Wzór modelu może być określony jak zwykle, jeśli te same zmienne mają być uwzględnione w obu równaniach. W przeciwnym razie można podać dwa zestawy predyktorów oddzielone pionowym paskiem, wpisz ?zeroinfl, aby dowiedzieć się więcej.
Tutaj model z zerową inflacją ze wszystkimi zmiennymi w obu równaniach:
> library(pscl)> mzip summary(mzip)Call:zeroinfl(formula = art ~ fem + mar + kid5 + phd + ment, data = ab)Pearson residuals: Min 1Q Median 3Q Max -2.3253 -0.8652 -0.2826 0.5404 7.2976 Count model coefficients (poisson with log link): Estimate Std. Error z value Pr(>|z|) (Intercept) 0.640838 0.121307 5.283 1.27e-07 ***femWomen -0.209145 0.063405 -3.299 0.000972 ***marMarried 0.103751 0.071111 1.459 0.144565 kid5 -0.143320 0.047429 -3.022 0.002513 ** phd -0.006166 0.031008 -0.199 0.842378 ment 0.018098 0.002294 7.888 3.07e-15 ***Zero-inflation model coefficients (binomial with logit link): Estimate Std. Error z value Pr(>|z|) (Intercept) -0.577059 0.509386 -1.133 0.25728 femWomen 0.109746 0.280082 0.392 0.69518 marMarried -0.354014 0.317611 -1.115 0.26502 kid5 0.217097 0.196482 1.105 0.26919 phd 0.001274 0.145263 0.009 0.99300 ment -0.134114 0.045243 -2.964 0.00303 **---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Number of iterations in BFGS optimization: 21 Log-likelihood: -1605 on 12 Df
Patrząc na równanie inflate widzimy, że jedynym znaczącym predyktorem bycia w klasie "zawsze zerowej" jest liczba artykułów opublikowanych przez mentora, z każdym artykułem mentora związanym z 12,6% niższym prawdopodobieństwem niepublikowania.
Patrząc na równanie dla średniej liczby artykułów wśród tych, którzy nie znaleźli się w klasie zawsze zerowej, znajdujemy znaczące niedogodności dla kobiet i naukowców z dziećmi poniżej piątego roku życia oraz duży pozytywny efekt liczby publikacji mentora, z każdym artykułem związanym z 1,8% wzrostem oczekiwanej liczby publikacji.
Aby sprawdzić, czy model rozwiązuje problem nadmiaru zer, przewidujemy π i μ oraz obliczamy łączne prawdopodobieństwo braku publikacji. Do ich uzyskania służą opcje w funkcji predict() o nazwach "zero" i "count". Istnieje również opcja "prob" do obliczania przewidywanej gęstości, ale jest to zbyteczne, ponieważ chcemy tylko prawdopodobieństwo zera.
> pr mu zip mean(zip) 0.2985679
Model rozwiązuje więc problem nadmiaru zer, przewidując, że 29,9% biochemików nie opublikuje żadnych artykułów, co jest znacznie bliższe obserwowanej wartości 30,0%.
Istnieje również model ujemnego dwumianu z wypełnieniem zerowym, który wykorzystuje ujemny dwumian dla liczby w klasie "nie zawsze zero". Model ten można dopasować za pomocą zeroinfl() z parametrem dist="negbin". Alternatywne połączenia dla równania inflate obejmują probit, który można określić przy użyciu link="probit". 
Porównanie modeli z AIC
Tak się składa, że dla tych danych ujemny dwumian również rozwiązuje problem! Oto prawdopodobieństwo zera artykułów w ujemnym dwumianowym. Postępujemy z pierwszych zasad, ale można też użyć wbudowanej gęstości dwumianu ujemnego
> munb theta znb # also dnbinom(0, mu=munb, size=theta)> mean(znb) 0.3035957
Model przewiduje, że 30,4% biochemików nie opublikuje żadnego artykułu w ciągu ostatnich trzech lat doktoratu, co jest bardzo bliskie zaobserwowanej wartości 30,0%.
Aby wybrać między dwumianem ujemnym a modelem z zerowym prawdopodobieństwem, musimy uciec się do innych kryteriów. Bardzo prostym sposobem na porównanie modeli o różnej liczbie parametrów jest obliczenie kryterium informacyjnego Akaike'a (AIC), które definiujemy jako
AIC = -2logL + 2p
gdzie p jest liczbą parametrów w modelu. Pierwszy człon jest zasadniczo odchyleniem, a drugi karą za liczbę parametrów, i można go obliczyć dla większości modeli za pomocą funkcji AIC. Można go obliczyć dla większości modeli za pomocą funkcji AIC. Uzyskuję go również "ręcznie", aby móc zweryfikować obliczenia. Dla naszych danych
> c(nb=AIC(mnb), zip=AIC(mzip)) nb zip 3135.917 3233.546 > mzip$rank aic sapply(list(mnb,mzip), aic) 3133.917 3233.546
 Dla tego zbioru danych model dwumianowy ujemny jest wyraźnym zwycięzcą pod względem parsymonii i dobroci dopasowania. 
Inne kryteria diagnostyczne, na które moglibyśmy spojrzeć, to rozkład krańcowy przewidywanych i obserwowanych liczebności oraz funkcje wariancji. 
Modele zero-trójkątne i płotkowe
Inne modele, których nie omówiliśmy, to zero-trójkątny Poisson i ujemny dwumian, zaprojektowane dla danych, które nie zawierają zer. Typowym przykładem jest długość pobytu w szpitalu, która wynosi co najmniej jeden dzień. Rozsądnym podejściem jest dopasowanie modelu Poissona lub ujemnego dwumianu, który wyklucza zero i przeskalowuje pozostałe prawdopodobieństwa tak, aby sumowały się do jednego. Należy zachować ostrożność przy interpretacji tych modeli, ponieważ μ nie jest oczekiwanym wynikiem, ale średnią rozkładu bazowego, który zawiera zera. Modele te są zaimplementowane w funkcji vglm() w pakiecie VGAM, używając rodzin pospoisson i posnegbinomial.
Alternatywnym podejściem do nadmiaru (lub niedoboru) zer jest użycie procesu dwustopniowego, z modelem logitowym do rozróżnienia między zliczeniami zerowymi i dodatnimi, a następnie z obciętym do zera modelem Poissona lub ujemnym modelem dwumianowym dla zliczeń dodatnich. W naszym przykładzie moglibyśmy użyć modelu logitowego do rozróżnienia tych, którzy publikują od tych, którzy nie publikują, a następnie obciętego modelu Poissona lub ujemnego modelu dwumianowego dla liczby artykułów tych, którzy publikują co najmniej jeden. Modele te są często nazywane modelami przeszkód. Można je dopasować w R używając oddzielnych modeli logitowych i zero-truncatedPoisson lub negative binomial i po prostu dodając log-likelihoods, lub używając funkcji hurdle() w pakiecie pscl.
Porównując modele płotkowe i modele z zerowym wypełnieniem, uważam, że rozróżnienie między zerem a jednym lub więcej jest bardziej przejrzyste w przypadku modeli płotkowych, ale interpretacja średniej jest bardziej przejrzysta w przypadku modeli z zerowym wypełnieniem.
Poniższy dokument JSS zawiera użyteczną dyskusję na temat wszystkich tych modeli: Regression Models for Count Data in R.
.