GR’s Website

4.A Modele pentru date de numărare supra-dispersate

Utilizăm datele din Long (1990) privind numărul de publicații produse de biochimiștii cu doctorat pentru a ilustra aplicarea modelelor Poisson, Poisson supra-dispersat, binomial negativ și Poisson cu inflație zero.

Variabilele din setul de date sunt

articole în ultimii trei ani de doctorat.

fem codificat unu pentru femei mar codificat unu dacă este căsătorit. kid5 numărul de copii sub șase ani

prestigiul doctoratului .D. program

articole de către mentor în ultimii trei ani

Variabilă	Descriere
art
phd
ment

Aceste date au fost, de asemenea, analizate de Long și Freese (2001) și sunt disponibile pe site-ul 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 

Numărul mediu de articole este de 1,69, iar varianța este de 3,71, puțin mai mult decât dublul mediei. Datele sunt supra-dispersate, dar, bineînțeles, nu am luat în considerare încă nicio covariantă.

Un model Poisson

Să potrivim modelul utilizat de Long și Freese(2001), un model aditiv simplu care utilizează toți cei cinci predictori.

> 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
și devianța și chi pătrat al lui Pearson sunt amândouă în jurul valorii de 1600.
 4.A.1. Variația extra-Poisson
Acum presupunem că varianța este mai degrabă proporțională decât egală cu media și estimăm parametrul de scară φ împărțind chi pătrat al lui Pearson la f.d. al acestuia:
> phi round(c(phi,sqrt(phi)),4) 1.8290 1.3524
Vezi că varianța este cu aproximativ 83% mai mare decât media. Acest lucru înseamnă că ar trebui să ajustăm erorile standard înmulțind cu 1,35, rădăcina pătrată a lui 1,83.
R poate face acest calcul pentru noi dacă folosim quasipoissonfamilia:
> 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
O abordare alternativă este de a ajusta un model Poisson și de a folosi estimatorul robust sau sandwich al erorilor standard. Acest lucru dă de obicei rezultate foarte asemănătoare cu modelul Poisson supra-dispersat.
 4.A.2 Regresia binomială negativă
Acum vom ajusta un model binomial negativ cu aceiași predictori.Pentru a face acest lucru avem nevoie de funcția glm.nb() din pachetulMASS.
> 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 thetaeste precizia efectului multiplicativerandom, și corespunde la 1/σ2 în cadrul notelor. Estimarea corespunde unei varianțe estimate de0,44 și este foarte semnificativă.
Pentru a testa semnificația acestui parametru, vă puteți gândi să calculați de două ori diferența de probabilități logaritmice între acest model și modelul Poisson, 180,2, și să o tratați ca un chi pătrat cu o f.d. Cu toate acestea, asimptotica obișnuită nu se aplică, deoarece ipoteza nulă se află pe o limită a spațiului parametrilor.
Există unele lucrări care arată că o aproximare mai bună este tratarea statisticii ca un amestec 50:50 de zero și un chi pătrat cu o f.d. Alternativ, tratarea statisticii ca un chi pătrat cu o f.d. oferă un test conservator. Oricum ar fi, avem dovezi copleșitoare de supradispersie.
Pentru testarea ipotezelor privind coeficienții de regresie putem folosi fie testele Wald, fie testele raportului de verosimilitate, care sunt posibile deoarece am făcut ipoteze distribuționale complete.
Există, de asemenea, o familie negative.binomial pentru glm și aceasta poate fi folosită cu condiția ca parametriitheta să fie dați. (Aceasta se bazează pe rezultatul conform căruia binomul negativ se află în familia glm pentru varianță fixă). Iată o verificare rapidă:
> 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)

