GR’s Website

4.A Modelli per dati di conteggio sovradispersi

Utilizziamo i dati di Long (1990) sul numero di pubblicazioni prodotte dai biochimici di dottorato per illustrare l’applicazione dei modelli Poisson, Poisson sovradispersi, binomiale negativo e Poisson inflazionato zero.

Le variabili del dataset sono

mar

articoli per mentore negli ultimi tre anni

Variabile	Descrizione
art	articoli negli ultimi tre anni di dottorato.
fem	codificato uno per le femmine
codificato uno se sposato
kid5	numero di figli sotto i sei anni
phd	prestigio del Ph.D. programma
ment

Questi dati sono stati analizzati anche da Long e Freese (2001), e sono disponibili sul sito di 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 

Il numero medio di articoli è 1,69 e la varianza è 3,71, un po’ più del doppio della media. I dati sono sovradimensionati, ma naturalmente non abbiamo ancora considerato nessuna covariata.

Un modello di Poisson

Adattiamo il modello usato da Long e Freese(2001), un modello additivo semplice che usa tutti e cinque i predittori.

> 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
e la devianza e il chi-quadrato di Pearson sono entrambi nel 1600.
 4.A.1. Variazione extra-Poisson
Assumiamo ora che la varianza sia proporzionale piuttosto che uguale alla media, e stimiamo il parametro di scala φ dividendo il chi-quadrato di Pearson per la sua d.f.:
> phi round(c(phi,sqrt(phi)),4) 1.8290 1.3524
Vediamo che la varianza è circa 83% più grande della media. Questo significa che dovremmo aggiustare gli errori standard moltiplicando per 1,35, la radice quadrata di 1,83.
R può fare questo calcolo per noi se usiamo la quasipoissonfamiglia:
> 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
Un approccio alternativo è quello di adattare un modello Poisson e usare lo stimatore robusto o sandwich degli errori standard. Questo di solito dà risultati molto simili al modello over-dispersed Poisson.
 4.A.2 Regressione binomiale negativa
Adattiamo ora un modello binomiale negativo con gli stessi predittori.Per fare questo abbiamo bisogno della funzione glm.nb() nel pacchettoMASS.
> 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  è la precisione dell'effetto multiplicativerandom, e corrisponde a 1/σ2 in poi. La stima corrisponde a una varianza stimata di 0,44 ed è altamente significativa.
Per testare la significatività di questo parametro si può pensare di calcolare due volte la differenza di log-likelihood tra questo modello e quello di Poisson, 180,2, e trattarla come un chi-quadrato con una d.f. La solita asintotica non si applica, tuttavia, perché l'ipotesi nulla è su un confine dello spazio dei parametri.
C'è un lavoro che mostra che una migliore approssimazione è trattare la statistica come una miscela 50:50 di zero e un chi-quadrato con una d.f. In alternativa, trattare la statistica come un chi-quadrato con una d.f. dà un test conservativo. In entrambi i casi, abbiamo una prova schiacciante di iperdispersione.
Per testare le ipotesi sui coefficienti di regressione possiamo usare sia i test di Wald che i test del rapporto di verosimiglianza, che sono possibili perché abbiamo fatto delle ipotesi distributive complete.
C'è anche una famiglia negative.binomial per glm e questa può essere usata a condizione che i parametritheta siano dati. (Questo si basa sul risultato che il binomio negativo è nella famiglia glm per la varianza fissa). Ecco un rapido controllo:
> 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)

