Quando si tratta di modellare i conteggi (cioè i numeri interi maggiori o uguali a 0), spesso si inizia con la regressione di Poisson. Questo è un modello lineare generalizzato in cui si assume che una risposta abbia una distribuzione di Poisson condizionata da una somma ponderata di predittori. Per esempio, potremmo modellare il numero di commozioni cerebrali documentate ai quarterback della NFL in funzione degli snap giocati e degli anni totali di esperienza della sua linea offensiva. Tuttavia un potenziale svantaggio della regressione di Poisson è che potrebbe non descrivere accuratamente la variabilità dei conteggi.

Una distribuzione di Poisson è parametrizzata da \(\lambda\), che è sia la sua media che la varianza. Anche se comodo da ricordare, spesso non è realistico. Una distribuzione di conteggi avrà di solito una varianza che non è uguale alla sua media. Quando vediamo questo accadere con dati che assumiamo (o speriamo) siano distribuiti a Poisson, diciamo che abbiamo una sotto- o sovradispersione, a seconda che la varianza sia più piccola o più grande della media. Eseguire la regressione di Poisson su dati di conteggio che mostrano questo comportamento risulta in un modello che non si adatta bene.

Un approccio che affronta questo problema è la regressione binomiale negativa. La distribuzione binomiale negativa, come la distribuzione di Poisson, descrive le probabilità del verificarsi di numeri interi maggiori o uguali a 0. A differenza della distribuzione di Poisson, la varianza e la media non sono equivalenti. Questo suggerisce che potrebbe servire come un’approssimazione utile per modellare i conteggi con variabilità diversa dalla sua media. La varianza di una distribuzione binomiale negativa è una funzione della sua media e ha un parametro aggiuntivo, k, chiamato parametro di dispersione. Diciamo che il nostro conteggio è una variabile casuale Y da una distribuzione binomiale negativa, allora la varianza di Y è

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

Come il parametro di dispersione diventa sempre più grande, la varianza converge allo stesso valore della media, e la binomiale negativa si trasforma in una distribuzione Poisson.

Per illustrare la distribuzione binomiale negativa, lavoriamo con alcuni dati dal libro, Categorical Data Analysis, di Alan Agresti (2002). I dati sono presentati nella tabella 13.6 nella sezione 13.4.3. I dati provengono da un sondaggio di 1308 persone in cui è stato chiesto loro quante vittime di omicidio conoscono. Le variabili sono resp, il numero di vittime che l’intervistato conosce, e race, la razza dell’intervistato (bianco o nero). La razza aiuta a spiegare quante vittime di omicidio una persona conosce? I dati devono prima essere inseriti in 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)

Prima di arrivare alla modellazione, esploriamo i dati. Per prima cosa notiamo che la maggior parte degli intervistati sono bianchi:

> table(race)racewhite black 1149 159 

I neri hanno una media più alta dei bianchi:

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

Per ogni razza la varianza del campione è circa il doppio della media. Sembra che abbiamo una sovradispersione.

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

Infine guardiamo la distribuzione dei conteggi per razza.

> 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

Per quanto riguarda l’adattamento del modello. Prima proviamo la regressione di Poisson usando la funzione glm() e mostriamo una parte dell’output riassuntivo.

> # 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 ***

La razza è molto significativa. Sembra che i neri abbiano molte più probabilità di conoscere qualcuno che è stato vittima di un omicidio. Ma cosa significa il coefficiente 1,73? In questo semplice modello con un predittore dicotomico, è la differenza nei conteggi log attesi. Se esponenziamo il coefficiente otteniamo un rapporto di medie campionarie:

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

Infatti se facciamo una previsione con questo modello ed esponenziamo i risultati, otteniamo le medie campionarie:

> 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 

Questo dice che il conteggio delle vittime note per i bianchi è distribuito come un Poisson con media e varianza uguali a 0.09, mentre il conteggio delle vittime conosciute per i neri è distribuito come un Poisson con media e varianza pari a 0,52. Sappiamo già dalla nostra analisi esplorativa che le varianze osservate erano molto più grandi, quindi non dovremmo essere troppo soddisfatti delle varianze stimate dal modello. Se esaminiamo i conteggi montati, vedremo ancora più prove della mancanza di adattamento:

> # 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

