GR’s Website

4.A Modellen voor overgedispergeerde telgegevens

We gebruiken gegevens van Long (1990) over het aantal publicaties van gepromoveerde biochemici om de toepassing van Poisson-, overgedispergeerde Poisson-, negatief binomiaal- en zero-inflated Poisson-modellen te illustreren.

De variabelen in de dataset zijn

mar

Variabele	Beschrijving
art	artikels in laatste drie jaar van Ph.D.
fem	één gecodeerd voor vrouwen
één gecodeerd indien gehuwd
kid5	aantal kinderen jonger dan zes jaar
phd	prestige van Ph.D. programma
ment	artikelen door mentor in laatste drie jaar

Deze gegevens zijn ook geanalyseerd door Long en Freese (2001), en zijn beschikbaar op de Stata website:

> 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 

Het gemiddelde aantal artikelen is 1,69 en de variantie is 3,71, iets meer dan tweemaal het gemiddelde. De gegevens zijn overgedispergeerd, maar we hebben natuurlijk nog geen rekening gehouden met covariaten.

Een Poisson Model

Laten we het model passen dat Long en Freese(2001) gebruiken, een eenvoudig additief model met alle vijf de voorspellers.

> 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
en de deviantie en Pearsons chi-kwadraat liggen beide in de 1600s.
 4.A.1. Extra-poissonvariatie
We nemen nu aan dat de variantie evenredig is met het gemiddelde in plaats van gelijk, en schatten de schaalparameter φ door Pearsons chi-kwadraat te delen door zijn d.f.:
> phi round(c(phi,sqrt(phi)),4) 1.8290 1.3524
We zien dat de variantie ongeveer 83% groter is dan het gemiddelde. Dit betekent dat we de standaardfouten moeten aanpassen door ze te vermenigvuldigen met 1,35, de vierkantswortel van 1,83.
R kan deze berekening voor ons uitvoeren als we de quasipoissonfamilie gebruiken:
> 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
Een alternatieve aanpak is het fitten van een Poisson-model en het gebruik van de robuuste of sandwich-schatter van de standaardfouten. Dit levert gewoonlijk resultaten op die sterk lijken op die van het overgedispergeerde Poisson-model.
 4.A.2 Negatief binomiale regressie
We passen nu een negatief binomiaal model toe met dezelfde voorspellers. Hiervoor hebben we de glm.nb() functie uit hetMASS pakket nodig.
> 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 is de precisie van het multiplicativerandom effect, en komt overeen met 1/σ2 in thenotes. De schatting komt overeen met een geschatte variantie van 0,44 en is zeer significant.
Om de significantie van deze parameter te testen kunt u denken aan een berekening van tweemaal het verschil in log-likelihoods tussen dit model en het Poisson-model, 180,2, en dit te behandelen als een chi-kwadraat met één d.f. De gebruikelijke asymptotiek is echter niet van toepassing, omdat de nulhypothese op een grens van de parameterruimte ligt.
Er is enig werk waaruit blijkt dat een betere benadering erin bestaat de statistiek te behandelen als een 50:50 mengsel van nul en een chi-kwadraat met één d.f. Een andere mogelijkheid is de statistiek te behandelen als een chi-kwadraat met één d.f., hetgeen een conservatieve test oplevert. Hoe dan ook, we hebben een overweldigend bewijs van overdispersie.
Voor het testen van hypothesen over de regressiecoëfficiënten kunnen we Wald-tests of likelihoodratio-tests gebruiken, die mogelijk zijn omdat we volledige verdelingsaannamen hebben gemaakt.
Er is ook een negative.binomial-familie voor glm en deze kan worden gebruikt mits de parameterstheta zijn gegeven. (Dit is gebaseerd op het resultaat dat de negatieve binomiaal in de glm-familie zit voor vaste variantie). Hier volgt een snelle controle:
> 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)