Possiamo anche calcolare i quantili usandoqgamma(p, shape, rate = 1, scale = 1/rate, lower.tail = TRUE). In particolare, quando l'effetto casuale ha varianza vi quartili sono 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
I biochimici al Q1 della distribuzione dell'eterogeneità non osservata pubblicano il 49% in meno di articoli rispetto alle loro caratteristiche osservate, mentre quelli alla mediana pubblicano il 14% in meno e quelli al Q3 il 33% in più del previsto.
Confrontiamo le stime e gli errori standard
Confrontiamo le stime dei parametri e gli errori standard secondo i modelli Poisson, Poisson iperdisperso e binomiale negativo:
> 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
Le stime del binomiale negativo non sono molto diverse da quelle basate sul modello Poisson, ed entrambe le serie porterebbero alle stesse conclusioni. 
Guardando gli errori standard, vediamo che entrambi gli approcci alla sovradispersione portano a errori standard stimati molto simili, e che la regressione ordinaria di Poisson sottostima gli errori standard.
Bontà dell'adattamento
Possiamo valutare la bontà dell'adattamento del modello binomiale negativo usando la devianza
> deviance(mnbg) 1004.281
Il modello binomiale negativo si adatta meglio di quello di Poisson, ma ha ancora una devianza superiore al valore critico del cinque per cento di 980.25.
La funzione di varianza
I modelli Poisson e binomiale negativo iperdisperso hanno funzioni di varianza diverse. Un modo per verificare quale può essere più appropriato è quello di creare gruppi basati sul predittore lineare, calcolare la media e la varianza per ogni gruppo, e infine tracciare la relazione media-varianza.
Ecco i gruppi basati sul predittore lineare binomiale negativo, creati usando cut con interruzioni al 5(5)95 percentile per produrre venti gruppi di dimensioni approssimativamente uguali. Notate che predict() calcola un predittore lineare, se non diversamente specificato. Per predire nella scala della risposta aggiungete l'opzione 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 

Il grafico traccia la media contro la varianza e sovrappone le curve corrispondenti al modello di Poisson iperdisperso, dove la varianza è φμ, e il modello binomiale negativo, dove la varianza è μ(1+μσ2).
La funzione di varianza di Poisson fa un buon lavoro per la maggior parte dei dati, ma non riesce a catturare le alte varianze degli studiosi più produttivi. La funzione di varianza binomiale negativa non è troppo diversa ma, essendo quadratica, può aumentare più velocemente e fa un lavoro migliore nella fascia alta. Concludiamo che il modello binomiale negativo fornisce una migliore descrizione dei dati rispetto al modello di Poisson eccessivamente disperso.
 4.A.3 Modelli di Poisson inflazionati da zero
 Un evento frequente con i dati di conteggio è un eccesso di zeri rispetto a quanto ci si aspetta in un modello Poisson. Questo è effettivamente un problema con i nostri dati: 
