GR’s Website

4.A Modely pro data s nadměrně rozptýlenými počty

Pro ilustraci použití Poissonova modelu, modelu s nadměrně rozptýlenými počty, záporného binomického modelu a Poissonova modelu s nulovou hodnotou použijeme údaje od Longa (1990) o počtu publikací vytvořených doktorandy biochemie.

Proměnné v souboru dat jsou

.

Proměnná	Popis
art	články za poslední tři roky doktorského studia.
fem	kódováno jedna u žen
mar	kódováno jedna, pokud jsou vdané
kid5	počet dětí do šesti let
phd	prestiž Ph.D. programu
ment	článků podle školitele v posledních třech letech

Tato data analyzovali také Long a Freese (2001) a jsou k dispozici na webových stránkách 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 

Průměrný počet článků je 1,69 a rozptyl je 3,71, tedy o něco více než dvojnásobek průměru. Data jsou nadměrně rozptýlená, ale samozřejmě jsme zatím neuvažovali žádné kovariáty.

Poissonův model

Posadíme model použitý Longem a Freesem(2001), jednoduchý aditivní model s použitím všech pěti prediktorů.

> 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
a odchylka i Pearsonův chí-kvadrát se pohybují v rozmezí 1600.
 4.A.1. Odchylka a Pearsonův chí-kvadrát se pohybují v rozmezí 1600. Mimopoissonová variance
Nyní předpokládáme, že rozptyl je úměrný spíše než roven průměru, a odhadneme parametr měřítka φ vydělením Pearsonova chí-kvadrátu jeho d.f.:
> phi round(c(phi,sqrt(phi)),4) 1.8290 1.3524
Vidíme, že rozptyl je asi o 83 % větší než průměr. To znamená, že bychom měli upravit standardní chyby vynásobením 1,35, což je odmocnina z 1,83.
R může tento výpočet provést za nás, pokud použijeme quasipoissonrodinu:
> 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
Alternativním přístupem je fitovat Poissonův model a použít therobust nebo sendvičový odhad standardních chyb. To obvykle dává výsledky velmi podobné nadsazenému Poissonovu modelu.
 4.A.2 Negativní binomická regrese
Nyní fitujeme negativní binomický model se stejnými prediktory. k tomu potřebujeme funkci glm.nb() v balíčkuMASS.
> 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 theta je přesnost multiplikativního náhodného efektu a odpovídá 1/σ2 v thenotes. Odhad odpovídá odhadovanému rozptylu0,44 a je vysoce významný.
Pro testování významnosti tohoto parametru lze uvažovat o výpočtu dvojnásobku rozdílu logaritmických pravděpodobností mezi tímto modelem a Poissonovým modelem, 180,2, a považovat jej za chí-kvadrát s jedním d.f. Obvyklá asymptotika však neplatí, protože nulová hypotéza je na hranici prostoru parametrů.
Existují práce, které ukazují, že lepší aproximací je považovat statistiku za směs 50:50 nuly a chí-kvadrátu s jedním d.f. Alternativou je považovat statistiku za chí-kvadrát s jedním d.f., což poskytuje konzervativní test. V každém případě máme zdrcující důkaz nadměrného rozptylu.
Pro testování hypotéz o regresních koeficientech můžeme použít buď Waldovy testy, nebo testy poměru pravděpodobnosti, které jsou možné, protože jsme přijali úplné distribuční předpoklady.
Pro glm existuje také rodina negative.binomial a tu lze použít za předpokladu, že jsou dány parametrytheta. (To vychází z výsledku, že záporný binomický je v rodině glm pro pevný rozptyl). Zde je rychlá kontrola:
> 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)

