GR’s hjemmeside

4.A Modeller for over-spredte tælledata

Vi bruger data fra Long (1990) om antallet af publikationer produceret af ph.d.-biokemikere til at illustrere anvendelsen af Poisson-, over-spredte Poisson-, negative binomial- og nul-inflaterede Poisson-modeller.

Variablerne i datasættet er

art

kodet en hvis gift

phd prestige af ph.d.-grad.D. program

Variabel	Beskrivelse
Artikler i de sidste tre år af ph.d.-studiet.
fem	kodet en for kvinder
mar
kid5	antal børn under seks år
ment	artikler af mentor i de sidste tre år

Disse data er også blevet analyseret af Long og Freese (2001), og er tilgængelige fra Stata-webstedet:

> 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 

Middelværdien af antallet af artikler er 1,69 og variansen er 3,71, lidt mere end det dobbelte af middelværdien. Dataene er overspredte, men vi har naturligvis ikke overvejet nogen kovariater endnu.

En Poisson-model

Lad os tilpasse den model, som Long og Freese(2001) har anvendt, en simpel additiv model, der anvender alle fem prædiktorer.

> 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
og afvigelsen og Pearsons chi-kvadrat er begge i 1600'erne.
 4.A.1. Ekstra-Poisson-variation
Vi antager nu, at variansen er proportional i stedet for lig med middelværdien, og estimerer skala-parameteren φ ved at dividere Pearsons chi-kvadrat med dens d.f.:
> phi round(c(phi,sqrt(phi)),4) 1.8290 1.3524
Vi ser, at variansen er ca. 83% større end middelværdien. Det betyder, at vi bør justere standardfejlene ved at gange med 1,35, kvadratroden af 1,83.
R kan foretage denne beregning for os, hvis vi bruger quasipoissonfamilien:
> 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
En alternativ tilgang er at tilpasse en Poisson-model og bruge therobust eller sandwich estimator af standardfejlene. Dette giver normalt resultater, der minder meget om den overspredte Poisson-model.
 4.A.2 Negativ binomialregression
Vi tilpasser nu en negativ binomialmodel med de samme prædiktorer.For at gøre dette har vi brug for glm.nb()-funktionen iMASS-pakken.
> 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 er præcisionen af den multiplikative binomiale effekt, og svarer til 1/σ2 i thenotes. Estimatet svarer til en estimeret varians på0,44 og er meget signifikant.
For at teste denne parameters signifikans kan man tænke på at beregneden dobbelte forskel i log-likelihoods mellem denne model og Poisson-modellen, 180,2, og behandle den som et chi-squared med en d.f. Den sædvanlige asymptotik gælder dog ikke, fordi nulhypotesen ligger på en grænse af parameterrummet.
Der er noget arbejde, der viser, at en bedre tilnærmelse er at behandle statistikken som en 50:50-blanding af nul og en chi-squared med en d.f. Alternativt giver en konservativ test, hvis man behandler statistikken som en chi-squared meden d.f.. Uanset hvad har vi overvældende beviser for overspredning.
For at teste hypoteser om regressionskoefficienterne kan vi bruge enten Wald-test eller likelihood ratio-test, hvilket er muligt, fordi vi har gjort fulde fordelingsantagelser.
Der er også en negative.binomial-familie for glm, og denne kan bruges, forudsat at parametrenetheta er givet. (Dette er baseret på det resultat, at det negative binomium er i glm-familien for fast varians). Her er en hurtig kontrol:
> 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)

