GR:s webbplats

4.A Modeller för överspridda räkningsdata

Vi använder data från Long (1990) om antalet publikationer som produceras av doktorander i biokemi för att illustrera tillämpningen av Poisson-, överspridda Poisson-, negativa binomial- och nollinflaterade Poisson-modeller.

Variablerna i datasetet är

kodat ett om gift

Variabel	Beskrivning
art	artiklar under de tre sista åren av doktorsexamen.
fem	kodat ett för kvinnor
mar
kid5	antal barn under sex år
phd	prestige för doktorsexamen.D. program
ment	artiklar av mentor under de senaste tre åren

Dessa data har också analyserats av Long och Freese (2001) och finns tillgängliga på Statas webbplats:

> 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 

Medeltalet för antalet artiklar är 1,69 och variansen är 3,71, lite mer än dubbelt så mycket som medelvärdet. Uppgifterna är överspridda, men vi har naturligtvis inte tagit hänsyn till några kovarianter ännu.

En Poissonmodell

Låt oss anpassa den modell som Long och Freese(2001) använde, en enkel additiv modell som använder alla fem prediktorer.

> 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
och avvikelsen och Pearsons chi-kvadrat ligger båda på 1600-talet.
 4.A.1. Extra-Poisson-variation
Vi antar nu att variansen är proportionell snarare än lika med medelvärdet och uppskattar skalparametern φ genom att dividera Pearsons chi-kvadrat med dess d.f.:
> phi round(c(phi,sqrt(phi)),4) 1.8290 1.3524
Vi ser att variansen är ungefär 83 % större än medelvärdet. Detta innebär att vi bör justera standardfelen genom att multiplicera med 1,35, kvadratroten av 1,83.
R kan göra denna beräkning åt oss om vi använder quasipoissonfamily:
> 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
Ett alternativt tillvägagångssätt är att anpassa en Poisson-modell och använda therobust- eller sandwich-estimatorn av standardfelen. Detta ger vanligtvis resultat som är mycket lika den överspridda Poissonmodellen.
 4.A.2 Negativ binomialregression
Vi anpassar nu en negativ binomialmodell med samma prediktorer.För att göra detta behöver vi glm.nb()-funktionen iMASS-paketet.
> 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 är precisionen för den multiplikativa slumpmässiga effekten, och motsvarar 1/σ2 i dåvarande noter. Uppskattningen motsvarar en uppskattad varians på0,44 och är mycket signifikant.
För att testa betydelsen av denna parameter kan man tänka sig att beräkna två gånger skillnaden i log-likelihoods mellan denna modell och Poisson-modellen, 180,2, och behandla den som en chi-squared med en d.f. Den vanliga asymptotiken gäller dock inte, eftersom nollhypotesen ligger på en gräns av parameterområdet.
Det finns en del arbete som visar att en bättre approximation är att behandla statistiken som en 50:50-blandning av noll och en chi-squared med en d.f. Alternativt ger behandlingen av statistiken som en chi-squared med en d.f. ett konservativt test. Oavsett vilket har vi överväldigande bevis för överspridning.
För att testa hypoteser om regressionskoefficienterna kan vi använda antingen Wald-test eller likelihood ratio-test, vilket är möjligt eftersom vi har gjort fullständiga fördelningsantaganden.
Det finns också en negative.binomial-familj för glm och denna kan användas under förutsättning att parametrarnatheta är givna. (Detta baseras på resultatet att det negativa binomialet ingår i glm-familjen för fast varians). Här är en snabb kontroll:
> 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 också beräkna kvantilar med hjälp avqgamma(p, shape, rate = 1, scale = 1/rate, lower.tail = TRUE). I synnerhet när den slumpmässiga effekten har varians v är kvartilerna 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
Biokemister som befinner sig på Q1 i fördelningen av icke-observerad heterogenitet publicerar 49 % färre artiklar än vad som förväntas utifrån deras observerade egenskaper, medan de som befinner sig på medianen publicerar 14 % färre och de som befinner sig på Q3 publicerar 33 % mer än vad som förväntas.
Genomförande av skattningar och standardavvikelser
Låt oss jämföra parameterskattningar och standardavvikelser enligt Poisson-, överspridda Poisson- och negativa binomialmodellerna:
> 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 negativa binomialskattningarna skiljer sig inte särskilt mycket från de som baseras på Poissonmodellen, och båda uppsättningarna skulle leda till samma slutsatser. 
Om vi tittar på standardfelen ser vi att båda tillvägagångssätten för överspridning leder till mycket likartade uppskattade standardfel, och att den vanliga Poissonregressionen underskattar standardfelen.
Goodness of Fit
Vi kan bedöma hur passande den negativa binomialmodellen är med hjälp av avvikelsen
> deviance(mnbg) 1004.281
Den negativa binomialmodellen passar bättre än Poissonmodellen, men har fortfarande en avvikelse som ligger över det femprocentiga kritiska värdet på980,25.
Variansfunktionen
Den överspridda Poissonmodellen och negativa binomialmodellen har olika variansfunktioner. Ett sätt att kontrollera vilken som kan vara lämpligast är att skapa grupper baserade på den linjära prediktoren, beräkna medelvärde och varians för varje grupp och slutligen plotta medelvärde-variansförhållandet.
Här är grupper baserade på den negativa binomiala linjära prediktoren, som skapats med hjälp av cut med avbrott vid 5(5)95-percentilerna för att skapa tjugo grupper av ungefär samma storlek. Observera att predict() beräknar en linjär prediktor, om inget annat anges. Om du vill förutsäga i svarets skala, lägg till alternativet 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 