We kunnen de kwantielen ook berekenen metqgamma(p, shape, rate = 1, scale = 1/rate, lower.tail = TRUE). In het bijzonder, wanneer het willekeurige effect variantie v heeft, zijn de kwartielen 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 op Q1 van de verdeling van niet-waargenomen heterogeniteit publiceren 49% minder papers dan verwacht op basis van hun waargenomen kenmerken, terwijl degenen op de mediaan 14% minder publiceren en degenen op Q3 33% meer publiceren dan verwacht.
Vergelijking van schattingen en standaardfouten
Laten we de parameterschattingen en standaardfouten onder het Poisson-model, het overgedispergeerde Poisson-model en het negatieve binomiale model vergelijken:
> 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 negatieve binomiale schattingen verschillen niet veel van die op basis van het Poisson-model, en beide reeksen zouden tot dezelfde conclusies leiden. 
Kijken we naar de standaardfouten, dan zien we dat beide benaderingen van overspreiding tot zeer vergelijkbare geschatte standaardfouten leiden, en dat de gewone Poisson-regressie de standaardfouten onderschat.
Goodness of Fit
We kunnen de goodness of fit van het negatieve binomiale model beoordelen met behulp van de deviantie
> deviance(mnbg) 1004.281
Het negatieve binomiale model past beter dan het Poisson-model, maar heeft nog steeds een deviantie boven de kritische waarde van vijf procent van980,25.
De variantiefunctie
De overgedispergeerde Poisson- en negatieve binomiale modellen hebben verschillende variantiefuncties. Een manier om na te gaan welk model het meest geschikt is, is groepen te maken op basis van de lineaire voorspeller, het gemiddelde en de variantie voor elke groep te berekenen en ten slotte de verhouding tussen gemiddelde en variantie te plotten.
Hier volgen groepen op basis van de negatieve binomiale lineaire voorspeller, gemaakt met cut met breuken bij de 5(5)95-percentielen om twintig groepen van ongeveer gelijke grootte te maken. Merk op dat predict() een lineaire voorspeller berekent, tenzij anders gespecificeerd. Voeg de optie type="response" toe om in de schaal van de respons te voorspellen. 
> 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 

De grafiek zet het gemiddelde uit tegen de variantie en toont de krommen die overeenkomen met het overgedispergeerde Poisson-model, waarbij de variantie φμ is, en het negatieve binomiale model, waarbij de variantie μ(1+μσ2) is.
De Poisson-variantiefunctie doet haar werk vrij goed voor het grootste deel van de gegevens, maar slaagt er niet in de hoge varianties van de meest productieve geleerden weer te geven. De negatieve binomiale variantiefunctie wijkt niet al te veel af, maar is kwadratisch, kan sneller stijgen en doet het beter in de hoogste regionen. We concluderen dat het negatief binomiaal model een betere beschrijving van de gegevens geeft dan het overgedispergeerde Poisson model.
 4.A.3 Zero-Inflated Poisson Models
 Een vaak voorkomend verschijnsel bij telgegevens is een overmaat aan nullen in vergelijking met wat wordt verwacht onder een Poisson model. Dit is in feite een probleem met onze gegevens: 
> zobs zpoi c(obs=mean(zobs), poi=mean(zpoi)) obs poi 0.3005464 0.2092071 
We zien dat 30,0% van de wetenschappers in de steekproef geen artikelen heeft gepubliceerd in de laatste drie jaar van hun promotieonderzoek, maar het Poisson-model voorspelt dat slechts 20,9% geen publicaties zou hebben. Het is duidelijk dat het model de kans op nultellingen onderschat.
Een manier om dit soort situaties te modelleren is aan te nemen dat de gegevens afkomstig zijn van een mengsel van twee populaties, een waarbij de telling altijd nul is, en een andere waarbij de telling een Poisson-verdeling heeft met gemiddelde μ. In dit model kunnen nul tellingen uit beide populaties komen, terwijl positieve tellingen alleen uit de tweede populatie komen.
In de context van publicaties door gepromoveerde biochemici kunnen we ons voorstellen dat sommigen een baan in gedachten hadden waar publicaties niet belangrijk zouden zijn, terwijl anderen een academische baan ambieerden waar een record aan publicaties werd verwacht. Leden van de eerste groep zouden nul artikelen publiceren, terwijl leden van de tweede groep 0,1,2,... zouden publiceren, een aantal waarvan mag worden aangenomen dat het Poisson verdeeld is.
De verdeling van de uitkomst kan dan worden gemodelleerd in termen van twee parameters, π de kans op "altijd nul", en μ, het gemiddelde aantal publicaties voor degenen die niet tot de "altijd nul"-groep behoren. Een natuurlijke manier om covariaten te introduceren is het modelleren van de logit van de kans π op altijd nul en de log van het gemiddelde μ voor degenen die niet tot de altijd nul-groep behoren.
Dit type model kan in R worden gefit met behulp van de functie zeroinfl() in het pakket pscl. De modelformule kan op de gebruikelijke wijze worden gespecificeerd als dezelfde variabelen in beide vergelijkingen moeten worden opgenomen. Anders kan men twee reeksen voorspellers opgeven, gescheiden door een verticale balk, type ?zeroinfl om meer te weten te komen.
Hier volgt een zero-inflated model met alle covariaten in beide vergelijkingen:
> 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
Kijkend naar de inflate vergelijking zien we dat de enige significante voorspeller van het behoren tot de klasse 'altijd nul' het aantal door de mentor gepubliceerde artikelen is, waarbij elk artikel van de mentor geassocieerd is met 12,6% minder kans om nooit te publiceren.
Kijken we naar de vergelijking voor het gemiddelde aantal artikelen onder degenen die niet in de "altijd nul"-klasse vallen, dan vinden we significante nadelen voor vrouwen en wetenschappers met kinderen jonger dan vijf jaar, en een groot positief effect van het aantal publicaties door de mentor, waarbij elk artikel geassocieerd is met een toename van 1,8% in het verwachte aantal publicaties.
Om te verifiëren dat het model het probleem van de overtollige nullen oplost, voorspellen we π en μ, en berekenen we de gecombineerde kans op geen publicaties. Er zijn opties in de predict()functie genaamd "zero" en "count" om deze te verkrijgen. Er is ook een optie "prob" om de voorspelde dichtheid te berekenen, maar dit is overkill omdat we alleen de kans op nul willen.
> pr mu zip mean(zip) 0.2985679
Dus het model lost het probleem van de overmaat aan nullen op door te voorspellen dat 29,9% van de biochemici geen artikelen zal publiceren, wat veel dichter bij de waargenomen waarde van 30,0% ligt.
Er is ook een nul-geïnflateerd negatief binomiaal model, dat een negatief binomiaal gebruikt voor de telling in de klasse "niet altijd nul". Dit model kan worden aangepast met zeroinfl() met de parameter dist="negbin". Alternatieve koppelingen voor de opblaasvergelijking zijn onder meer de probit, die kan worden gespecificeerd met link="probit". 
Modelvergelijking met AIC
 Toevallig lost voor deze gegevens de negatieve binomiaal het probleem ook op! Hier is de waarschijnlijkheid van nul artikelen in het negatieve binomium. We gaan uit van eerste principes, maar men zou ook de ingebouwde negatieve binomiale dichtheid kunnen gebruiken
