Lorsqu’il s’agit de modéliser des comptages (c’est-à-dire des nombres entiers supérieurs ou égaux à 0), nous commençons souvent par une régression de Poisson. Il s’agit d’un modèle linéaire généralisé où l’on suppose qu’une réponse a une distribution de Poisson conditionnelle à une somme pondérée de prédicteurs. Par exemple, nous pourrions modéliser le nombre de commotions cérébrales documentées chez les quarterbacks de la NFL comme une fonction des sessions de jeu et du nombre total d’années d’expérience de sa ligne offensive. Cependant un inconvénient potentiel de la régression de Poisson est qu’elle peut ne pas décrire avec précision la variabilité des comptages.

Une distribution de Poisson est paramétrée par \(\lambda\), qui se trouve être à la fois sa moyenne et sa variance. Bien que pratique à retenir, ce n’est pas souvent réaliste. Une distribution de nombres aura généralement une variance qui n’est pas égale à sa moyenne. Lorsque cela se produit avec des données que nous supposons (ou espérons) être distribuées selon la loi de Poisson, nous disons qu’il y a sous-dispersion ou surdispersion, selon que la variance est plus petite ou plus grande que la moyenne. Effectuer une régression de Poisson sur des données de comptage qui présentent ce comportement aboutit à un modèle qui ne s’ajuste pas bien.

Une approche qui aborde ce problème est la régression binomiale négative. La distribution binomiale négative, comme la distribution de Poisson, décrit les probabilités d’occurrence des nombres entiers supérieurs ou égaux à 0. Contrairement à la distribution de Poisson, la variance et la moyenne ne sont pas équivalentes. Cela suggère qu’elle pourrait servir d’approximation utile pour modéliser des nombres dont la variabilité est différente de sa moyenne. La variance d’une distribution binomiale négative est une fonction de sa moyenne et possède un paramètre supplémentaire, k, appelé paramètre de dispersion. Disons que notre comptage est une variable aléatoire Y issue d’une distribution binomiale négative, alors la variance de Y est

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

A mesure que le paramètre de dispersion devient de plus en plus grand, la variance converge vers la même valeur que la moyenne, et la binomiale négative se transforme en une distribution de Poisson.

Pour illustrer la distribution binomiale négative, travaillons avec des données tirées du livre, Categorical Data Analysis, d’Alan Agresti (2002). Les données sont présentées dans le tableau 13.6 de la section 13.4.3. Les données proviennent d’une enquête menée auprès de 1308 personnes auxquelles on a demandé combien de victimes d’homicide elles connaissaient. Les variables sont resp, le nombre de victimes que la personne interrogée connaît, et race, la race de la personne interrogée (noire ou blanche). La race permet-elle d’expliquer le nombre de victimes d’homicide qu’une personne connaît ? Les données doivent d’abord être saisies dans 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)

Avant de passer à la modélisation, explorons les données. Tout d’abord, nous remarquons que la plupart des répondants sont blancs:

> table(race)racewhite black 1149 159 

Les Noirs ont un nombre moyen plus élevé que les Blancs :

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

Pour chaque race, la variance de l’échantillon est à peu près le double de la moyenne. Il semble que nous ayons une surdispersion.

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

Enfin, nous examinons la distribution des comptages par race.

> 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

On passe à l’ajustement du modèle. Nous essayons d’abord la régression de Poisson en utilisant la fonction glm() et montrons une partie de la sortie de résumé.

> # 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 race est très significative. Il semble que les Noirs soient beaucoup plus susceptibles de connaître quelqu’un qui a été victime d’un homicide. Mais que signifie le coefficient 1,73 ? Dans ce modèle simple avec un prédicteur dichotomique, il s’agit de la différence entre les logarithmes des effectifs attendus. Si nous exponentions le coefficient, nous obtenons un rapport de moyennes d’échantillon :

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

En fait, si nous faisons une prédiction avec ce modèle et que nous exponentions les résultats, nous obtenons les moyennes d’échantillon :

> 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 

Cela dit que le compte de victimes connues pour les blancs est distribué comme un Poisson avec une moyenne et une variance égale à 0.09, tandis que le nombre de victimes connues pour les Noirs est distribué comme une Poisson avec une moyenne et une variance égales à 0,52. Nous savons déjà, grâce à notre analyse exploratoire, que les variances observées étaient beaucoup plus grandes, nous ne devrions donc pas être trop satisfaits des variances estimées du modèle. Si nous examinons les comptes ajustés, nous verrons encore plus de preuves du manque d’ajustement :

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

