Site web de GR

4.A Modèles pour les données de comptage surdispersées

Nous utilisons les données de Long (1990) sur le nombre de publications produites par les biochimistes titulaires d’un doctorat pour illustrer l’application des modèles de Poisson, de Poisson surdispersé, binomial négatif et de Poisson zéro-gonflé.

Les variables de l’ensemble de données sont

articles dans les trois dernières années du doctorat.

.

prestige du Ph.D. programme

articles par mentor dans les trois dernières années

Variable	Description
art
fem	codé un pour les femmes
mar	codé un si marié
kid5	nombre d’enfants de moins de six ans
phd
ment

Ces données ont également été analysées par Long et Freese (2001), et sont disponibles sur le site Web de 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 

Le nombre moyen d’articles est de 1,69 et la variance est de 3,71, un peu plus du double de la moyenne. Les données sont surdispersées, mais bien sûr, nous n’avons pas encore considéré de covariables.

Un modèle de Poisson

Envisageons le modèle utilisé par Long et Freese(2001), un modèle additif simple utilisant les cinq prédicteurs.

> 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
et la déviance et le chi-carré de Pearson sont tous deux dans les 1600.
 4.A.1. Variation extra-poisson
Nous supposons maintenant que la variance est proportionnelle plutôt qu'égale à la moyenne, et estimons le paramètre d'échelle φ en divisant le chi carré de Pearson par sa d.f.:
> phi round(c(phi,sqrt(phi)),4) 1.8290 1.3524
Nous voyons que la variance est environ 83% plus grande que la moyenne. Cela signifie que nous devrions ajuster les erreurs standard en multipliant par 1,35, la racine carrée de 1,83.
R peut faire ce calcul pour nous si nous utilisons la quasipoissonfamille:
> 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
Une approche alternative consiste à ajuster un modèle de Poisson et à utiliser l'estimateur robuste ou sandwich des erreurs standard. Cela donne généralement des résultats très similaires au modèle de Poisson surdispersé.
 4.A.2 Régression binomiale négative
Nous ajustons maintenant un modèle binomial négatif avec les mêmes prédicteurs.Pour ce faire, nous avons besoin de la fonction glm.nb() du paquetMASS.
> 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 est la précision de l'effet multiplicativerandomique, et correspond à 1/σ2 dans thenotes. L'estimation correspond à une variance estimée de0,44 et est hautement significative.
Pour tester la signification de ce paramètre, vous pouvez penser à calculerdeux fois la différence de log-vraisemblance entre ce modèle et le modèle de Poisson, 180,2, et à la traiter comme un chi carré avec un d.f.. Les asymptotiques habituelles ne s'appliquent pas, cependant, parce que l'hypothèse nulle est sur une frontière de l'espace des paramètres.
Certains travaux montrent qu'une meilleure approximation consiste à traiter la statistique comme un mélange 50:50 de zéro et d'un chi-carré avec un d.f. Alternativement, traiter la statistique comme un chi-carré avec un d.f. donne un test conservateur. Quoi qu'il en soit, nous avons une preuve écrasante de surdispersion.
Pour tester les hypothèses sur les coefficients de régression, nous pouvons utiliser soit des tests de Wald, soit des tests de rapport de vraisemblance, qui sont possibles parce que nous avons fait des hypothèses de distribution complètes.
Il y a aussi une famille negative.binomial pour glm et cela peut être utilisé à condition que les paramètrestheta soient donnés. (Ceci est basé sur le résultat que la binomiale négative est dans la famille glm pour une variance fixe). Voici une vérification rapide:
> 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)