Vi kan også beregne kvantiser ved hjælp afqgamma(p, shape, rate = 1, scale = 1/rate, lower.tail = TRUE). Især når den tilfældige effekt har varians v, er kvartilerne 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
Biokemikere på Q1 i fordelingen af uobserveret heterogenitet udgiver 49% færre artikler end forventet ud fra deres observerede karakteristika, mens dem på medianen udgiver 14% færre og dem på Q3 udgiver 33% flere end forventet.
Sammenligning af estimater og standardafvigelser
Lad os sammenligne parameterestimater og standardafvigelser under Poisson-modellen, den overspredte Poisson-model og den negative binomialmodel:
> 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
De negative binomiale estimater er ikke meget forskellige fra dem, der er baseret på Poisson-modellen, og begge sæt ville føre til de samme konklusioner. 
Hvis vi ser på standardfejlene, ser vi, at begge tilgange til overspredning fører til meget ensartede estimerede standardfejl, og at den almindelige Poisson-regression undervurderer standardfejlene.
Godhed af tilpasning
Vi kan vurdere godheden af tilpasning af den negative binomialmodel ved hjælp af afvigelsen
> deviance(mnbg) 1004.281
Den negative binomialmodel passer bedre end Poisson-modellen, men har stadig en afvigelse over den kritiske værdi på fem procent på980,25.
Variansfunktionen
Den overspredte Poisson- og negative binomialmodel har forskellige variansfunktioner. En måde at kontrollere, hvilken der kan være mest hensigtsmæssig, er at oprette grupper baseret på den lineære prædiktor, beregne middelværdien og variansen for hver gruppe og endelig plotte forholdet mellem middelværdi og varians.
Her er grupper baseret på den negative binomiale lineære prædiktor, oprettet ved hjælp af cut med afbrydelser ved 5(5)95-percentilernefor at producere tyve grupper af omtrent samme størrelse. Bemærk, at predict() beregner en lineær prædiktor, medmindre andet er angivet. Hvis du vil forudsige i svarets skala, skal du tilføje indstillingen 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 

Diagrammet viser middelværdien i forhold til variansen og overlejrer de kurver, der svarer til den overspredte Poisson-model, hvor variansen er φμ, og den negative binomialmodel, hvor variansen er μ(1+μσ2).
Poisson-variansfunktionen gør et ret godt stykke arbejde for størstedelen af dataene, men formår ikke at indfange de høje varianser for de mest produktive forskere. Den negative binomiale variansfunktion er ikke meget anderledes, men da den er kvadratisk, kan den stige hurtigere og gør et bedre stykke arbejde i den høje ende. Vi konkluderer, at den negative binomialmodel giver en bedre beskrivelse af dataene end den overdrevent spredte Poisson-model.
 4.A.3 Nul-inflaterede Poisson-modeller
 En hyppig forekomst med tælledata er et overskud af nuller i forhold til, hvad der forventes under en Poisson-model. Dette er faktisk et problem med vores data: 
> zobs zpoi c(obs=mean(zobs), poi=mean(zpoi)) obs poi 0.3005464 0.2092071 
Vi ser, at 30,0 % af forskerne i stikprøven ikke udgav nogen artikler i de sidste tre år efter deres ph.d., men Poisson-modellen forudsiger, at kun 20,9 % ikke ville have nogen publikationer. Modellen undervurderer tydeligvis sandsynligheden for nultællinger.
En måde at modellere denne type situation på er at antage, at dataene kommer fra en blanding af to populationer, en hvor tællingerne altid er nul, og en anden hvor tællingerne har en Poisson-fordeling med middelværdi μ. I denne model kan nultællinger komme fra begge populationer, mens positive tællinger kun kommer fra den anden population.
I forbindelse med publikationer fra ph.d.-biokemikere kan vi forestille os, at nogle havde i tankerne job, hvor publikationer ikke ville være vigtige, mens andre sigtede mod akademiske job, hvor man forventede en række publikationer. Medlemmer af den første gruppe ville publicere nul artikler, mens medlemmer af den anden gruppe ville publicere 0,1,2,..., et antal, der kan antages at have en Poisson-fordeling.
Fordelingen af resultatet kan derefter modelleres i form af to parametre, π sandsynligheden for "altid nul", og μ, det gennemsnitlige antal publikationer for dem, der ikke er i "altid nul"-gruppen. En naturlig måde at indføre kovariater på er at modellere logit af sandsynligheden π for altid nul og logit af middelværdien μ for dem, der ikke er i gruppen af altid nul.
Denne type model kan tilpasses i R ved hjælp af zeroinfl()funktionen i pscl-pakken. Modelformlen kan specificeres som sædvanlig, hvis de samme variabler skal indgå i begge ligninger. Ellers kan man angive to sæt prædiktorer adskilt af en lodret streg, skriv ?zeroinfl for at få mere at vide.
Her er en nul-infleret model med alle kovariater i begge ligninger:
> 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
Kigger vi på den inflaterede ligning, ser vi, at den eneste signifikante prædiktor for at være i klassen "altid nul" er antallet af artikler udgivet af mentoren, idet hver artikel af mentoren er forbundet med 12,6 % lavere odds for aldrig at udgive.
Kigger vi på ligningen for det gennemsnitlige antal artikler blandt dem, der ikke er i klassen "altid nul", finder vi signifikante ulemper for kvinder og forskere med børn under fem år og en stor positiv effekt af antallet af mentorens publikationer, idet hver artikel er forbundet med en stigning på 1,8 % i det forventede antal publikationer.
For at verificere, at modellen løser problemet med overskydende nuller, forudsiger vi π og μ og beregner den samlede sandsynlighed for ikke at publicere. Der er muligheder i predict()funktionen kaldet "zero" og "count" for at opnå disse. Der er også en mulighed "prob" for at beregne den forudsagte tæthed, men det er overkill, da vi kun ønsker sandsynligheden for nul.
> pr mu zip mean(zip) 0.2985679
Så modellen løser problemet med overskydende nuller og forudsiger, at 29,9 % af biokemikerne ikke vil udgive nogen artikler, hvilket er meget tættere på den observerede værdi på 30,0 %.
Der findes også en nul-inflated negativ binomialmodel, som bruger en negativ binomial for tallet i klassen "ikke altid nul". Denne model kan tilpasses ved hjælp af zeroinfl() med parameteren dist="negbin". Alternative forbindelser til den oppustede ligning omfatter probit,som kan specificeres ved hjælp af link="probit". 
Modelsammenligning med AIC
For disse data løser den negative binomialmodel tilfældigvis også problemet! Her er sandsynligheden for nul artikler i det negative binomial. Vi går ud fra første principper, men man kunne også bruge den indbyggede negative binomialdensitet
> munb theta znb # also dnbinom(0, mu=munb, size=theta)> mean(znb) 0.3035957
Modellen forudsiger, at 30,4% af biokemikerne ikke ville publicere nogen artikler i de sidste tre år af deres ph.d., hvilket er meget tæt på den observerede værdi på 30,0%.
For at vælge mellem den negative binomialmodel og den nulopblæste model må vi ty til andre kriterier. En meget enkel måde at sammenligne modeller med forskellige antal parametre på er at beregne Akaike's Informationskriterium (AIC), som vi definerer som
AIC = -2logL + 2p
hvor p er antallet af parametre i modellen. Det første udtryk er i det væsentlige afvigelsen og det andet en straf for antallet af parametre, og det kan beregnes for de fleste modeller ved hjælp af AIC-funktionen. Jeg får det også "i hånden", så vi kan verificere beregningen. For vores 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
 For dette datasæt er den negative binomialmodel en klar vinder med hensyn til parsimony og goodness of fit. 