Grafen visar medelvärdet i förhållande till variansen och överlagrar kurvorna som motsvarar den överspridda Poissonmodellen, där variansen är φμ, och den negativa binomialmodellen, där variansen är μ(1+μσ2).
Poissonvariansfunktionen gör ett ganska bra jobb för huvuddelen av data, men misslyckas med att fånga upp de höga varianserna hos de mest produktiva forskarna. Den negativa binomiala variansfunktionen är inte alltför annorlunda, men eftersom den är kvadratisk kan den stiga snabbare och gör ett bättre jobb i den höga delen. Vi drar slutsatsen att den negativa binomialmodellen ger en bättre beskrivning av data än den överspridda Poissonmodellen.
 4.A.3 Noll-inflaterade Poissonmodeller
 En vanlig företeelse med räkneuppgifter är ett överskott av nollor jämfört med vad som förväntas enligt en Poissonmodell. Detta är faktiskt ett problem med våra data: 
> zobs zpoi c(obs=mean(zobs), poi=mean(zpoi)) obs poi 0.3005464 0.2092071 
Vi ser att 30,0 % av vetenskapsmännen i urvalet inte publicerade några artiklar under de tre sista åren efter sin doktorsexamen, men Poisson-modellen förutspår att endast 20,9 % skulle ha inga publikationer. Det är uppenbart att modellen underskattar sannolikheten för nollräkningar.
Ett sätt att modellera denna typ av situation är att anta att data kommer från en blandning av två populationer, en där räkningen alltid är noll och en annan där räkningen har en Poissonfördelning med medelvärde μ. I denna modell kan nollräkningar komma från båda populationerna, medan positiva räkningar endast kommer från den andra populationen.
I samband med publikationer av disputerade biokemister kan vi föreställa oss att vissa hade i åtanke arbeten där publikationer inte skulle vara viktiga, medan andra syftade på akademiska arbeten där man förväntade sig ett antal publikationer. Medlemmarna i den första gruppen skulle publicera noll artiklar, medan medlemmarna i den andra gruppen skulle publicera 0,1,2,..., ett antal som kan antas ha en Poissonfördelning.
Fördelningen av utfallet kan då modelleras i termer av två parametrar, π sannolikheten för "alltid noll", och μ, medelantalet publikationer för dem som inte ingår i "alltid noll"-gruppen. Ett naturligt sätt att införa kovariater är att modellera logit av sannolikheten π för alltid noll och logit av medelvärdet μ för dem som inte ingår i klassen alltid noll.
Denna typ av modell kan anpassas i R med hjälp av funktionen zeroinfl() i paketet pscl. Modellformeln kan specificeras som vanligt om samma variabler ska ingå i båda ekvationerna. Annars kan man ange två uppsättningar av prediktorer som separeras av ett vertikalt streck, skriv ?zeroinfl för att få veta mer.
Här är en nollinflaterad modell med alla kovariater i båda ekvationerna:
> 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
Om vi tittar på den inflaterade ekvationen ser vi att den enda signifikanta prediktoren för att tillhöra klassen "alltid noll" är antalet artiklar som publicerats av mentorn, där varje artikel av mentorn är förknippad med 12,6 % lägre odds för att aldrig publicera.
Om vi tittar på ekvationen för det genomsnittliga antalet artiklar bland dem som inte tillhör klassen "alltid noll" finner vi betydande nackdelar för kvinnor och forskare med barn under fem år, och en stor positiv effekt av antalet publikationer av mentorn, där varje artikel är förknippad med en ökning av det förväntade antalet publikationer med 1,8 %.
För att verifiera att modellen löser problemet med överflödiga nollor förutspår vi π och μ och beräknar den kombinerade sannolikheten att aldrig publicera. Det finns alternativ i predict()funktionen som heter "zero" och "count" för att få fram dessa. Det finns också ett alternativ "prob" för att beräkna den förutspådda tätheten, men detta är överflödigt eftersom vi bara vill ha sannolikheten för noll.
> pr mu zip mean(zip) 0.2985679
Modellen löser alltså problemet med alltför många nollor och förutspår att 29,9 % av biokemisterna inte kommer att publicera några artiklar, vilket är mycket närmare det observerade värdet 30,0 %.
Det finns också en nollinflaterad negativ binomialmodell, där man använder en negativ binomial för räkningen i klassen "inte alltid noll". Denna modell kan anpassas med hjälp av zeroinfl() med parametern dist="negbin". Alternativa länkar för den uppblåsta ekvationen inkluderar probit,som kan specificeras med link="probit". 
Modelljämförelse med AIC
För dessa data löser det negativa binomialet också problemet! Här är sannolikheten för nollartiklar i det negativa binomialet. Vi utgår från första principer men man skulle också kunna använda den inbyggda negativa binomialdensiteten
> munb theta znb # also dnbinom(0, mu=munb, size=theta)> mean(znb) 0.3035957
Modellen förutsäger att 30,4 % av biokemisterna inte skulle publicera några artiklar under de tre sista åren av sin doktorsexamen, vilket ligger mycket nära det observerade värdet 30,0 %.
För att kunna välja mellan den negativa binomialmodellen och den nolluppblåsta modellen måste vi tillgripa andra kriterier. Ett mycket enkelt sätt att jämföra modeller med olika antal parametrar är att beräkna Akaikes informationskriterium (AIC), som vi definierar som
AIC = -2logL + 2p
där p är antalet parametrar i modellen. Den första termen är i huvudsak avvikelsen och den andra ett straff för antalet parametrar, och den kan beräknas för de flesta modeller med hjälp av funktionen AIC. Jag erhåller den också "för hand" så att vi kan verifiera beräkningen. För våra 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
 För detta dataset är den negativa binomialmodellen en klar vinnare när det gäller sparsamhet och passningsnoggrannhet. 