De asemenea, putem calcula cuantele folosindqgamma(p, shape, rate = 1, scale = 1/rate, lower.tail = TRUE). În special, atunci când efectul aleator are varianța v, cuartilele sunt 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
Biochimiștii de la Q1 din distribuția eterogenității neobservate publică cu 49% mai puține lucrări decât se așteaptă din caracteristicile lor observate, în timp ce cei de la mediană publică cu 14% mai puține și cei de la Q3 publică cu 33% mai mult decât se așteaptă.
Compararea estimărilor și a erorilor standard
Să comparăm estimările parametrilor și erorile standard în cadrul modelelor Poisson, Poisson supradispersat și binomial negativ:
> 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
Stimările binomialului negativ nu sunt foarte diferite de cele bazate pe modelul Poisson, iar ambele seturi ar duce la aceleașiconcluzii. 
Cu privire la erorile standard, observăm că ambele abordări ale supra-dispersiei conduc la erori standard estimate foarte asemănătoare,și că regresia Poisson obișnuită subestimează erorile standard.
Bunătatea de potrivire
Potem evalua bunătatea de potrivire a modelului binomial negativfolosind devianța
> deviance(mnbg) 1004.281
Modelul binomial negativ se potrivește mai bine decât cel Poisson, dar are totuși o devianță peste valoarea critică de cinci procente de 980,25.
Funcția de varianță
Modelul Poisson supra-dispersat și modelul binomial negativ au funcții de varianță diferite. O modalitate de a verifica care dintre ele poate fi mai potrivită este de a crea grupuri bazate pe predictorul liniar, de a calcula media și varianța pentru fiecare grup și, în final, de a reprezenta grafic relația medie-varianță.
Iată grupurile bazate pe predictorul liniar binomial negativ, create folosind cut cu întreruperi la percentilele 5(5)95pentru a produce douăzeci de grupuri de dimensiuni aproximativ egale. Rețineți că predict() calculează un predictor liniar,dacă nu se specifică altfel. Pentru a prezice în scara răspunsului, adăugați opțiunea 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 

Graficul prezintă media în raport cu varianța și suprapune curbele corespunzătoare modelului Poisson supradispersat, în care varianța este φμ, și modelului binomial negativ, în care varianța este μ(1+μσ2).
Funcția de varianță Poisson face o treabă destul de bună pentru majoritatea datelor, dar nu reușește să surprindă varianțele ridicate ale celor mai productivi cercetători. Funcția de varianță binomială negativă nu este prea diferită, dar, fiind o funcție pătratică, poate crește mai repede și face o treabă mai bună la limita superioară. Concluzionăm că modelul binomial negativ oferă o descriere mai bună a datelor decât modelul Poisson supradispersat. 
 4.A.3 Modele Poisson cu inflație de zerouri
 O apariție frecventă în cazul datelor de numărare este un exces de zerouri în comparație cu ceea ce este de așteptat în cadrul unui model Poisson. Aceasta este, de fapt, o problemă cu datele noastre: 
