GR honlapja

4.A Modellek túlszórt számlálási adatokra

Long (1990) adatait használjuk a PhD biokémikusok által készített publikációk számáról, hogy szemléltessük a Poisson, a túlszórt Poisson, a negatív binomiális és a nulla-inflált Poisson modellek alkalmazását.

Az adathalmazban szereplő változók a következők

art articles in last three years of Ph.D.

személyek esetén eggyel kódolva

mar személyek esetén eggyel kódolva. kid5 hat év alatti gyermekek száma phd presztízs a Ph.D. program

Variable	Description
fem
ment	cikkek mentor által az elmúlt három évben

Ezeket az adatokat Long és Freese (2001) is elemezte, és a Stata honlapjáról elérhetőek:

> 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 

A cikkek átlagos száma 1,69, a szórás pedig 3,71, ami valamivel több mint kétszerese az átlagnak. Az adatok túlszórtak, de természetesen még nem vettünk figyelembe semmilyen kovariátort.

A Poisson-modell

Feltesszük a Long és Freese(2001) által használt modellt, egy egyszerű additív modellt, amely mind az öt prediktort használja.

> 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
és a deviancia és a Pearson-féle chi-négyzet egyaránt az 1600-as értékek közé esik.
 4.A.1. Extra-Poisson-variáció
Most feltételezzük, hogy a variancia nem egyenlő, hanem arányos az átlaggal, és a φ skálaparamétert úgy becsüljük meg, hogy a Pearson-féle chi-négyzetet elosztjuk annak d.f.-jével:
> phi round(c(phi,sqrt(phi)),4) 1.8290 1.3524
Azt látjuk, hogy a variancia körülbelül 83%-kal nagyobb, mint az átlag. Ez azt jelenti, hogy a standard hibákat 1,35-tel, azaz 1,83 négyzetgyökével megszorozva kell kiigazítanunk.
R elvégezheti helyettünk ezt a számítást, ha a quasipoissoncsaládot használjuk:
> 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
Egy alternatív megközelítés a Poisson-modell illesztése és a standard hibák therobust vagy szendvicsbecslőjének használata. Ez általában nagyon hasonló eredményeket ad, mint a túlszórt Poisson-modell.
 4.A.2 Negatív binomiális regresszió
Most egy negatív binomiális modellt illesztünk ugyanezekkel a prediktorokkal.Ehhez a MASS csomag glm.nb() függvényére van szükségünk.
> 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 theta a multiplikatív véletlen hatás pontossága, és 1/σ2-nek felel meg az akkoriban. A becslés0,44-es becsült varianciának felel meg, és rendkívül szignifikáns.
Ez a paraméter szignifikanciájának teszteléséhez úgy gondolhatod, hogy kiszámítod a log-likelihoodok különbségének kétszeresét e modell és a Poisson-modell között, 180,2-t, és ezt egy d.f.-es chi-négyzetként kezeled. A szokásos aszimptotika azonban nem alkalmazható, mivel a nullhipotézis a paramétertér egy határán van.
Van néhány munka, amely azt mutatja, hogy jobb közelítés, ha a statisztikát a nulla és a chi négyzet egy d.f.-es chi-négyzet 50:50 arányú keverékeként kezeljük. Akárhogy is, elsöprő bizonyítékunk van a túldiszperzióra.
A regressziós együtthatókra vonatkozó hipotézisek tesztelésére használhatunk Wald-teszteket vagy likelihood ratio teszteket, amelyek azért lehetségesek, mert teljes eloszlási feltevéseket tettünk.
A glm-nak van egy negative.binomial családja is, és ez is használható, feltéve, hogy a paraméterektheta adottak. (Ez azon az eredményen alapul, hogya negatív binomiális fix variancia esetén a glm családban van.) Íme egy gyors ellenőrzés:
> 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)

