GR:n verkkosivut

4.A Mallit ylihajotetuille laskentatiedoille

Käytämme Longin (1990) tietoja tohtorintutkinnon suorittaneiden biokemistien tuottamien julkaisujen määrästä havainnollistaaksemme Poissonin, ylihajotetun Poissonin, negatiivisen binomiaalin ja nollahajotetun Poissonin mallien soveltamista.

Aineiston muuttujat ovat

art artikkelit kolmen viimeisen vuoden aikana väitöskirjan tekemisestä.

koodattu yksi naisille

mar koodattu yksi jos naimisissa. kid5 alle kuusivuotiaiden lasten lukumäärä phd prestiisi Ph.D. program

Muuttuja	Kuvaus
fem
ment	artikkelit mentorin toimesta viimeisten kolmen vuoden aikana

Tämä aineisto on analysoitu myös Longin ja Freesen (2001) toimesta, ja se on saatavissa Statan verkkosivuilta:

> 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 

Artikkeleiden lukumäärän keskiarvoksi saatiin 1,69 artikkelia, ja varianssi on 3,71, eli se on hiukan yli kaksi kertaa keskiarvo. Aineisto on ylihajonnut, mutta emme tietenkään ole vielä huomioineet mitään kovariaatteja.

A Poisson-malli

Sovitetaan Longin ja Freesen(2001) käyttämä malli, yksinkertainen additiivinen malli, jossa käytetään kaikkia viittä ennustetta.

> 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
ja poikkeama ja Pearsonin khiin neliö ovat molemmat 1600:n luokkaa.
 4.A.1. Extra-Poisson-vaihtelu
 Oletamme nyt, että varianssi on verrannollinen eikä yhtä suuri kuin keskiarvo, ja estimoimme skaalausparametrin φ jakamalla Pearsonin khiin neliö sen d.f:llä:
> phi round(c(phi,sqrt(phi)),4) 1.8290 1.3524
Havaitsemme, että varianssi on noin 83 % suurempi kuin keskiarvo. Tämä tarkoittaa, että meidän pitäisi mukauttaa keskivirheet kertomalla ne 1,35:llä, joka on 1,83:n neliöjuuri.
R voi tehdä tämän laskutoimituksen puolestamme, jos käytämme quasipoissonperhettä:
> 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
Vaihtoehtoinen lähestymistapa on sovittaa Poisson-malli ja käyttää keskivirheiden therobust- tai sandwich-estimaattoria. Tällöin saadaan yleensä hyvin samankaltaisia tuloksia kuin ylihajotetulla Poisson-mallilla.
 4.A.2 Negatiivinen binomiaaliregressio
Sovitamme nyt negatiivisen binomiaalimallin samoilla ennustajilla.tähän tarvitsemme MASS-paketin glm.nb()-funktiotaMASS.
> 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:n theta tarkkuus on multiplikatiivisen satunnaisvaikutuksen tarkkuus, ja se vastaa tällöin 1/σ2:ta. Estimaatti vastaa estimoitua varianssia0.44 ja on erittäin merkitsevä.
Testataksesi tämän parametrin merkitsevyyttä voit ajatella laskevasi tämän mallin ja Poisson-mallin log-likelihoodien erotuksen kaksinkertaisena, 180.2, ja käsitellä sitä khiin neliönä, jossa on yksi d.f.. Tavanomainen asymptotiikka ei kuitenkaan päde, koska nollahypoteesi on parametriavaruuden rajalla.
Joitakin töitä on tehty, jotka osoittavat, että parempi approksimaatio on käsitellä tilastoa nollan ja khiin neliön 50:50 sekoituksena, jossa on yksi d.f. Vaihtoehtoisesti tilastoa käsittelemällä sitä khiin neliönä, jossa on yksi d.f., saadaan konservatiivinen testi. Kummassakin tapauksessa meillä on ylivoimainen todiste ylidispersiosta.
Regressiokertoimia koskevien hypoteesien testaamiseen voimme käyttää joko Wald-testejä tai likelihood ratio -testejä, jotka ovat mahdollisia, koska olemme tehneet täydelliset jakaumaoletukset.
Lisäksi glm:lle on olemassa negative.binomial-perhe glm:lle, ja tätä voidaan käyttää edellyttäen, että parametrittheta on annettu. (Tämä perustuu tulokseen, jonka mukaan negatiivinen binomi on glm-perheessä kiinteälle varianssille). Tässä on nopea tarkistus:
> 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)