> munb theta znb # also dnbinom(0, mu=munb, size=theta)> mean(znb) 0.3035957
Het model voorspelt dat 30,4% van de biochemici in de laatste drie jaar van hun doctoraat geen artikels zouden publiceren, wat zeer dicht ligt bij de waargenomen waarde van 30,0%.
 Om te kiezen tussen het negatieve binomiale model en het nulopgeblazen model moeten we onze toevlucht nemen tot andere criteria. Een zeer eenvoudige manier om modellen met verschillende aantallen parameters te vergelijken is de berekening van het Akaike's informatiecriterium (AIC), dat we definiëren als
AIC = -2logL + 2p
waarbij p het aantal parameters in het model is. De eerste term is in wezen de deviantie en de tweede een straf voor het aantal parameters, en hij kan voor de meeste modellen worden berekend met de functie AIC. Ik bereken het ook "met de hand", zodat we de berekening kunnen verifiëren. Voor onze gegevens
> 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
 is voor deze dataset het negatieve binomiale model een duidelijke winnaar in termen van parsimony en goodness of fit. 
Andere diagnostische criteria waarnaar we kunnen kijken, zijn de marginale verdeling van voorspelde en waargenomen tellingen en de variantiefuncties. 
Nul-getruncete en Hurdle-modellen
Andere modellen die we nog niet hebben behandeld, zijn de nul-getruncete Poisson en negatieve binomial, ontworpen voor gegevens die geen nullen bevatten. Een veelvoorkomend voorbeeld is de verblijfsduur in een ziekenhuis, die ten minste één dag bedraagt. Een verstandige aanpak is het fitten van een Poisson of negatief binomiaal model waarin geen nullen voorkomen en de andere kansen worden geschaald tot één. Deze modellen moeten voorzichtig worden geïnterpreteerd omdat μ niet de verwachte uitkomst is, maar het gemiddelde van een onderliggende verdeling waarin de nullen zijn opgenomen. Deze modellen zijn geïmplementeerd in de functie vglm() in het pakket VGAM, met gebruikmaking van de families pospoisson en posnegbinomial.
Een alternatieve benadering van een overmaat (of een tekort) aan nullen is het gebruik van een tweefasenproces, met een logit-model om onderscheid te maken tussen nul- en positieve tellingen en vervolgens een zero-truncated Poisson of negatief binomiaal model voor de positieve tellingen. In ons voorbeeld zouden we een logit-model kunnen gebruiken om onderscheid te maken tussen degenen die publiceren en degenen die dat niet doen, en vervolgens een afgekapt Poisson- of negatief binomiaal model voor het aantal artikelen van degenen die ten minste één artikel publiceren. Deze modellen worden vaak hurdle-modellen genoemd. Zij kunnen in R worden ingepast met gebruikmaking van de afzonderlijke logit- en zero-getruncated Poisson of negatief binomiaal modellen en door eenvoudig de log-likelihoods op te tellen, of met gebruikmaking van de hurdle()-functie in het pscl-pakket.
In vergelijking met hurdle- en zero-inflated-modellen vind ik het onderscheid tussen nul en één of meer duidelijker bij hurdle-modellen, maar de interpretatie van het gemiddelde is duidelijker bij zero-inflated-modellen.
Het volgende JSS-artikel bevat een nuttige bespreking van al deze modellen: Regressiemodellen voor telgegevens in R.