Kvantiliseket is kiszámíthatunk aqgamma(p, shape, rate = 1, scale = 1/rate, lower.tail = TRUE) segítségével. Különösen, ha a véletlen hatás varianciája va kvartilisek 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
A megfigyeletlen heterogenitás eloszlásának Q1-ében lévő biokémikusok 49%-kal kevesebb publikációt publikálnak, mint a megfigyelt jellemzőik alapján várható, míg a mediánban lévők 14%-kal kevesebbet, a Q3-ban lévők pedig 33%-kal többet publikálnak a vártnál.
A becslések és standard hibák összehasonlítása
Hasonlítsuk össze a Poisson-, a túlszórt Poisson- és a negatív binomiális modell szerinti paraméterbecsléseket és standard hibákat:
> 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
A negatív binomiális becslések nem nagyon különböznek a Poisson-modellen alapuló becslésektől, és mindkét halmaz ugyanarra a következtetésre vezetne. 
A standard hibákat vizsgálva azt látjuk, hogy a túlszórás mindkét megközelítése nagyon hasonló becsült standard hibákhoz vezet,és hogy a közönséges Poisson-regresszió alulbecsüli a standard hibákat.
Az illeszkedés jósága
A negatív binomiális modell illeszkedésének jóságát a deviancia
> deviance(mnbg) 1004.281
A negatív binomiális modell jobban illeszkedik, mint a Poisson, de még mindig az öt százalékos kritikus érték980,25 feletti devianciával rendelkezik.
A szórásfüggvény
A túlszórt Poisson és a negatív binomiális modellnek különböző szórásfüggvényei vannak. Az egyik módja annak, hogy ellenőrizzük, melyik lehet megfelelőbb, hogy csoportokat hozunk létre a lineáris prediktor alapján, kiszámítjuk az egyes csoportok átlagát és szórását, és végül ábrázoljuk az átlag-variáció összefüggést.
Itt vannak a negatív binomiális lineáris prediktoron alapuló csoportok, amelyeket a cut segítségével hoztunk létre, az 5(5)95 percentilisnél lévő törésekkel, hogy húsz, közel azonos méretű csoportot hozzunk létre. Vegye figyelembe, hogy a predict() lineáris prediktort számol,hacsak nincs másképp megadva. A válasz skáláján történő előrejelzéshezadjuk hozzá a type="response" opciót. 
> 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 

A grafikon az átlagot ábrázolja a varianciával szemben, és egymásra helyezi a túlszórt Poisson-modellnek megfelelő görbéket,ahol a variancia φμ, és a negatív binomiálismodellt, ahol a variancia μ(1+μσ2).
A Poisson varianciafüggvény elég jó munkát végez az adatok nagy részénél, de nem képes megragadni a legtermékenyebb tudósok magas varianciáját. A negatív binomiális varianciafüggvény nem különbözik túlságosan, de mivel kvadratikus, gyorsabban emelkedhet, és jobb munkát végez a magas értékeknél. Arra a következtetésre jutunk, hogy a negatív binomiális modell jobban leírja az adatokat, mint a túlságosan szórt Poisson-modell. 
 4.A.3. Nullával felfújt Poisson-modellek
 A számlálási adatoknál gyakori jelenség a nullák túlsúlya ahhoz képest, ami a Poisson-modell szerint várható. Ez tulajdonképpen a mi adataink esetében jelent problémát: 
> zobs zpoi c(obs=mean(zobs), poi=mean(zpoi)) obs poi 0.3005464 0.2092071 
Láthatjuk, hogy a mintában szereplő tudósok 30,0%-a nem publikált cikket a doktori cím megszerzésének utolsó három évében, de a Poisson-modell azt jósolja, hogy csak 20,9%-nak nem lenne publikációja. A modell nyilvánvalóan alulbecsüli a nullás számok valószínűségét.
Az ilyen típusú helyzetek modellezésének egyik módja, hogy feltételezzük, hogy az adatok két populáció keverékéből származnak, az egyikből, ahol a számok mindig nullásak, a másikból, ahol a számoknak μ átlagú Poisson-eloszlása van. Ebben a modellben a nulla számlálás bármelyik populációból származhat, míg a pozitív számlálás csak a másodikból.
A doktori fokozatot szerzett biokémikusok publikációival összefüggésben elképzelhető, hogy egyesek olyan állásokra gondoltak, ahol a publikációk nem lennének fontosak, míg mások olyan tudományos állásokra törekedtek, ahol elvárt a publikációs rekord. Az első csoport tagjai nulla cikket publikálnának, míg a második csoport tagjai 0,1,2,... cikket publikálnának, amely számról feltételezhető, hogy Poisson-eloszlású.
A végeredmény eloszlása ekkor két paraméterrel modellezhető: π a "mindig nulla" valószínűsége, μ pedig a publikációk átlagos száma azok számára, akik nem tartoznak a "mindig nulla" csoportba. A kovariánsok bevezetésének természetes módja, ha a mindig nulla π valószínűségének és a mindig nulla csoportba nem tartozók μ átlagának logaritmusát modellezzük.
Egy ilyen típusú modell az R-ben a pscl csomag zeroinfl()függvényével illeszthető. A modellképletet a szokásos módon lehet megadni, ha ugyanazokat a változókat kell bevonni mindkét egyenletbe. Ellenkező esetben meg lehet adni két prediktorcsoportot, amelyeket egy függőleges sáv választ el egymástól, írja be a ?zeroinfl-t, hogy többet tudjon meg.
Itt van egy nullára felfújt modell, amelyben az összes kovariáns szerepel mindkét egyenletben:
> 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
A felfújt egyenletet vizsgálva azt látjuk, hogy a "mindig nulla" osztályba tartozás egyetlen jelentős prediktora a mentor által publikált cikkek száma, a mentor minden egyes cikke 12,6%-kal alacsonyabb esélyt jelent a soha nem publikálásra.
Ha megnézzük az átlagos cikkszámra vagy cikkekre vonatkozó egyenletet azok körében, akik nem tartoznak a mindig nulla osztályba, jelentős hátrányt találunk a nők és az öt év alatti gyermeket nevelő tudósok esetében, valamint a mentor által publikált cikkek számának nagy pozitív hatását: minden egyes cikk 1,8%-kal növeli a publikációk várható számát.
Azt ellenőrizendő, hogy a modell megoldja a túlzott nullák problémáját, megjósoljuk π és μ értékét, és kiszámítjuk a publikálás elmaradásának együttes valószínűségét. A predict()funkcióban vannak "zero" és "count" nevű opciók, amelyekkel ezeket megkaphatjuk. Van egy "prob" opció is a megjósolt sűrűség kiszámításához, de ez túlzás, mivel csak a nulla valószínűségére vagyunk kíváncsiak.
> pr mu zip mean(zip) 0.2985679
A modell tehát megoldja a többlet nullák problémáját, és azt jósolja, hogy a biokémikusok 29,9%-a nem publikál cikket, ami sokkal közelebb van a megfigyelt 30,0%-os értékhez.
Létezik egy nulla-inflált negatív binomiális modell is, amely a "nem mindig nulla" osztályba tartozó számhoz negatív binomot használ. Ez a modell a zeroinfl() segítségével illeszthető a dist="negbin" paraméterrel. A felfújt egyenlet alternatív kapcsolatai közé tartozik a probit,amely a link="probit" használatával adható meg. 
Modell-összehasonlítás AIC-vel
Ezekre az adatokra történetesen a negatív binomiális is megoldja a problémát! Íme a nulla cikkek valószínűsége a negatív binomiálisban. Első elvekből indulunk ki, de használhatnánk a beépített negatív binomiális sűrűséget is
> munb theta znb # also dnbinom(0, mu=munb, size=theta)> mean(znb) 0.3035957
A modell azt jósolja, hogy a biokémikusok 30,4%-a nem publikál cikket a doktorátusa utolsó három évében, ami nagyon közel van a megfigyelt 30,0%-os értékhez.
A negatív binomiális és a nullával felfújt modell közötti választáshoz más kritériumokhoz kell folyamodnunk. Egy nagyon egyszerű módja a különböző paraméterszámú modellek összehasonlításának az Akaike információs kritérium (AIC) kiszámítása, amelyet a következőképpen határozunk meg:
AIC = -2logL + 2p
ahol p a modell paramétereinek száma. Az első kifejezés lényegében az eltérés, a második pedig a paraméterek számának büntetése, és a legtöbb modellre kiszámítható a AIC függvény segítségével. Én is "kézzel" kapom meg, hogy ellenőrizni tudjuk a számítást. A mi adatainkra
> 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
 Erre az adatkészletre a negatív binomiális modell egyértelmű győztes a takarékosság és az illeszkedés jósága szempontjából. 