Voidaan myös laskea kvantiileja käyttämälläqgamma(p, shape, rate = 1, scale = 1/rate, lower.tail = TRUE). Erityisesti kun satunnaisvaikutuksen varianssi on vkvarttiilit ovat 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
Havaitsemattoman heterogeenisuuden jakauman Q1:ssä olevat biokemistit julkaisevat 49 % vähemmän artikkeleita kuin heidän havaittujen ominaisuuksiensa perusteella odotetaan, kun taas mediaanissa olevat julkaisevat 14 % vähemmän ja Q3:ssa olevat 33 % enemmän kuin odotetaan.
Estimaattien ja keskivirheiden vertailu
Vertaillaan parametriestimaatteja ja keskivirheitä Poissonin, ylihajotetun Poissonin ja negatiivisen binomiaalimallin mukaan:
> 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
Negatiivisen binomiaalimallin estimaatit eivät poikkea juurikaan Poissonin malliin perustuvista estimaateista, ja molemmat joukot johtaisivat samoihin päätelmiin. 
Katsottaessa keskivirheitä havaitaan, että molemmat lähestymistavat ylihajontaan johtavat hyvin samankaltaisiin estimoituihin keskivirheisiin,ja että tavallinen Poisson-regressio aliarvioi keskivirheet.
Soveltuvuuden hyvyys
Voidaan arvioida negatiivisen binomimallin soveltuvuutta käyttämällä poikkeamaa
> deviance(mnbg) 1004.281
Negatiivinen binomimalli sopii paremmin kuin Poisson-malli, mutta sen poikkeama on silti yli viiden prosentin kriittisen arvon980,25. 
Varianssifunktio
Ylihajonnalla varustetuilla Poisson- ja negatiivisella binomimallilla on erilaiset varianssifunktiot. Yksi tapa tarkistaa, kumpi niistä on sopivampi, on luoda lineaariseen ennustajaan perustuvia ryhmiä, laskea kullekin ryhmälle keskiarvo ja varianssi ja lopuksi piirtää keskiarvon ja varianssin suhde.
Tässä on negatiiviseen binomiseen lineaariseen ennustajaan perustuvia ryhmiä, jotka on luotu cut-ohjelmalla, jossa on taukoja 5(5)95-prosenttiluvun kohdalla, jotta saadaan aikaan kaksikymmentä suunnilleen samankokoista ryhmää. Huomaa, että predict() laskee lineaarisen ennustajan, ellei toisin mainita. Jos haluat ennustaa vastauksen asteikolla, lisää vaihtoehto 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 

Kuvaajassa esitetään keskiarvo suhteessa varianssiin ja asetetaan päällekkäin käyrät, jotka vastaavat ylihajotettua Poisson-mallia, jossa varianssi on φμ, ja negatiivista binomimallia, jossa varianssi on μ(1+μσ2).
Poisson-varianssifunktio tekee melko hyvää työtä suurimmalle osalle aineistosta, mutta ei onnistu kuvaamaan kaikkein tuottavimpien tiedemiesten korkeita variansseja. Negatiivinen binomivarianssifunktio ei ole kovin erilainen, mutta koska se on kvadraattinen, se voi nousta nopeammin ja tekee parempaa työtä suuressa ääripäässä. Päättelemme, että negatiivinen binomimalli kuvaa aineistoa paremmin kuin ylihajonnut Poisson-malli. 
 4.A.3 Nollatäytteiset Poisson-mallit
 Usein laskentatiedoissa esiintyy nollien ylimäärää verrattuna siihen, mitä Poisson-mallissa odotetaan. Tämä on itse asiassa ongelma aineistossamme: 