Nous pouvons également calculer des quantiles en utilisantqgamma(p, shape, rate = 1, scale = 1/rate, lower.tail = TRUE). En particulier, lorsque l'effet aléatoire a une variance vles quartiles sont 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
Les biochimistes situés à Q1 de la distribution de l'hétérogénéité non observée publient 49% moins d'articles que prévu à partir de leurs caractéristiques observées, tandis que ceux situés à la médiane publient 14% de moins et ceux situés à Q3 publient 33% de plus que prévu.
Comparaison des estimations et des erreurs standard
Comparons les estimations des paramètres et les erreurs standard sous les modèles de Poisson, de Poisson surdispersé et binomial négatif:
> 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
Les estimations binomiales négatives ne sont pas très différentes de celles basées sur le modèle de Poisson, et les deux ensembles conduiraient aux mêmes conclusions. 
En regardant les erreurs standard, nous voyons que les deux approches de la surdispersion conduisent à des erreurs standard estimées très similaires,et que la régression de Poisson ordinaire sous-estime les erreurs standard.
Bonne qualité d'ajustement
Nous pouvons évaluer la bonne qualité d'ajustement du modèle binomial négatif en utilisant la déviance
> deviance(mnbg) 1004.281
Le modèle binomial négatif s'ajuste mieux que le modèle de Poisson, mais a toujours une déviance supérieure à la valeur critique de cinq pour cent de980,25.
La fonction de variance
Les modèles de Poisson et binomial négatif surdispersés ont des fonctions de variance différentes. Une façon de vérifier lequel peut être le plus approprié est de créer des groupes basés sur le prédicteur linéaire, de calculer la moyenne et la variance pour chaque groupe, et enfin de tracer la relation moyenne-variance.
Voici des groupes basés sur le prédicteur linéaire binomial négatif, créés à l'aide de cut avec des ruptures aux 5(5)95 percentilespour produire vingt groupes de taille approximativement égale. Notez que predict() calcule un prédicteur linéaire,sauf indication contraire. Pour prédire dans l'échelle de la réponseajoutez l'option 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 

Le graphique trace la moyenne en fonction de la variance et superpose les courbes correspondant au modèle de Poisson surdispersé,où la variance est φμ, et au modèle binomial négatif, où la variance est μ(1+μσ2).
La fonction de variance de Poisson fait un assez bon travail pour la majeure partie des données, mais ne parvient pas à capturer les variances élevées des universitaires les plus productifs. La fonction de variance binomiale négative n'est pas trop différente mais, étant quadratique, peut monter plus rapidement et fait un meilleur travail à l'extrémité supérieure. Nous concluons que le modèle binomial négatif fournit une meilleure description des données que le modèle de Poisson surdispersé.
 4.A.3 Modèles de Poisson à zéro gonflé
 Une occurrence fréquente avec les données de comptage est un excès de zéros par rapport à ce qui est attendu sous un modèle de Poisson. C'est en fait un problème avec nos données : 
> zobs zpoi c(obs=mean(zobs), poi=mean(zpoi)) obs poi 0.3005464 0.2092071 
Nous voyons que 30,0% des scientifiques de l'échantillon n'ont publié aucun article au cours des trois dernières années de leur doctorat, mais le modèle de Poisson prédit que seuls 20,9% n'auraient aucune publication. Il est clair que le modèle sous-estime la probabilité des comptes nuls.
Une façon de modéliser ce type de situation est de supposer que les données proviennent d'un mélange de deux populations, une où les comptes sont toujours nuls, et une autre où le compte a une distribution de Poisson avec une moyenne μ. Dans ce modèle, les comptes nuls peuvent provenir de l'une ou l'autre des populations, tandis que les comptes positifs ne proviennent que de la seconde.
Dans le contexte des publications des biochimistes titulaires d'un doctorat, on peut imaginer que certains avaient en tête des emplois où les publications ne seraient pas importantes, tandis que d'autres visaient des emplois universitaires où l'on s'attendait à un record de publications. Les membres du premier groupe publieraient zéro article, tandis que les membres du second groupe publieraient 0,1,2,..., un décompte que l'on peut supposer avoir une distribution de Poisson.
La distribution du résultat peut alors être modélisée en termes de deux paramètres, π la probabilité de " toujours zéro ", et μ, le nombre moyen de publications pour ceux qui ne font pas partie du groupe " toujours zéro ". Une façon naturelle d'introduire des covariables est de modéliser le logit de la probabilité π de toujours zéro et le log de la moyenne μ pour ceux qui ne sont pas dans la classe toujours zéro.
Ce type de modèle peut être ajusté dans R en utilisant la fonction zeroinfl()du paquet pscl. La formule du modèle peut être spécifiée comme d'habitude si les mêmes variables doivent être incluses dans les deux équations. Sinon, on peut fournir deux ensembles de prédicteursséparés par une barre verticale, tapez ?zeroinfl pour en savoir plus.
Voici un modèle gonflé à zéro avec toutes les covariables dans les deux équations:
> 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
En regardant l'équation gonflée, nous voyons que le seul prédicteur significatif d'être dans la classe 'toujours zéro' est le nombre d'articles publiés par le mentor, chaque article du mentor étant associé à 12,6% de chances inférieures de ne jamais publier.
En regardant l'équation pour le nombre moyen d'articles parmi ceux qui ne sont pas dans la classe toujours zéro, nous trouvons des désavantages significatifs pour les femmes et les scientifiques avec des enfants de moins de cinq ans, et un grand effet positif du nombre de publications par le mentor, chaque article étant associé à une augmentation de 1,8% du nombre attendu de publications.
Pour vérifier que le modèle résout le problème de l'excès de zéros, nous prédisons π et μ, et calculons la probabilité combinée de ne jamais publier. Il existe des options dans la predict()fonction appelées "zero" et "count" pour les obtenir. Il y a aussi une option "prob" pour calculer la densité prédite, mais c'est superflu car nous voulons seulement la probabilité de zéro.
> pr mu zip mean(zip) 0.2985679
Le modèle résout donc le problème de l'excès de zéros, en prédisant que 29,9% des biochimistes ne publieront aucun article, ce qui est beaucoup plus proche de la valeur observée de 30,0%.
Il y a aussi un modèle binomial négatif gonflé à zéro, qui utilise un binôme négatif pour le compte dans la classe 'pas toujours zéro'. Ce modèle peut être ajusté en utilisant zeroinfl() avec le paramètre dist="negbin". Les liens alternatifs pour l'équation de gonflement incluent le probit,qui peut être spécifié en utilisant link="probit". 
Comparaison des modèles avec AIC
Comme il se trouve, pour ces données, la binomiale négative résout aussi le problème ! Voici la probabilité d'articles nuls dans la binomiale négative. Nous procédons à partir des premiers principes mais on pourrait aussi utiliser la densité binomiale négative intégrée
> munb theta znb # also dnbinom(0, mu=munb, size=theta)> mean(znb) 0.3035957
Le modèle prédit que 30,4% des biochimistes ne publieraient aucun article dans les trois dernières années de leur doctorat, très proche de la valeur observée de 30,0%.
Pour choisir entre le binôme négatif et les modèles gonflés à zéro, nous devons recourir à d'autres critères. Une façon très simple de comparer des modèles avec différents nombres de paramètres est de calculer le critère d'information d'Akaike (AIC), que nous définissons comme
AIC = -2logL + 2p
où p est le nombre de paramètres du modèle. Le premier terme est essentiellement la déviance et le second une pénalité pour le nombre de paramètres, et il peut être calculé pour la plupart des modèles en utilisant la fonction AIC. Je l'obtiens également "à la main" afin de pouvoir vérifier le calcul. Pour nos données
> 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
 Pour cet ensemble de données, le modèle binomial négatif est clairement gagnant en termes de parcimonie et de qualité d'ajustement. 