> zobs zpoi c(obs=mean(zobs), poi=mean(zpoi)) obs poi 0.3005464 0.2092071 
Vediamo che il 30,0% degli scienziati del campione non ha pubblicato alcun articolo negli ultimi tre anni del loro dottorato, ma il modello di Poisson prevede che solo il 20,9% non abbia pubblicazioni. Chiaramente il modello sottostima la probabilità di zero conteggi.
Un modo per modellare questo tipo di situazione è assumere che i dati provengano da una miscela di due popolazioni, una dove il conteggio è sempre zero, e un'altra dove il conteggio ha una distribuzione di Poisson con media μ. In questo modello i conteggi zero possono provenire da entrambe le popolazioni, mentre i conteggi positivi provengono solo dalla seconda.
Nel contesto delle pubblicazioni dei biochimici di dottorato, possiamo immaginare che alcuni avessero in mente lavori in cui le pubblicazioni non sarebbero state importanti, mentre altri miravano a lavori accademici in cui ci si aspettava un record di pubblicazioni. I membri del primo gruppo avrebbero pubblicato zero articoli, mentre i membri del secondo gruppo avrebbero pubblicato 0,1,2,..., un conteggio che può essere assunto avere una distribuzione di Poisson.
La distribuzione del risultato può quindi essere modellata in termini di due parametri, π la probabilità di 'sempre zero', e μ, il numero medio di pubblicazioni per quelli non nel gruppo 'sempre zero'. Un modo naturale per introdurre le covariate è quello di modellare il logit della probabilità π di sempre zero e il log della media μ per coloro che non sono nella classe sempre zero.
Questo tipo di modello può essere adattato in R usando la funzione zeroinfl() del pacchetto pscl. La formula del modello può essere specificata come al solito se le stesse variabili devono essere incluse in entrambe le equazioni. Altrimenti si possono fornire due serie di predittori separati da una barra verticale, digitando ?zeroinfl per saperne di più.
Ecco un modello gonfiato a zero con tutte le covariate in entrambe le equazioni:
> 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
Guardando l'equazione gonfiata vediamo che l'unico predittore significativo di essere nella classe 'sempre zero' è il numero di articoli pubblicati dal mentore, con ogni articolo del mentore associato al 12,6% di probabilità di non pubblicare mai.
Guardando l'equazione per il numero medio di articoli tra coloro che non sono nella classe sempre zero, troviamo svantaggi significativi per le femmine e gli scienziati con figli sotto i cinque anni, e un grande effetto positivo del numero di pubblicazioni da parte del mentore, con ogni articolo associato a un aumento dell'1,8% nel numero previsto di pubblicazioni.
Per verificare che il modello risolva il problema dell'eccesso di zero prevediamo π e μ, e calcoliamo la probabilità combinata di non pubblicare. Ci sono opzioni nella predict()funzione chiamate "zero" e "count" per ottenere queste. C'è anche un'opzione "prob" per calcolare la densità prevista, ma questo è eccessivo dato che vogliamo solo la probabilità di zero.
> pr mu zip mean(zip) 0.2985679
Quindi il modello risolve il problema dell'eccesso di zero, prevedendo che il 29,9% dei biochimici non pubblicherà alcun articolo, molto più vicino al valore osservato del 30,0%.
C'è anche un modello binomiale negativo inflazionato di zero, che usa un binomio negativo per il conteggio nella classe 'non sempre zero'. Questo modello può essere adattato usando zeroinfl() con il parametro dist="negbin". Collegamenti alternativi per l'equazione di inflazione includono il probit, che può essere specificato usando link="probit". 
Confronto del modello con AIC
Come succede, anche per questi dati il binomiale negativo risolve il problema! Ecco la probabilità di zero articoli nel binomio negativo. Procediamo dai primi principi, ma si potrebbe anche usare la densità binomiale negativa incorporata
> munb theta znb # also dnbinom(0, mu=munb, size=theta)> mean(znb) 0.3035957
Il modello predice che il 30,4% dei biochimici non pubblicherebbe alcun articolo negli ultimi tre anni di dottorato, molto vicino al valore osservato del 30,0%.
Per scegliere tra i modelli binomiale negativo e zero gonfiato dobbiamo ricorrere ad altri criteri. Un modo molto semplice per confrontare modelli con diversi numeri di parametri è quello di calcolare il Criterio di Informazione di Akaike (AIC), che definiamo come
AIC = -2logL + 2p
dove p è il numero di parametri nel modello. Il primo termine è essenzialmente la devianza e il secondo una penalità per il numero di parametri, e può essere calcolato per la maggior parte dei modelli usando la funzione AIC. Lo ottengo anche "a mano" in modo da poter verificare il calcolo. Per i nostri dati
> 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
 Per questo set di dati il modello binomiale negativo è un chiaro vincitore in termini di parsimonia e bontà di adattamento. 
Altri criteri diagnostici che potremmo guardare sono la distribuzione marginale dei conteggi previsti e osservati e le funzioni di varianza. 
Modelli troncati a zero e a ostacoli
Altri modelli che non abbiamo coperto sono il Poisson troncato a zero e il binomiale negativo, progettati per dati che non includono zeri. Un esempio comune è la durata della permanenza in un ospedale, che è di almeno un giorno. Un approccio ragionevole è quello di adattare un modello Poisson o binomiale negativo che esclude lo zero e ridimensiona le altre probabilità per sommarle a uno. Si dovrebbe fare attenzione ad interpretare questi modelli perché μ non è il risultato atteso, ma la media di una distribuzione sottostante che include gli zeri. Questi modelli sono implementati nella funzione vglm() del pacchetto VGAM, usando le famiglie pospoisson e posnegbinomial.
Un approccio alternativo all'eccesso (o alla scarsità) di zeri è quello di usare un processo a due stadi, con un modello logit per distinguere tra conteggi zero e positivi e poi un modello Poisson o binomiale negativo troncato a zero per i conteggi positivi. Nel nostro esempio potremmo usare un modello logit per differenziare quelli che pubblicano da quelli che non lo fanno, e poi un modello Poisson troncato o binomiale negativo per il numero di articoli di quelli che ne pubblicano almeno uno. Questi modelli sono spesso chiamati modelli a ostacoli. Possono essere applicati in R usando i modelli logit e zero-truncatedPoisson o binomiale negativo separati e semplicemente aggiungendo le probabilità logiche, o usando la funzione hurdle() del pacchetto pscl.
Confrontando i modelli hurdle e zero-inflated, trovo che la distinzione tra zero e uno o più sia più chiara con i modelli hurdle, ma l'interpretazione della media è più chiara con i modelli zero-inflated.
Il seguente documento JSS ha una discussione utile di tutti questi modelli: Regression Models for Count Data in R.
.