Când vine vorba de modelarea numerelor (adică a numerelor întregi mai mari sau egale cu 0), începem adesea cu regresia Poisson. Acesta este un model liniar generalizat în care se presupune că un răspuns are o distribuție Poisson condiționată de o sumă ponderată de predictori. De exemplu, am putea modela numărul de contuzii documentate la fundașii din NFL în funcție de numărul de meciuri jucate și de numărul total de ani de experiență a liniei sale ofensive. Cu toate acestea, un potențial dezavantaj al regresiei Poisson este că s-ar putea să nu descrie cu exactitate variabilitatea numărului de contuzii.

O distribuție Poisson este parametrizată prin \(\lambda\), care se întâmplă să fie atât media, cât și varianța sa. Deși este convenabil de reținut, nu este adesea realistă. O distribuție de numere va avea, de obicei, o varianță care nu este egală cu media sa. Atunci când vedem că acest lucru se întâmplă cu date pe care le presupunem (sau sperăm) că sunt distribuite Poisson, spunem că avem o subdispersie sau o supradispersie, în funcție de faptul că varianța este mai mică sau mai mare decât media. Efectuarea regresiei Poisson pe date de numărare care prezintă acest comportament are ca rezultat un model care nu se potrivește bine.

O abordare care abordează această problemă este Regresia Binomială Negativă. Distribuția binomială negativă, ca și distribuția Poisson, descrie probabilitățile de apariție a numerelor întregi mai mari sau egale cu 0. Spre deosebire de distribuția Poisson, varianța și media nu sunt echivalente. Acest lucru sugerează că ar putea servi ca o aproximare utilă pentru modelarea numerelor cu variabilitate diferită de media sa. Varianța unei distribuții binomiale negative este o funcție a mediei sale și are un parametru suplimentar, k, numit parametru de dispersie. Să spunem că numărătoarea noastră este variabila aleatoare Y dintr-o distribuție binomială negativă, atunci varianța lui Y este

$$ var(Y) = \mu + \mu^{2}/k $$

Pe măsură ce parametrul de dispersie devine din ce în ce mai mare, varianța converge la aceeași valoare cu media, iar binomul negativ se transformă într-o distribuție Poisson.

Pentru a ilustra distribuția binomială negativă, să lucrăm cu câteva date din cartea, Categorical Data Analysis, de Alan Agresti (2002). Datele sunt prezentate în tabelul 13.6 din secțiunea 13.4.3. Datele provin dintr-un sondaj efectuat pe un eșantion de 1308 persoane, în care acestea au fost întrebate câte victime de omucidere cunosc. Variabilele sunt resp, numărul de victime pe care le cunoaște respondentul, și race, rasa respondentului (negru sau alb). Ajută rasa să explice câte victime de omucidere cunoaște o persoană? Datele trebuie mai întâi să fie introduse în R:

> # Table 13.6> # Agresti, p. 561> black <- c(119,16,12,7,3,2,0)> white <- c(1070,60,14,4,0,0,1)> resp <- c(rep(0:6,times=black), rep(0:6,times=white))> race <- factor(c(rep("black", sum(black)), rep("white", sum(white))),+ levels = c("white","black"))> victim <- data.frame(resp, race)

Înainte de a ajunge la modelare, să explorăm datele. În primul rând, observăm că majoritatea respondenților sunt albi:

> table(race)racewhite black 1149 159 

Negrii au un număr mediu mai mare decât albii:

> with(victim, tapply(resp, race, mean)) white black 0.09225413 0.52201258 

Pentru fiecare rasă, varianța eșantionului este aproximativ dublă față de medie. Se pare că avem o supradispersie.

> with(victim, tapply(resp, race, var)) white black 0.1552448 1.1498288 

În cele din urmă ne uităm la distribuția numărătorii pe rase.

> table(resp, race) raceresp white black 0 1070 119 1 60 16 2 14 12 3 4 7 4 0 3 5 0 2 6 1 0

Pentru a ajusta modelul. Mai întâi încercăm regresia Poisson folosind funcția glm() și arătăm o parte din rezultatul sumar.