Les autres critères de diagnostic que nous pourrions examiner sont la distribution marginale des effectifs prédits et observés et les fonctions de variance. 
Modèles à troncature zéro et à haies
Les autres modèles que nous n'avons pas couverts sont le Poisson à troncature zéro et la binomiale négative, conçus pour les données qui ne comprennent pas de zéros. Un exemple courant est la durée de séjour dans un hôpital, qui est d'au moins un jour. Une approche judicieuse consiste à ajuster un modèle de Poisson ou binomial négatif qui exclut le zéro et modifie l'échelle des autres probabilités pour que leur somme soit égale à un. Il faut être prudent dans l'interprétation de ces modèles car μ n'est pas le résultat attendu, mais la moyenne d'une distribution sous-jacente qui inclut les zéros. Ces modèles sont implémentés dans la fonction vglm() du paquet VGAM, en utilisant les familles pospoisson et posnegbinomial.
Une approche alternative à l'excès (ou à la pénurie) de zéros est d'utiliser un processus à deux étapes, avec un modèle logit pour distinguer les comptes zéro et positifs, puis un modèle de Poisson ou binomial négatif tronqué à zéro pour les comptes positifs. Dans notre exemple, nous pourrions utiliser un modèle logit pour différencier ceux qui publient de ceux qui ne le font pas, puis un modèle de Poisson tronqué ou binomial négatif pour le nombre d'articles de ceux qui en publient au moins un. Ces modèles sont souvent appelés modèles à seuil. Ils peuvent être ajustés dans R en utilisant les modèles logit et Poisson tronqué zéro ou binomial négatif séparés et en ajoutant simplement les log-vraisemblances, ou en utilisant la fonction hurdle() dans le paquet pscl.
Comparant les modèles à haies et les modèles à dégonflement zéro, je trouve que la distinction entre zéro et un ou plus est plus claire avec les modèles à haies, mais l'interprétation de la moyenne est plus claire avec les modèles à dégonflement zéro.
Le document JSS suivant présente une discussion utile de tous ces modèles : Modèles de régression pour les données de comptage en R.
.