Sopra abbiamo prima salvato le medie previste in un oggetto chiamato fmeans. Abbiamo poi generato dei conteggi montati usando la funzione dpois insieme ai mezzi stimati per prevedere la probabilità di ottenere da 0 a 6. Abbiamo poi moltiplicato queste probabilità per il numero di intervistati per ottenere i conteggi adattati. Infine abbiamo combinato il tutto in un frame di dati per confrontare facilmente i valori osservati e quelli adattati. Due delle cose più drammatiche da notare è che stiamo sottoadattando i conteggi 0 e sovraadattando i conteggi 1.

Possiamo usare un rootogramma per visualizzare l’adattamento di un modello di regressione di conteggio. La funzione rootogram() nel pacchetto countreg lo rende facile.

> countreg::rootogram(pGLM)

La linea rossa curva è l’adattamento teorico di Poisson. “Appesa” ad ogni punto della linea rossa c’è una barra, la cui altezza rappresenta la differenza tra i conteggi previsti e quelli osservati. Una barra appesa sotto lo 0 indica un adattamento insufficiente. Una barra appesa sopra lo 0 indica un overfitting. I conteggi sono stati trasformati con una trasformazione della radice quadrata per evitare che i conteggi più piccoli vengano oscurati e sopraffatti da quelli più grandi. Vediamo una grande quantità di underfitting per i conteggi 2 e superiori e un massiccio overfitting per il conteggio 1.

Ora proviamo ad applicare un modello binomiale negativo. Abbiamo notato che la variabilità dei conteggi è maggiore per entrambe le razze. Sembrerebbe che la distribuzione binomiale negativa approssimi meglio la distribuzione dei conteggi.

Per adattare un modello binomiale negativo in R ci rivolgiamo alla funzione glm.nb() nel pacchetto MASS (un pacchetto che viene installato con R). Ancora una volta mostriamo solo una parte dell’output riassuntivo:

> # 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 

Prima di tutto notiamo che i coefficienti sono gli stessi di prima. Ancora una volta possiamo esponenziare il coefficiente di gara per ottenere un rapporto di medie campionarie e fare previsioni per ottenere le medie campionarie originali.

> # 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 

Ma si noti che l’errore standard per il coefficiente di razza è più grande, indicando più incertezza nella nostra stima (0,24 contro 0,15). Questo ha senso data la variabilità osservata nei nostri conteggi. Notate anche la stima di Theta. Questo è il nostro parametro di dispersione. Possiamo accedervi direttamente dal nostro oggetto modello come segue:

> nbGLM$theta 0.2023119

E possiamo usarlo per ottenere varianze stimate per i conteggi:

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

Queste sono molto più vicine alle varianze osservate di quelle date dal modello Poisson.

Ancora una volta visualizziamo il fit usando un rootogramma:

> countreg::rootogram(nbGLM)

Questo sembra molto meglio del rootogramma del modello di Poisson. C’è un leggero underfitting/overfitting per i conteggi da 1 a 3, ma altrimenti sembra abbastanza buono.

Per ottenere ulteriori informazioni sul nostro modello binomiale negativo, usiamo i suoi parametri per simulare i dati e confrontare i dati simulati con quelli osservati. Di seguito carichiamo il pacchetto magrittr per accedere all’operatore %>% che ci permette di “concatenare” le funzioni.

Prima tabuliamo i conteggi e creiamo un grafico a barre per i partecipanti bianchi e neri, rispettivamente. Poi usiamo i parametri del modello per simulare i dati di una distribuzione binomiale negativa. I due parametri sono mu e size (cioè il parametro di dispersione). Notate che usiamo la funzione coef() per estrarre i coefficienti appropriati per ogni razza. Per i bianchi è solo l’intercetta, per i neri è l’intercetta e la pendenza (quindi li sommiamo). Poi esponiamo per convertire dalla scala log alla scala originale. Simuliamo lo stesso numero di osservazioni che abbiamo nei nostri dati originali.

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)

I dati simulati sono molto simili ai dati osservati, dandoci nuovamente fiducia nella scelta della regressione binomiale negativa per modellare questi dati. Questi grafici dimostrano anche la natura condizionale del nostro modello. La distribuzione binomiale negativa dei conteggi dipende, o è condizionata, dalla razza. Ogni razza ha una media diversa ma un parametro di dispersione comune.

• 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

Per domande o chiarimenti su questo articolo, contattare la UVA Library StatLab: [email protected]

Vedi l’intera collezione di articoli UVA Library StatLab.

Clay Ford
Consulente di ricerca statistica
University of Virginia Library