> # Poisson model> pGLM <- glm(resp ~ race, data=victim, family = poisson)> summary(pGLM)Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.38321 0.09713 -24.54 <2e-16 ***raceblack 1.73314 0.14657 11.82 <2e-16 ***

Rasa este foarte semnificativă. Se pare că negrii sunt mult mai predispuși să cunoască pe cineva care a fost victima unei omucideri. Dar ce înseamnă coeficientul 1,73? În acest model simplu cu un singur predictor dihotomic, este diferența de numere logaritmice așteptate. Dacă exponențiem coeficientul, obținem un raport al mediilor eșantionului:

> exp(coef(pGLM))raceblack 5.658419 > # same thing> mean(victim$resp)/mean(victim$resp) 5.658419

De fapt, dacă facem o predicție cu acest model și exponențiem rezultatele, obținem mediile eșantionului:

> exp(predict(pGLM, newdata = data.frame(race=c("white","black")))) 1 2 0.09225413 0.52201258 > # same thing> with(victim, tapply(resp, race, mean)) white black 0.09225413 0.52201258 

Acest lucru spune că numărul de victime cunoscute pentru albi este distribuit ca un Poisson cu media și varianța egale cu 0.09, în timp ce numărul de victime cunoscute pentru negri este distribuit ca un Poisson cu media și varianța egale cu 0,52. Știm deja din analiza noastră exploratorie că varianțele observate au fost mult mai mari, așa că nu ar trebui să fim prea mulțumiți de varianțele estimate de model. Dacă examinăm valorile ajustate, vom vedea și mai multe dovezi ale lipsei de potrivire:

> # fitted counts for Poisson GLM:> fmeans <- exp(predict(pGLM, newdata = data.frame(race = c("white","black"))))> fmeans 1 2 0.09225413 0.52201258 > fittedW <- dpois(0:6,lambda = fmeans) * sum(victim$race=="white") > fittedB <- dpois(0:6,lambda = fmeans) * sum(victim$race=="black") > data.frame(Response=0:6,BlackObs=black, BlackFit=round(fittedB,1), + WhiteObs=white, WhiteFit=round(fittedW,1)) Response BlackObs BlackFit WhiteObs WhiteFit1 0 119 94.3 1070 1047.72 1 16 49.2 60 96.73 2 12 12.9 14 4.54 3 7 2.2 4 0.15 4 3 0.3 0 0.06 5 2 0.0 0 0.07 6 0 0.0 1 0.0

După aceea am salvat mai întâi mediile prezise într-un obiect numit fmeans. Apoi am generat numere ajustate folosind funcția dpois împreună cu mediile estimate pentru a prezice probabilitatea de a obține de la 0 la 6. Apoi am înmulțit aceste probabilități cu numărul de respondenți pentru a obține numere ajustate. În cele din urmă, am combinat totul într-un cadru de date pentru a compara cu ușurință valorile observate și cele ajustate. Două dintre cele mai spectaculoase lucruri care trebuie observate sunt faptul că subapreciem numărul 0 și supraapreciem numărul 1.

Potem folosi o rootogramă pentru a vizualiza ajustarea unui model de regresie a numărului de persoane. Funcția rootogram() din pachetul countreg ușurează acest lucru.

> countreg::rootogram(pGLM)

Linia roșie curbă este ajustarea Poisson teoretică. “Atârnând” de fiecare punct de pe linia roșie este o bară, a cărei înălțime reprezintă diferența dintre numărătorile așteptate și cele observate. O bară care atârnă sub 0 indică o subadaptare. O bară care atârnă peste 0 indică o supraadaptare. Numărătorile au fost transformate cu o transformare a rădăcinii pătrate pentru a preveni ca numărătorile mai mici să fie întunecate și copleșite de numărători mai mari. Observăm o mare subadaptare pentru numerele 2 și mai mari și o supraadaptare masivă pentru numărul 1.

Acum să încercăm să ajustăm un model binomial negativ. Am observat că variabilitatea numărătorilor este mai mare pentru ambele rase. S-ar părea că distribuția binomială negativă ar aproxima mai bine distribuția numerelor.

Pentru a ajusta un model binomial negativ în R ne adresăm funcției glm.nb() din pachetul MASS (un pachet care vine instalat împreună cu R). Din nou, arătăm doar o parte din rezumatul ieșirii:

> # Fit a negative binomial model> library(MASS)> nbGLM <- glm.nb(resp ~ race, data=victim)> summary(nbGLM)Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.3832 0.1172 -20.335 < 2e-16 ***raceblack 1.7331 0.2385 7.268 3.66e-13 *** Theta: 0.2023 

În primul rând, observați că coeficienții sunt aceiași ca înainte. Încă o dată, putem exponențializa coeficientul de rasă pentru a obține un raport al mediilor eșantionului și putem face predicții pentru a obține mediile eșantionului original.

> # fitted counts for Negative Binomial GLM:> fmeans <- exp(predict(pGLM, newdata = data.frame(race = c("white","black"))))> fmeans # same as pGLM 1 2 0.09225413 0.52201258 

Dar observați că eroarea standard pentru coeficientul de rasă este mai mare, ceea ce indică mai multă incertitudine în estimarea noastră (0,24 față de 0,15). Acest lucru are sens, având în vedere variabilitatea observată în numărătorile noastre. De asemenea, observați estimarea lui Theta. Acesta este parametrul nostru de dispersie. Îl putem accesa direct din obiectul nostru model după cum urmează:

> nbGLM$theta 0.2023119

Și îl putem folosi pentru a obține varianțe estimate pentru numărători:

> fmeans + fmeans^2 * (1/nbGLM$theta) 1 2 0.134322 1.868928 

Acestea sunt mult mai apropiate de varianțele observate decât cele date de modelul Poisson.

Încă o dată vizualizăm ajustarea folosind o rootogramă:

> countreg::rootogram(nbGLM)

Aceasta arată mult mai bine decât rootograma modelului Poisson. Există o ușoară subadaptare/supraadaptare pentru numărătorile de la 1 la 3, dar în rest arată destul de bine.

Pentru a obține mai multe informații despre modelul nostru binomial negativ, să folosim parametrii acestuia pentru a simula date și să comparăm datele simulate cu datele observate. Mai jos vom încărca pachetul magrittr pentru a avea acces la operatorul %>% care ne permite să “înlănțuim” funcții.

În primul rând, tabulăm numărătorile și creăm un barplot pentru participanții albi și, respectiv, negri. Apoi folosim parametrii modelului pentru a simula datele dintr-o distribuție binomială negativă. Cei doi parametri sunt mu și size (adică parametrul de dispersie). Observați că folosim funcția coef() pentru a extrage coeficienții corespunzători pentru fiecare rasă. Pentru alb este doar intercepția, pentru negru este vorba de intercepție și pantă (deci le adunăm). Apoi facem o exponențializare pentru a converti de la scara logaritmică la scara originală. Simulăm același număr de observații pe care îl avem în datele noastre originale.

library(magrittr) # for %>% operatorop <- par(mfrow=c(1,2))set.seed(1)victim$resp %>% `)) %>% table() %>% barplot(main = "Simulated White")victim$resp %>% `[`(victim$race=="black") %>% table() %>% barplot(main = "Observed Black")rnbinom(n = 159, size = nbGLM$theta, mu = exp(sum(coef(nbGLM)))) %>% table() %>% barplot(main = "Simulated Black")par(op)

Datele simulate sunt foarte asemănătoare cu datele observate, ceea ce ne dă din nou încredere în alegerea regresiei binomiale negative pentru a modela aceste date. Aceste diagrame demonstrează, de asemenea, natura condiționată a modelului nostru. Distribuția binomială negativă a numerelor depinde, sau este condiționată de rasă. Fiecare rasă are o medie diferită, dar un parametru de dispersie comun.

• Agresti, Alan (2002), Categorical Data Analysis, Wiley.
• Kleiber, Christian & Zeileis, Achim (2016): Visualizing Count Data Regressions Using Rootograms, The American Statistician, DOI: 10.1080/00031305.2016.1173590

Pentru întrebări sau clarificări cu privire la acest articol, contactați UVA Library StatLab: [email protected]

Vezi întreaga colecție de articole din UVA Library StatLab.

Clay Ford
Consilier de cercetare statistică
University of Virginia Library
.