> zobs zpoi c(obs=mean(zobs), poi=mean(zpoi)) obs poi 0.3005464 0.2092071 
Me näemme, että 30,0 % otokseen kuuluvista tutkijoista ei julkaissut yhtään artikkelia väitöskirjansa kolmen viimeisen vuoden aikana, mutta Poisson-malli ennustaa, että vain 20,9 %:lla ei olisi yhtään julkaisua. Malli selvästi aliarvioi nollalukujen todennäköisyyttä.
Yksi tapa mallintaa tämäntyyppisiä tilanteita on olettaa, että aineisto on peräisin kahden populaation sekoituksesta, joista toisessa lukumäärät ovat aina nolla ja toisessa lukumäärillä on Poisson-jakauma, jonka keskiarvo μ. Tässä mallissa nollaluvut voivat tulla kummastakin populaatiosta, kun taas positiiviset luvut tulevat vain jälkimmäisestä.
Voidaan kuvitella, että tohtoriksi väitelleiden biokemistien julkaisujen yhteydessä joillakin oli mielessä työpaikat, joissa julkaisuilla ei olisi merkitystä, kun taas toiset tavoittelivat akateemisia työtehtäviä, joissa odotettiin julkaisujen määrää. Ensimmäisen ryhmän jäsenet julkaisisivat nolla artikkelia, kun taas toisen ryhmän jäsenet julkaisisivat 0,1,2,..., lukumäärän, jolla voidaan olettaa olevan Poisson-jakauma.
Tuloksen jakauma voidaan tällöin mallintaa kahdella parametrilla, π todennäköisyydellä "aina nolla" ja μ julkaisujen keskimääräisellä lukumäärällä niille, jotka eivät kuulu "aina nolla"-ryhmään. Luonnollinen tapa ottaa käyttöön kovariaatit on mallintaa "aina nolla"-luokan todennäköisyyden π logit ja "aina nolla"-luokkaan kuulumattomien julkaisujen keskiarvon μ logit.
Tämmöinen malli voidaan sovittaa R:ssä käyttämällä pscl-paketin zeroinfl()funktiota. Mallin kaava voidaan määritellä tavalliseen tapaan, jos samat muuttujat halutaan sisällyttää molempiin yhtälöihin. Muussa tapauksessa voidaan antaa kaksi ennustejoukkoa, jotka on erotettu pystypalkilla, kirjoita ?zeroinfl saadaksesi lisätietoja.
Tässä on nollaan puhallettu malli, jossa kaikki muuttujat ovat mukana molemmissa yhtälöissä:
> 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
Katsomalla puhallettua yhtälöä huomaamme, että ainoa merkitsevä ennustaja "aina nollaa" -luokkaan kuulumiselle on mentorin julkaisemien artikkeleiden määrä, ja jokainen mentorin julkaisema artikkeli liittyy 12,6 % alhaisempaan todennäköisyyteen olla koskaan julkaisematta.
Katsottaessa yhtälöä, joka koskee keskimääräistä artikkelien lukumäärää niiden joukossa, jotka eivät kuulu "aina nolla"-luokkaan, havaitaan, että naisilla ja tutkijoilla, joilla on alle viiden vuoden ikäisiä lapsia, on merkittäviä haittoja ja että mentorin julkaisujen lukumäärällä on suuri positiivinen vaikutus, sillä jokainen artikkeli liittyy 1,8 %:n lisäykseen odotettavissa olevaan julkaisujen lukumäärään.
Varmistuaksemme siitä, että malli ratkaisee liiallisista nollakertoimista aiheutuvan ongelman ennustamme π:n ja μ:n ja laskemme yhteenlasketut todennäköisyydet jäädä ilman julkaisuja. Näiden saamiseksi predict()funktiossa on vaihtoehdot "zero" ja "count". On myös vaihtoehto "prob", jolla lasketaan ennustettu tiheys, mutta tämä on turhaa, koska haluamme vain nollan todennäköisyyden.
> pr mu zip mean(zip) 0.2985679
Malli ratkaisee siis liiallisen nollan ongelman ja ennustaa, että 29,9 % biokemisteistä ei julkaise yhtään artikkelia, mikä on paljon lähempänä havaittua arvoa, joka on 30,0 %.
Mallilla on myös nollaan puhallutettu negatiivinen binomimalli, joka käyttää negatiivista binomia laskennassa luokassa 'ei aina nolla'. Tämä malli voidaan sovittaa käyttämällä zeroinfl() parametria dist="negbin". Vaihtoehtoisia linkkejä inflate-yhtälölle ovat probit,joka voidaan määrittää käyttämällä link="probit". 
Mallien vertailu AIC:n avulla
Sattumoisin myös negatiivinen binomi ratkaisee ongelman tässä aineistossa! Tässä on nolla-artikkelien todennäköisyys negatiivisella binomilla. Lähdemme liikkeelle ensimmäisistä periaatteista, mutta voitaisiin myös käyttää sisäänrakennettua negatiivista binomitiheyttä
> munb theta znb # also dnbinom(0, mu=munb, size=theta)> mean(znb) 0.3035957
Malli ennustaa, että 30,4 % biokemisteistä ei julkaisisi yhtään artikkelia tohtorintutkintonsa kolmen viimeisen vuoden aikana, mikä on hyvin lähellä havaittua arvoa 30,0 %.
Voidaksemme tehdä valinnan negatiivisen binomi- ja nollatäydennysmallien väliltä joudumme turvautumaan muihin kriteereihin. Erittäin yksinkertainen tapa vertailla malleja, joissa on eri määrä parametreja, on laskea Akaiken informaatiokriteeri (AIC, Akaike's Information Criterion), joka määritellään seuraavasti
AIC = -2logL + 2p
jossa p on mallin parametrien määrä. Ensimmäinen termi on pohjimmiltaan poikkeama ja toinen parametrien lukumäärää koskeva rangaistus, ja se voidaan laskea useimmille malleille AIC-funktion avulla. Saan sen myös "käsin", jotta voimme tarkistaa laskennan. Aineistomme osalta
> 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
 Tämän aineiston osalta negatiivinen binomimalli on selvä voittaja yksinkertaisuuden ja sopivuuden suhteen. 
Muut diagnostiset kriteerit, joita voisimme tarkastella, ovat ennustettujen ja havaittujen lukumäärien marginaalijakaumat ja varianssifunktiot. 
Nollakatkoviivoitetut ja hurdle-mallit
Muut mallit, joita emme ole käsitelleet, ovat nollakatkoviivoitetut Poisson- ja negatiivinen binomimalli, jotka on suunniteltu aineistoille, jotka eivät sisällä nollia. Yleinen esimerkki on sairaalassaoloaika, joka on vähintään yksi päivä. Järkevä lähestymistapa on sovittaa Poisson- tai negatiivinen binomimalli, joka sulkee pois nollan ja skaalaa muut todennäköisyydet siten, että niiden summa on yksi. Näiden mallien tulkinnassa on oltava varovainen, koska μ ei ole odotettu lopputulos vaan nollat sisältävän jakauman keskiarvo. Nämä mallit on toteutettu VGAM-paketin vglm()-funktiossa käyttäen perheitä pospoisson ja posnegbinomial.
Vaihtoehtoinen lähestymistapa nollien ylimäärään (tai puutteeseen) on käyttää kaksivaiheista prosessia, jossa käytetään logit-mallia nolla- ja positiivisten lukumäärien erottamiseksi toisistaan ja sen jälkeen nollaan karsittua Poisson- tai negatiivista binomimallia positiivisille lukumäärille. Esimerkissämme voisimme käyttää logit-mallia erottaaksemme toisistaan ne, jotka julkaisevat, ja ne, jotka eivät julkaise, ja sitten typistettyä Poisson- tai negatiivista binomimallia niiden artikkelien lukumäärälle, jotka julkaisevat vähintään yhden. Näitä malleja kutsutaan usein esteellisyysmalleiksi. Ne voidaan sovittaa R:ssä käyttämällä erillisiä logit- ja nollakarsittujaPoisson- tai negatiivisia binomimalleja ja yksinkertaisesti lisäämällä log-likelihoodit tai käyttämällä pscl-paketin hurdle()-funktiota.
Vertaillessani hurdle- ja nollatyhjennettyjä malleja pidän eroa nollan ja yhden tai useamman välillä selkeämpänä hurdle-malleilla, mutta keskiarvon tulkinta on selkeämpi nollatyhjennetyillä malleilla.
Oheisessa JSS-artikkelissa on hyödyllinen keskustelu kaikista näistä malleista: Regression Models for Count Data in R.