Andra diagnostiska kriterier som vi kan titta på är marginalfördelningen av förutsagda och observerade antal och variansfunktionerna. 
Zerotrunkerade modeller och Hurdle-modeller
Andra modeller som vi inte har behandlat är den nolltrunkerade Poissonmodellen och den negativa binomialmodellen, som är utformade för data som inte innehåller nollor. Ett vanligt exempel är vistelsetid på sjukhus, som är minst en dag. Ett förnuftigt tillvägagångssätt är att anpassa en Poisson- eller negativ binomialmodell som utesluter noll och ändrar skalan på de andra sannolikheterna så att de summerar till ett. Man bör vara försiktig när man tolkar dessa modeller eftersom μ inte är det förväntade resultatet utan medelvärdet av en underliggande fördelning som inkluderar nollorna. Dessa modeller implementeras i funktionen vglm() i paketet VGAM, med hjälp av familjerna pospoisson och posnegbinomial.
Ett alternativt tillvägagångssätt vid överskott (eller brist) på nollor är att använda en tvåstegsprocess, med en logit-modell för att skilja mellan nollor och positiva räkningar och sedan en nolltrunkerad Poisson- eller negativ binomialmodell för de positiva räkningarna. I vårt exempel skulle vi kunna använda en logitmodell för att skilja dem som publicerar från dem som inte gör det, och sedan en trunkerad Poisson- eller negativ binomialmodell för antalet artiklar för dem som publicerar minst en. Dessa modeller kallas ofta för "hurdle-modeller". De kan anpassas i R med hjälp av de separata logitmodellerna och de nolltrunkerade Poisson- eller negativa binomialmodellerna och genom att helt enkelt addera log-likelihoods, eller med hjälp av funktionen hurdle() i paketet pscl.
Vid jämförelse av hurdle- och nollinflaterade modeller anser jag att skillnaden mellan noll och en eller flera är tydligare med hurdle-modeller, men tolkningen av medelvärdet är tydligare med nollinflaterade modeller.
I följande JSS-artikel finns en användbar diskussion av alla dessa modeller: Regression Models for Count Data in R.