Můžeme také vypočítat kvantily pomocíqgamma(p, shape, rate = 1, scale = 1/rate, lower.tail = TRUE). Konkrétně, když má náhodný efekt rozptyl v, jsou kvartily 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
Biochemici na Q1 rozdělení nepozorované heterogenity publikují o 49 % méně prací, než se očekává z jejich pozorovaných charakteristik, zatímco ti na mediánu publikují o 14 % méně a ti na Q3 o 33 % více, než se očekává.
Srovnání odhadů a standardních chyb
Srovnejme odhady parametrů a standardní chyby podle Poissonova modelu, modelu s nadměrným rozptylem Poisson a negativního binomického modelu:
> 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
Odhady podle negativního binomického modelu se příliš neliší od odhadů založených na Poissonově modelu a obě sady by vedly ke stejným závěrům. 
Podíváme-li se na standardní chyby, vidíme, že oba přístupy k nadměrnému rozptylu vedou k velmi podobným odhadům standardních chyb,a že běžná Poissonova regrese standardní chyby podhodnocuje.
Dobrá shoda
Dobrou shodu záporného binomického modelu můžeme posoudit pomocí odchylky
> deviance(mnbg) 1004.281
Záporný binomický model sedí lépe než Poissonův, ale stále má odchylku nad pětiprocentní kritickou hodnotou980,25.
Funkce rozptylu
Model s nadměrným rozptylem Poisson a záporný binomický model mají různé funkce rozptylu. Jedním ze způsobů, jak ověřit, který z nich může být vhodnější, je vytvořit skupiny na základě lineárního prediktoru, vypočítat průměr a rozptyl pro každou skupinu a nakonec vykreslit vztah mezi průměrem a rozptylem.
Tady jsou skupiny založené na lineárním prediktoru záporného binomu, vytvořené pomocí cut s přerušením na 5(5)95 percentilech, aby vzniklo dvacet přibližně stejně velkých skupin. Všimněte si, že predict() počítá lineární prediktor,pokud není uvedeno jinak. Chcete-li predikovat v měřítku odpovědi, přidejte možnost 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 

Kraf vykresluje průměr versus rozptyl a překrývá křivky odpovídající nadměrně rozptýlenému Poissonovu modelu, kde je rozptyl φμ, a zápornému binomickému modelu, kde je rozptyl μ(1+μσ2).
Poissonova rozptylová funkce odvádí docela dobrou práci pro většinu dat, ale nedokáže zachytit vysoké rozptyly nejproduktivnějších vědců. Záporná binomická funkce rozptylu se příliš neliší, ale protože je kvadratická, může rychleji narůstat a odvádí lepší práci na vysokém konci. Dospěli jsme k závěru, že záporný binomický model poskytuje lepší popis dat než příliš rozptýlený Poissonův model.
 4.A.3 Poissonovy modely s nulovým rozptylem
 Častým jevem u počítaných dat je nadbytek nul ve srovnání s tím, co se očekává u Poissonova modelu. To je u našich dat skutečně problém: 