Más diagnosztikai kritériumok, amelyeket megvizsgálhatunk, a megjósolt és a megfigyelt számlálások marginális eloszlása és a varianciafüggvények. 
Zéró-csonkított és hurdle-modellek
Más modellek, amelyekkel még nem foglalkoztunk, a nulla-csonkított Poisson- és a negatív binomiális modellek, amelyeket olyan adatokra terveztek, amelyek nem tartalmaznak nullákat. Gyakori példa a kórházi tartózkodás hossza, amely legalább egy nap. Ésszerű megközelítés egy olyan Poisson- vagy negatív binomiális modell illesztése, amely kizárja a nullát, és a többi valószínűséget átméretezi, hogy azok összege egy legyen. Óvatosan kell értelmezni ezeket a modelleket, mert μ nem a várható eredmény, hanem egy olyan alapeloszlás átlaga, amely tartalmazza a nullákat. Ezeket a modelleket a VGAM csomag vglm() függvénye implementálja, a pospoisson és posnegbinomial családok használatával.
A nullák többletének (vagy hiányának) alternatív megközelítése a kétlépcsős folyamat használata, egy logit modellel a nulla és pozitív számok megkülönböztetésére, majd egy nullával lemetszett Poisson vagy negatív binomiális modellel a pozitív számokra. Példánkban használhatnánk egy logit-modellt, hogy megkülönböztessük azokat, akik publikálnak, azoktól, akik nem publikálnak, majd egy nullára csonkolt Poisson- vagy negatív binomiális modellt azok cikkszámára, akik legalább egyet publikálnak. Ezeket a modelleket gyakran nevezik akadálymodelleknek. Az R-ben a külön logit és a nullára csonkolt Poisson vagy negatív binomiális modellekkel és a log-likelihoodok egyszerű összeadásával, vagy a pscl csomag hurdle() függvényének használatával illeszthetők R-ben.
A gátló és a nulla-csonkított modellek összehasonlítása során úgy találom, hogy a nulla és az egy vagy több közötti különbségtétel a gátló modellek esetében világosabb, de az átlag értelmezése a nulla-csonkított modellek esetében egyértelműbb.
A következő JSS-dolgozat hasznos tárgyalást tartalmaz mindezen modellekről: Regression Models for Count Data in R.