> zobs zpoi c(obs=mean(zobs), poi=mean(zpoi)) obs poi 0.3005464 0.2092071 
Vezi că 30,0% dintre oamenii de știință din eșantion nu au publicat niciun articol în ultimii trei ani de la doctorat, dar modelul Poisson prezice că doar 20,9% nu ar avea nicio publicație. În mod clar, modelul subestimează probabilitatea de numărare zero.
O modalitate de a modela acest tip de situație este să presupunem că datele provin dintr-un amestec de două populații, una în care numărarea este întotdeauna zero, iar cealaltă în care numărarea are o distribuție Poisson cu media μ. În acest model, numărătorile zero pot proveni din oricare dintre populații, în timp ce numărătorile pozitive provin doar din cea de-a doua. 
În contextul publicațiilor doctoranzilor în biochimie, ne putem imagina că unii aveau în vedere locuri de muncă în care publicațiile nu ar fi fost importante, în timp ce alții vizau locuri de muncă academice în care se aștepta un record de publicații. Membrii primului grup ar fi publicat zero articole, în timp ce membrii celui de-al doilea grup ar fi publicat 0,1,2,..., un număr despre care se poate presupune că are o distribuție Poisson.
Distribuția rezultatului poate fi apoi modelată în funcție de doi parametri, π probabilitatea de "întotdeauna zero" și μ, numărul mediu de publicații pentru cei care nu fac parte din grupul "întotdeauna zero". Un mod natural de a introduce covariate este de a modela logit al probabilității π de a avea întotdeauna zero și logaritmul mediei μ pentru cei care nu fac parte din clasa "întotdeauna zero".
Acest tip de model poate fi ajustat în R utilizând funcția zeroinfl() din pachetul pscl. Formula modelului poate fi specificată ca de obicei dacă aceleași variabile urmează să fie incluse în ambele ecuații. În caz contrar, se pot furniza două seturi de predictoriseparate de o bară verticală, tastați ?zeroinfl pentru a afla mai multe.
Iată un model umflat cu zero cu toate covariatele în ambele ecuații:
> 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
Urmărind ecuația umflată, vedem că singurul predictor semnificativ al apartenenței la clasa "întotdeauna zero" este numărul de articole publicate de mentor, fiecare articol al mentorului fiind asociat cu șanse cu 12,6% mai mici de a nu publica niciodată.
Urmărind ecuația pentru numărul mediu sau numărul mediu de articole în rândul celor care nu se află în clasa "întotdeauna zero", găsim dezavantaje semnificative pentru femeile și oamenii de știință cu copii sub cinci ani, precum și un efect pozitiv mare al numărului de publicații ale mentorului, fiecare articol fiind asociat cu o creștere de 1,8% a numărului așteptat de publicații.
Pentru a verifica dacă modelul rezolvă problema excesului de zerouri, prezicem π și μ și calculăm probabilitatea combinată de a nu publica. Există opțiuni în predict()funcție numite "zero" și "count" pentru a le obține. Există, de asemenea, o opțiune "prob" pentru a calcula densitatea prezisă, dar aceasta este exagerată, deoarece dorim doar probabilitatea de zero.
> pr mu zip mean(zip) 0.2985679
Acum, modelul rezolvă problema excesului de zerouri, prezicând că 29,9% dintre biochimiști nu vor publica niciun articol, mult mai aproape de valoarea observată de 30,0%.
Există, de asemenea, un model binomial negativ umflat cu zero, care utilizează un binom negativ pentru numărul din clasa "nu întotdeauna zero". Acest model poate fi ajustat folosind zeroinfl() cu parametrul dist="negbin". Legăturile alternative pentru ecuația de inflație includ probit,care poate fi specificat folosind link="probit". 
Compararea modelelor cu AIC
Cum se întâmplă, pentru aceste date, binomul negativ rezolvă și el problema! Iată probablitatea articolelor zero în binomul negativ. Noi procedăm de la primele principii, dar s-ar putea folosi și densitatea binomului negativ încorporat
> munb theta znb # also dnbinom(0, mu=munb, size=theta)> mean(znb) 0.3035957
Modelul prezice că 30,4% dintre biochimiști nu vor publica niciun articol în ultimii trei ani de doctorat, foarte aproape de valoarea observată de 30,0%.
Pentru a alege între modelul binomului negativ și cel cu umflătură zero trebuie să recurgem la alte criterii. Un mod foarte simplu de a compara modele cu numere diferite de parametri este calcularea criteriului informațional Akaike (AIC), pe care îl definim ca
AIC = -2logL + 2p
unde p este numărul de parametri din model. Primul termen este, în esență, devianța, iar al doilea este o penalizare pentru numărul de parametri și poate fi calculat pentru majoritatea modelelor cu ajutorul funcției AIC. De asemenea, îl obțin "de mână" pentru a putea verifica calculul. Pentru datele noastre
> 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
 Pentru acest set de date, modelul binomial negativ este un câștigător clar în ceea ce privește parsimonia și buna potrivire. 
Alte criterii de diagnostic pe care le-am putea analiza sunt distribuția marginală a numerelor prezise și observate și funcțiile de varianță. 
Modelurile Zero-Truncated și Hurdle
Alte modele pe care nu le-am abordat sunt Poisson trunchiat cu zero și binomial negativ, concepute pentru date care nu includ zerouri. Un exemplu obișnuit este durata de ședere într-un spital, care este de cel puțin o zi. O abordare rezonabilă este să se potrivească un model Poisson sau un model binomial negativ care exclude zero și redimensionează celelalte probabilități pentru a le însuma la unu. Trebuie să fim atenți la interpretarea acestor modele, deoarece μ nu este rezultatul așteptat, ci media unei distribuții subiacente care include zerourile. Aceste modele sunt implementate în funcția vglm() din pachetul VGAM, utilizând familiile pospoisson și posnegbinomial.
O abordare alternativă a excesului (sau a lipsei) de zerouri este de a utiliza un proces în două etape, cu un model logit pentru a face distincția între numerele zero și cele pozitive și apoi un model Poisson trunchiat cu zero sau un model binomial negativ pentru numerele pozitive. În exemplul nostru, am putea utiliza un model logit pentru a-i diferenția pe cei care publică de cei care nu publică și apoi un model Poisson trunchiat sau un model binomial negativ pentru numărul de articole ale celor care publică cel puțin unul. Aceste modele sunt adesea numite modele hurdle. Ele pot fi ajustate în R folosind modelele logit separate și modelele Poisson sau binomiale negative trunchiate cu zero și adăugând pur și simplu probabilitățile logaritmice, sau folosind funcția hurdle() din pachetul pscl.
Comparând modelele cu hurdle și modelele cu inflexiune zero, consider că distincția dintre zero și unul sau mai multe este mai clară cu modelele cu hurdle, dar interpretarea mediei este mai clară cu modelele cu inflexiune zero.
 Următoarea lucrare JSS conține o discuție utilă despre toate aceste modele: Regression Models for Count Data in R.
.