Au-dessus, nous avons d’abord enregistré les moyennes prédites dans un objet appelé fmeans. Nous avons ensuite généré des comptes ajustés en utilisant la fonction dpois avec les moyennes estimées pour prédire la probabilité d’obtenir 0 à 6. Nous avons ensuite multiplié ces probabilités par le nombre de répondants pour obtenir des comptes ajustés. Enfin, nous avons combiné le tout dans un cadre de données afin de comparer facilement les valeurs observées et ajustées. Deux des choses les plus spectaculaires à noter sont que nous sous-adaptons les comptes de 0 et sur-adaptons les comptes de 1.

Nous pouvons utiliser un rootogramme pour visualiser l’ajustement d’un modèle de régression de comptage. La fonction rootogram() du package countreg rend cela facile.

> countreg::rootogram(pGLM)

La ligne courbe rouge est l’ajustement théorique de Poisson. “Suspendu” à chaque point de la ligne rouge est une barre, dont la hauteur représente la différence entre les comptages attendus et observés. Une barre suspendue en dessous de 0 indique un sous-ajustement. Une barre supérieure à 0 indique un ajustement excessif. Les comptages ont été transformés avec une transformation de racine carrée pour éviter que les plus petits comptages soient obscurcis et écrasés par les plus grands. Nous voyons beaucoup de sous-adaptation pour les comptes 2 et plus et une sur-adaptation massive pour le compte 1.

Essayons maintenant d’ajuster un modèle binomial négatif. Nous avons remarqué que la variabilité des comptes était plus grande pour les deux races. Il semblerait que la distribution binomiale négative s’approcherait mieux de la distribution des comptages.

Pour ajuster un modèle binomial négatif dans R, nous nous tournons vers la fonction glm.nb() du package MASS (un package qui est installé avec R). Encore une fois, nous ne montrons qu’une partie de la sortie sommaire:

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

Premièrement, remarquez que les coefficients sont les mêmes que précédemment. Une fois de plus, nous pouvons exponentialiser le coefficient de la course pour obtenir un rapport des moyennes de l’échantillon et faire des prédictions pour obtenir les moyennes de l’échantillon 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 

Mais remarquez que l’erreur standard pour le coefficient de race est plus grande, indiquant plus d’incertitude dans notre estimation (0,24 contre 0,15). Ceci est logique étant donné la variabilité observée dans nos comptages. Remarquez également l’estimation de Thêta. C’est notre paramètre de dispersion. Nous pouvons y accéder directement depuis notre objet modèle comme suit :

> nbGLM$theta 0.2023119

Et nous pouvons l’utiliser pour obtenir des variances estimées pour les comptages :

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

Elles sont beaucoup plus proches des variances observées que celles données par le modèle de Poisson.

Une fois encore, nous visualisons l’ajustement en utilisant un rootogramme:

> countreg::rootogram(nbGLM)

Cela semble bien meilleur que le rootogramme du modèle de Poisson. Il y a un léger sous-ajustement/surajustement pour les comptes 1 à 3, mais sinon, cela semble assez bon.

Pour mieux comprendre notre modèle binomial négatif, utilisons ses paramètres pour simuler des données et comparons les données simulées aux données observées. Ci-dessous, nous chargeons le paquet magrittr pour accéder à l’opérateur %>% qui nous permet de “chaîner” des fonctions.

D’abord, nous tabulons les comptes et créons un diagramme à barres pour les participants blancs et noirs, respectivement. Ensuite, nous utilisons les paramètres du modèle pour simuler les données d’une distribution binomiale négative. Les deux paramètres sont mu et la taille (c’est-à-dire le paramètre de dispersion). Remarquez que nous utilisons la fonction coef() pour extraire les coefficients appropriés pour chaque race. Pour les blancs, c’est juste l’intercept, pour les noirs, c’est l’intercept et la pente (nous les additionnons donc). Ensuite, nous exponentions pour convertir l’échelle logarithmique à l’échelle originale. Nous simulons le même nombre d’observations que nous avons dans nos données originales.

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)

Les données simulées sont très similaires aux données observées, ce qui nous donne à nouveau confiance dans le choix de la régression binomiale négative pour modéliser ces données. Ces graphiques démontrent également la nature conditionnelle de notre modèle. La distribution binomiale négative des comptes dépend, ou est conditionnée, par la race. Chaque race a une moyenne différente mais un paramètre de dispersion commun.

• 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

Pour toute question ou clarification concernant cet article, contactez le StatLab de la bibliothèque de l’UVA : [email protected]

Voir la collection complète des articles du StatLab de la bibliothèque de l’UVA.

Clay Ford
Consultant en recherche statistique
Bibliothèque de l’Université de Virginie
.