Andre diagnostiske kriterier, som vi kunne se på, er den marginale fordeling af forudsagte og observerede tællinger og variansfunktionerne. 
Zero-Truncated- og Hurdle-modeller
Andre modeller, som vi ikke har dækket, er den nultrunkerede Poisson- og den negative binomialmodel, der er designet til data, der ikke indeholder nuller. Et almindeligt eksempel er længden af opholdet på et hospital, som er mindst én dag. En fornuftig fremgangsmåde er at tilpasse en Poisson- eller negativ binomialmodel, der udelukker nul og omskalerer de andre sandsynligheder, så de summerer til én. Man skal være forsigtig med at fortolke disse modeller, fordi μ ikke er det forventede resultat, men middelværdien af en underliggende fordeling, der omfatter nuller. Disse modeller er implementeret i funktionen vglm() i pakken VGAM ved hjælp af familierne pospoisson og posnegbinomial.
En alternativ tilgang til overskud (eller mangel) af nuller er at anvende en to-trins proces med en logit-model til at skelne mellem nul- og positive tællinger og derefter en nul-truncated Poisson- eller negativ binomialmodel for de positive tællinger. I vores eksempel kunne vi bruge en logit-model til at skelne mellem dem, der publicerer, og dem, der ikke gør det, og derefter en truncated Poisson- eller negativ binomialmodel for antallet af artikler for dem, der publicerer mindst én. Disse modeller kaldes ofte hurdle-modeller. De kan tilpasses i R ved hjælp af de separate logit- og nul-truncatedPoisson- eller negative binomialmodeller og blot ved at lægge log-likelihoods sammen, eller ved hjælp af hurdle()-funktionen i pscl-pakken.
Ved sammenligning af hurdle- og nul-inflaterede modeller finder jeg, at sondringen mellem nul og en eller flere er tydeligere med hurdle-modeller, men fortolkningen af middelværdien er tydeligere med nul-inflaterede modeller.
Der følgende JSS-papir indeholder en nyttig diskussion af alle disse modeller: Regression Models for Count Data in R.