> zobs zpoi c(obs=mean(zobs), poi=mean(zpoi)) obs poi 0.3005464 0.2092071 
Vidíme, že 30,0 % vědců ve vzorku nepublikovalo v posledních třech letech svého doktorátu žádný článek, ale Poissonův model předpovídá, že pouze 20,9 % nebude mít žádné publikace. Je zřejmé, že model podhodnocuje pravděpodobnost nulových počtů.
Jedním ze způsobů, jak modelovat tento typ situace, je předpokládat, že data pocházejí ze směsi dvou populací, jedné, kde je počet vždy nulový, a druhé, kde má počet Poissonovo rozdělení se střední hodnotou μ . V tomto modelu mohou nulové počty pocházet z obou populací, zatímco kladné počty pocházejí pouze z druhé populace.
V souvislosti s publikacemi doktorandů biochemie si můžeme představit, že někteří měli na mysli zaměstnání, kde by publikace nebyly důležité, zatímco jiní usilovali o akademická místa, kde se očekával rekordní počet publikací. Členové první skupiny by publikovali nula článků, zatímco členové druhé skupiny by publikovali 0,1,2,..., přičemž lze předpokládat, že tento počet má Poissonovo rozdělení.
Rozdělení výsledku pak lze modelovat pomocí dvou parametrů, π pravděpodobnost "vždy nula" a μ, průměrný počet publikací pro ty, kteří nejsou ve skupině "vždy nula". Přirozeným způsobem, jak zavést kovariáty, je modelovat logit pravděpodobnosti π vždy nula a logaritmus průměrného počtu μ pro ty, kteří nejsou ve skupině vždy nula.
Tento typ modelu lze fitovat v R pomocí funkce zeroinfl() v balíčku pscl. Vzorec modelu lze zadat jako obvykle, pokud mají být do obou rovnic zahrnuty stejné proměnné. V opačném případě lze zadat dvě sady prediktorůoddělených svislým sloupcem, pro více informací zadejte ?zeroinfl.
Tady je model s nulovou inflací se všemi kovariátami v obou rovnicích:
> 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
Podíváme-li se na inflační rovnici, vidíme, že jediným významným prediktorem pro zařazení do třídy "vždy nula" je počet článků publikovaných mentorem, přičemž každý článek mentora je spojen s 12,6% nižší pravděpodobností, že nikdy nepublikuje.
Podíváme-li se na rovnici pro průměrný počet nebo články mezi těmi, kteří nejsou ve třídě "vždy nula", zjistíme významné znevýhodnění žen a vědců s dětmi mladšími pěti let a velký pozitivní vliv počtu publikací mentora, přičemž každý článek je spojen s 1,8% zvýšením očekávaného počtu publikací.
Pro ověření, že model řeší problém přebytku nul, predikujeme π a μ a vypočítáme kombinovanou pravděpodobnost nepublikování. Pro jejich získání existují ve funkci predict()volby nazvané "zero" a "count". Existuje také volba "prob" pro výpočetpředpovídané hustoty, ale ta je nadbytečná, protože chceme pouzepravděpodobnost nuly.
> pr mu zip mean(zip) 0.2985679
Model tedy řeší problém přebytku nul a předpovídá, že 29,9 % biochemiků nebude publikovat žádné články, což je mnohem blíže pozorované hodnotě 30,0 %.
Existuje také model s nulovým infiltrátem negativního binomu, který používázáporný binom pro počet ve třídě "ne vždy nula". Tento model lze fitovat pomocí zeroinfl() s parametrem dist="negbin". Alternativní vazby pro nafouknutou rovnici zahrnují probit,který lze zadat pomocí link="probit". 
Srovnání modelu s AIC
Jak to tak bývá, pro tato data řeší problém i záporný binomický model! Zde je pravděpodobnost nulových článků v záporném binomickém modelu. Vycházíme z prvních principů, ale bylo by možné použít i vestavěnou hustotu záporného binomu
> munb theta znb # also dnbinom(0, mu=munb, size=theta)> mean(znb) 0.3035957
Model předpovídá, že 30,4 % biochemiků by v posledních třech letech doktorského studia nepublikovalo žádný článek, což je velmi blízko pozorované hodnotě 30,0 %.
Pro volbu mezi modelem se záporným binomem a modelem s nafouknutou nulou se musíme uchýlit k dalším kritériím. Velmi jednoduchým způsobem, jak porovnat modely s různým počtem parametrů, je výpočet Akaikeho informačního kritéria (AIC), které definujeme jako
AIC = -2logL + 2p
kde p je počet parametrů v modelu. První člen je v podstatě odchylka a druhý penalizace za počet parametrů a lze jej vypočítat pro většinu modelů pomocí funkce AIC. Získávám ho také "ručně", abychom si mohli výpočet ověřit. Pro naše data
> 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
 Pro tento soubor dat je negativní binomický model jasným vítězem z hlediska úspornosti a dobré shody. 
Další diagnostická kritéria, na která bychom se mohli podívat, jsou mezní rozdělení předpovídaných a pozorovaných počtů a funkce rozptylu. 
Nulově zkrácené a překážkové modely
Další modely, kterými jsme se nezabývali, jsou Poissonův a záporný binomický model s nulovým zkrácením, určené pro data, která neobsahují nuly. Běžným příkladem je délka pobytu v nemocnici, která činí nejméně jeden den. Rozumným přístupem je použít Poissonův nebo záporný binomický model, který vylučuje nulu a ostatní pravděpodobnosti přeškáluje tak, aby jejich součet byl roven jedné. Při interpretaci těchto modelů je třeba být opatrný, protože μ není očekávaný výsledek, ale střední hodnota základního rozdělení, které zahrnuje nuly. Tyto modely jsou implementovány ve funkci vglm() v balíčku VGAM s použitím rodin pospoisson a posnegbinomial.
Alternativním přístupem k nadbytku (nebo nedostatku) nul je použití dvoustupňového procesu s logitovým modelem pro rozlišení nulových a kladných počtů a poté s Poissonovým nebo záporným binomickým modelem s ořezanou nulou pro kladné počty. V našem příkladu bychom mohli použít logitový model pro rozlišení těch, kteří publikují, od těch, kteří nepublikují, a pak zkrácený Poissonův nebo negativní binomický model pro počet článků těch, kteří publikují alespoň jeden. Tyto modely se často nazývají překážkové modely. Lze je vypočítat v R pomocí samostatných logitových a nulových ořezaných Poissonových nebo záporných binomických modelů a jednoduše sečíst logaritmické pravděpodobnosti nebo použít funkci hurdle() v balíčku pscl.
Při porovnání překážkových a nulově zkrácených modelů se mi zdá, že rozlišení mezi nulou a jedničkou nebo více je jasnější u překážkových modelů, ale interpretace průměru je jasnější u nulově zkrácených modelů.
V následujícím článku JSS je užitečná diskuse o všech těchto modelech: Regresní modely pro početní data v jazyce R.
.