Régression de Poisson

Introduction à la régression de Poisson

La régression de Poisson sert à modéliser le nombre de fois qu'un événement survient au cours d'une période donnée ou dans une zone donnée. Le critère de jugement est donc un comptage (angl. count) : 0, 1, 2, 3, etc. Pour utiliser correctement cette méthode, il faut comprendre quelques notions centrales ainsi que la loi de probabilité sur laquelle elle repose, à savoir la loi de Poisson. Ce chapitre présente ces notions, la loi de Poisson et le modèle, puis se termine par un exemple complet dans R.

Supposons que l'on étudie le nombre de patients qui consultent chaque jour dans un centre de soins primaires. On souhaite savoir quels facteurs (prédicteurs) influencent le nombre de consultations, et pouvoir prévoir (prédire) combien de consultations sont attendues un jour donné. Chaque consultation est un événement (angl. event), et le nombre de consultations par jour est un comptage. Si l'on rapporte ce nombre à la durée concernée, on obtient un taux (événements par unité de temps) : si 10 patients consultent en une journée, le taux est de 10 patients par jour. Le taux peut aussi se rapporter à une surface, par exemple le nombre de personnes recourant aux soins par kilomètre carré.

La régression de Poisson peut être utilisée dès lors que les événements sont générés par un processus de Poisson, c'est-à-dire un processus qui donne naissance à une loi de Poisson (voir ci-dessous). La nature de l'événement importe peu. La méthode sert à étudier le nombre de décès, d'hospitalisations, d'infections, de clients ou de courriels.

La loi de Poisson

La loi de Poisson a été décrite par le mathématicien français Siméon Denis Poisson. Elle donne la probabilité qu'un certain nombre d'événements survienne dans un intervalle de temps ou d'espace donné, à condition que les événements surviennent indépendamment les uns des autres et à un rythme moyen constant. La loi ne comporte qu'un seul paramètre, λ (lambda), qui est le nombre attendu d'événements dans l'intervalle. La probabilité d'observer exactement $k$ événements est

$$P(Y = k) = \frac{\lambda^{k} e^{-\lambda}}{k!}, \qquad k = 0, 1, 2, \ldots$$

Une propriété caractéristique est que la moyenne et la variance sont toutes deux égales à λ. C'est précisément cette propriété qui, souvent, n'est pas vérifiée dans les données réelles (voir la section sur la surdispersion).

La figure montre la loi de Poisson pour trois valeurs de λ. La loi est discrète : seuls les entiers ont une probabilité. Lorsque λ est petit, l'essentiel de la probabilité se situe près de zéro, et l'absence d'événement est fréquente. Lorsque λ augmente, la distribution se déplace vers la droite, s'élargit et devient de plus en plus symétrique, et l'absence d'événement devient improbable. Le mode de la distribution se situe à λ ou juste en dessous (pour un λ entier, λ et λ − 1 sont également probables), tandis que la moyenne vaut toujours exactement λ.

{
  "$schema": "https://vega.github.io/schema/vega-lite/v5.json",
  "description": "Loi de Poisson pour λ = 1, 4 et 10 : probabilité de chaque nombre d'événements k.",
  "width": "container",
  "height": 280,
  "config": {"view": {"stroke": null}},
  "data": {
    "values": [
      {"k": 0, "lambda": "λ = 1", "p": 0.3679}, {"k": 1, "lambda": "λ = 1", "p": 0.3679}, {"k": 2, "lambda": "λ = 1", "p": 0.1839}, {"k": 3, "lambda": "λ = 1", "p": 0.0613}, {"k": 4, "lambda": "λ = 1", "p": 0.0153}, {"k": 5, "lambda": "λ = 1", "p": 0.0031}, {"k": 6, "lambda": "λ = 1", "p": 0.0005}, {"k": 7, "lambda": "λ = 1", "p": 0.0001},
      {"k": 0, "lambda": "λ = 4", "p": 0.0183}, {"k": 1, "lambda": "λ = 4", "p": 0.0733}, {"k": 2, "lambda": "λ = 4", "p": 0.1465}, {"k": 3, "lambda": "λ = 4", "p": 0.1954}, {"k": 4, "lambda": "λ = 4", "p": 0.1954}, {"k": 5, "lambda": "λ = 4", "p": 0.1563}, {"k": 6, "lambda": "λ = 4", "p": 0.1042}, {"k": 7, "lambda": "λ = 4", "p": 0.0595}, {"k": 8, "lambda": "λ = 4", "p": 0.0298}, {"k": 9, "lambda": "λ = 4", "p": 0.0132}, {"k": 10, "lambda": "λ = 4", "p": 0.0053}, {"k": 11, "lambda": "λ = 4", "p": 0.0019}, {"k": 12, "lambda": "λ = 4", "p": 0.0006},
      {"k": 1, "lambda": "λ = 10", "p": 0.0005}, {"k": 2, "lambda": "λ = 10", "p": 0.0023}, {"k": 3, "lambda": "λ = 10", "p": 0.0076}, {"k": 4, "lambda": "λ = 10", "p": 0.0189}, {"k": 5, "lambda": "λ = 10", "p": 0.0378}, {"k": 6, "lambda": "λ = 10", "p": 0.0631}, {"k": 7, "lambda": "λ = 10", "p": 0.0901}, {"k": 8, "lambda": "λ = 10", "p": 0.1126}, {"k": 9, "lambda": "λ = 10", "p": 0.1251}, {"k": 10, "lambda": "λ = 10", "p": 0.1251}, {"k": 11, "lambda": "λ = 10", "p": 0.1137}, {"k": 12, "lambda": "λ = 10", "p": 0.0948}, {"k": 13, "lambda": "λ = 10", "p": 0.0729}, {"k": 14, "lambda": "λ = 10", "p": 0.0521}, {"k": 15, "lambda": "λ = 10", "p": 0.0347}, {"k": 16, "lambda": "λ = 10", "p": 0.0217}, {"k": 17, "lambda": "λ = 10", "p": 0.0128}, {"k": 18, "lambda": "λ = 10", "p": 0.0071}, {"k": 19, "lambda": "λ = 10", "p": 0.0037}, {"k": 20, "lambda": "λ = 10", "p": 0.0019}
    ]
  },
  "encoding": {
    "x": {"field": "k", "type": "quantitative", "title": "Nombre d'événements (k)", "scale": {"domain": [0, 20]}, "axis": {"tickMinStep": 1, "labelFontSize": 12, "titleFontSize": 13, "grid": false}},
    "y": {"field": "p", "type": "quantitative", "title": "Probabilité P(Y = k)", "axis": {"labelFontSize": 12, "titleFontSize": 13, "format": ".2f"}},
    "color": {"field": "lambda", "type": "nominal", "sort": ["λ = 1", "λ = 4", "λ = 10"], "scale": {"range": ["#0e7490", "#b0442c", "#149e61"]}, "legend": {"title": null, "orient": "top", "labelFontSize": 13}}
  },
  "layer": [
    {"mark": {"type": "line", "strokeWidth": 1, "opacity": 0.35}},
    {"mark": {"type": "point", "filled": true, "size": 60}, "encoding": {"tooltip": [{"field": "lambda", "title": "Nombre attendu"}, {"field": "k", "title": "Nombre d'événements"}, {"field": "p", "title": "Probabilité", "format": ".4f"}]}}
  ]
}

Exemple. Vous travaillez dans un centre de soins primaires où, en moyenne, 10 patients consultent chaque jour, donc λ = 10. Ce nombre varie d'un jour à l'autre. Si les patients arrivent indépendamment les uns des autres, la loi de Poisson permet de calculer la probabilité qu'exactement 5 patients se présentent un jour donné, ou que plus de 5 se présentent. Dans R, on utilise dpois() (probabilité d'exactement k) et ppois() (probabilité de k ou moins) :

# Probabilité qu'exactement 5 patients se présentent lorsque le nombre attendu est 10
dpois(5, lambda = 10)
#> [1] 0.03783327

# Probabilité que plus de 5 patients se présentent
ppois(5, lambda = 10, lower.tail = FALSE)
#> [1] 0.932914

Exactement 5 patients se présentent donc environ un jour sur 26 (3,8 %), tandis que plus de 5 patients se présentent un peu plus de 93 % des jours.

Hypothèses de la loi de Poisson

• Le résultat $k$ est un nombre d'événements dans un intervalle donné et est un entier non négatif (0, 1, 2, 3, etc.).
• Les événements surviennent indépendamment les uns des autres : la survenue d'un événement n'influence pas le moment où survient le suivant.
• Le rythme moyen de survenue des événements est constant dans l'intervalle.
• Deux événements ne peuvent pas survenir exactement au même instant.

La régression de Poisson est un modèle linéaire généralisé

Comme son nom l'indique, la régression de Poisson est un type d'analyse de régression. Elle est utilisée lorsque la variable dépendante (le critère de jugement, Y) est un comptage. La régression de Poisson appartient aux modèles linéaires généralisés (angl. generalized linear models, GLM).

Qu'est-ce qu'un GLM ?

Les modèles linéaires généralisés constituent un cadre qui englobe plusieurs modèles statistiques, dont la régression linéaire classique. En régression linéaire, Y est modélisé directement comme une fonction linéaire des prédicteurs ($X_1, X_2, X_3, \ldots$) plus un terme d'erreur normalement distribué. Cela convient mal aux comptages : un comptage ne peut pas être négatif, et la dispersion augmente souvent avec la moyenne. Dans un GLM, le critère de jugement peut au contraire suivre d'autres lois de la famille dite exponentielle (notamment les lois de Poisson, binomiale, gamma et normale), et la relation entre les prédicteurs et la valeur attendue du critère de jugement passe par une fonction de lien.

La régression de Poisson utilise le logarithme naturel comme fonction de lien. C'est le logarithme du nombre attendu, $\mu = E(Y)$, qui est modélisé comme une fonction linéaire des prédicteurs :

$$\log(\mu_i) = \beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + \cdots + \beta_p X_{pi}$$

Contrairement à la régression linéaire, l'équation ne comporte pas de terme d'erreur distinct : la variation aléatoire autour de $\mu_i$ est décrite par la loi de Poisson. Le logarithme garantit que le nombre attendu est toujours positif, $\mu_i = e^{\beta_0 + \beta_1 X_{1i} + \cdots}$, et rend les effets des prédicteurs multiplicatifs.

Interprétation des coefficients : rapports des taux d'incidence (IRR)

Le modèle étant logarithmique, les coefficients s'interprètent après exponentiation. $e^{\beta}$ est appelé rapport des taux d'incidence (angl. incidence rate ratio, IRR) et indique par quel facteur le nombre attendu (ou le taux) est multiplié lorsque le prédicteur augmente d'une unité, ou pour une catégorie par rapport à la catégorie de référence, les autres prédicteurs restant inchangés. Un IRR = 1 signifie l'absence de différence, un IRR = 2 un taux doublé et un IRR = 0,5 un taux divisé par deux.

Exemples de questions auxquelles répond la régression de Poisson

Exemple 1. Comment la température extérieure et la météo influencent-elles le nombre de rhumes diagnostiqués chaque semaine dans un centre de soins primaires ? Les données peuvent se présenter ainsi (les 13 premières semaines sont affichées) :

La semaine, la température et la météo sont des prédicteurs utilisés pour expliquer et prédire le nombre de rhumes. Le modèle s'écrit

$$\log(\text{nombre attendu de rhumes}_i) = \beta_0 + \beta_1,\text{semaine}_i + \beta_2,\text{température}_i + \beta_3,\text{météo}_i$$

où $\beta_0$ est l'ordonnée à l'origine (le logarithme du nombre attendu lorsque tous les prédicteurs valent zéro ou se trouvent dans leur catégorie de référence) et $\beta_1, \beta_2, \beta_3$ sont les coefficients des prédicteurs. La météo est une variable catégorielle et reçoit en pratique un coefficient par catégorie, à l'exception de la catégorie de référence.

Exemple 2. Combien de clients effectuent un achat dans un magasin chaque jour, et qu'est-ce qui influence ce nombre ?

Exemple 3. Combien de personnes empruntent une ligne de bus donnée pour se rendre au travail, et comment ce nombre est-il influencé par la météo et le jour de la semaine ?

Offset : analyser des taux

Souvent, ce n'est pas le nombre en lui-même qui est intéressant, mais le nombre rapporté à une exposition : temps, surface, distance ou taille de la population. Deux exemples :

• Vous étudiez des décès et connaissez la durée de suivi de chaque individu (le temps d'exposition). C'est alors l'incidence, le nombre d'événements divisé par le temps d'exposition, qui est intéressante. Le terme event rate est souvent utilisé pour désigner le nombre d'événements par unité de temps.
• Vous étudiez le nombre d'espèces animales dans différentes zones géographiques de tailles différentes. C'est alors le nombre d'espèces par unité de surface qui est intéressant.

L'exposition est introduite dans le modèle sous forme d'offset : le logarithme de l'exposition est ajouté dans le membre de droite avec un coefficient fixé à 1,

$$\log(\mu_i) = \log(t_i) + \beta_0 + \beta_1 X_{1i} + \cdots \quad\Longleftrightarrow\quad \log!\left(\frac{\mu_i}{t_i}\right) = \beta_0 + \beta_1 X_{1i} + \cdots$$

Le modèle décrit alors le taux $\mu_i / t_i$, et $e^{\beta}$ devient un rapport de taux. Dans R, l'offset est indiqué avec la fonction offset() dans la formule, ou avec l'argument offset :

# Modèle de Poisson utilisant la surface comme offset
resultat <- glm(events ~ väder + temperatur + vecka + offset(log(area)),
                family = poisson(link = "log"), data = data)

Notez que l'offset doit être indiqué sur l'échelle logarithmique (log(area)), puisque le modèle est logarithmique.

Hypothèses de la régression de Poisson

• Le critère de jugement est un comptage (entiers non négatifs).
• Les observations sont indépendantes les unes des autres.
• Le logarithme du nombre attendu est une fonction linéaire des prédicteurs. Les prédicteurs continus peuvent devoir être transformés ou modélisés de manière flexible (par exemple avec des splines) si la relation n'est pas log-linéaire.
• La variance est égale à la moyenne (équidispersion). Cette hypothèse est souvent violée, voir la section suivante.
• L'exposition (temps, surface, population) est correctement spécifiée comme offset lorsque c'est le taux qui est intéressant.

Surdispersion et excès de zéros

Le modèle de Poisson suppose que la variance est égale à la moyenne. Dans les données réelles, la variance est souvent supérieure à la moyenne, ce que l'on appelle surdispersion. Elle s'explique souvent par l'absence de prédicteurs importants dans le modèle, par la non-indépendance des observations (par exemple des mesures répétées chez un même individu) ou par des différences entre individus que le modèle ne prend pas en compte. Les estimations ponctuelles sont rarement affectées de façon notable, mais les erreurs standard sont trop petites, les intervalles de confiance trop étroits et les valeurs p trop faibles.

La surdispersion peut être évaluée en divisant la statistique du khi-deux de Pearson (ou la déviance résiduelle) par le nombre de degrés de liberté. Un rapport proche de 1 plaide en faveur de l'hypothèse de Poisson. Un rapport nettement supérieur à 1 (règle empirique : au-delà de 1,5) évoque une surdispersion.

# Rapport entre le khi-deux de Pearson et le nombre de degrés de liberté
sum(residuals(model, type = "pearson")^2) / df.residual(model)

En cas de surdispersion, plusieurs options existent :

# 1. Quasi-Poisson : mêmes estimations, erreurs standard multipliées en fonction de la dispersion
model.qp <- glm(y ~ x1 + x2 + offset(log(tid)), family = quasipoisson, data = data)

# 2. Régression binomiale négative : modélise la variation supplémentaire par un paramètre propre
library(MASS)
model.nb <- glm.nb(y ~ x1 + x2 + offset(log(tid)), data = data)

# 3. Erreurs standard robustes (sandwich) pour un modèle de Poisson classique
library(sandwich)
library(lmtest)
coeftest(model, vcov = vcovHC(model, type = "HC0"))

Un autre problème fréquent est l'excès de zéros (angl. zero inflation). Il survient lorsque les zéros proviennent de deux processus différents. Imaginons que l'on étudie le nombre de cigarettes fumées par jour dans une classe où la moitié des élèves fument. Chez les fumeurs, le nombre de cigarettes suit un processus de Poisson, tandis que les non-fumeurs fument toujours zéro cigarette. Les données contiennent alors nettement plus de zéros que ne le prévoit la loi de Poisson. Dans ces situations, un modèle à inflation de zéros est plus adapté, par exemple zeroinfl() du package pscl, qui modélise en deux parties la probabilité d'un zéro « structurel » et le nombre chez les autres :

library(pscl)
# Avant | : modèle du nombre. Après | : modèle des zéros structurels.
model.zip <- zeroinfl(cigaretter ~ alder + kon | foraldrar_roker, data = data, dist = "poisson")
model.zinb <- zeroinfl(cigaretter ~ alder + kon | foraldrar_roker, data = data, dist = "negbin")

Exemple : cancer de la peau à Minneapolis et à Dallas

Les données proviennent d'une étude classique sur les cancers cutanés non mélanocytaires chez les habitants blancs de quatre régions des États-Unis (Scotto J, Kopf AW, Urbach F. Non-melanoma skin cancer among Caucasians in four areas of the United States. Cancer. 1974;34(4):1333-8. PMID 4422100). Deux de ces villes sont utilisées ici, Minneapolis et Dallas. Les données comprennent le nombre de nouveaux cas (cases) et la taille de la population (n) dans huit classes d'âge pour chaque ville. Le jeu de données est petit et peut être lu directement dans R avec read.table(). Les variables u1 à u7 sont des variables indicatrices des classes d'âge et ne sont pas nécessaires dans l'analyse ci-dessous, où la classe d'âge est créée sous forme de facteur.

# Lecture des données directement dans R
nonmel <- read.table(header = TRUE,
text = "
cases city u1 u2 u3 u4 u5 u6 u7 n
1 1 0 1 0 0 0 0 0 0 172675
2 16 0 0 1 0 0 0 0 0 123065
3 30 0 0 0 1 0 0 0 0 96216
4 71 0 0 0 0 1 0 0 0 92051
5 102 0 0 0 0 0 1 0 0 72159
6 130 0 0 0 0 0 0 1 0 54722
7 133 0 0 0 0 0 0 0 1 32185
8 40 0 0 0 0 0 0 0 0 8328
9 4 1 1 0 0 0 0 0 0 181343
10 38 1 0 1 0 0 0 0 0 146207
11 119 1 0 0 1 0 0 0 0 121374
12 221 1 0 0 0 1 0 0 0 111353
13 259 1 0 0 0 0 1 0 0 83004
14 310 1 0 0 0 0 0 1 0 55932
15 226 1 0 0 0 0 0 0 1 29007
16 65 1 0 0 0 0 0 0 0 7583
")

# Création des variables de classe d'âge et de ville.
# age.range est un facteur à 8 niveaux dont la classe d'âge la plus élevée (85+) est le niveau de référence.
# within() évite d'écrire nonmel$ devant chaque variable.
nonmel <- within(nonmel, {
  age.range <- rep(c("15_24", "25_34", "35_44", "45_54", "55_64", "65_74", "75_84", "85+"), 2)
  age.range <- factor(age.range)
  age.range <- relevel(age.range, ref = "85+")
  city <- factor(city, 0:1, c("Minneapolis", "Dallas"))
})

Le modèle

Nous nous intéressons à l'incidence du cancer de la peau, c'est-à-dire au nombre de cas rapporté à la taille de la population, et à la façon dont elle diffère entre les villes et entre les classes d'âge. Le nombre de cas est la variable dépendante et la taille de la population est introduite comme offset (sur l'échelle logarithmique) :

# L'offset peut être indiqué dans la formule
model.1 <- glm(cases ~ city + age.range + offset(log(n)),
               family = poisson(link = "log"), data = nonmel)

# ou avec l'argument offset (même modèle)
model.1 <- glm(cases ~ city + age.range, offset = log(n),
               family = poisson(link = "log"), data = nonmel)

summary(model.1)

Call:
glm(formula = cases ~ city + age.range, family = poisson(link = "log"),
    data = nonmel, offset = log(n))

Deviance Residuals:
     Min        1Q    Median        3Q       Max
-1.50598  -0.48566   0.01639   0.36926   1.24763

Coefficients:
               Estimate Std. Error z value Pr(>|z|)
(Intercept)     -5.4834     0.1037 -52.890  < 2e-16 ***
cityDallas       0.8039     0.0522  15.399  < 2e-16 ***
age.range15_24  -6.1742     0.4577 -13.488  < 2e-16 ***
age.range25_34  -3.5440     0.1675 -21.160  < 2e-16 ***
age.range35_44  -2.3268     0.1275 -18.254  < 2e-16 ***
age.range45_54  -1.5790     0.1138 -13.871  < 2e-16 ***
age.range55_64  -1.0869     0.1109  -9.800  < 2e-16 ***
age.range65_74  -0.5288     0.1086  -4.868 1.13e-06 ***
age.range75_84  -0.1157     0.1109  -1.042    0.297
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 2789.6810  on 15  degrees of freedom
Residual deviance:    8.2585  on  7  degrees of freedom
AIC: 120.5

Number of Fisher Scoring iterations: 4

Les coefficients de summary() sont sur l'échelle logarithmique. L'ordonnée à l'origine (−5,4834) est le logarithme de l'incidence dans le groupe de référence, soit les personnes de 85 ans et plus à Minneapolis : $e^{-5{,}4834} = 0{,}0042$ cas par personne, soit 4,2 cas pour 1 000 personnes pendant la période d'étude.

Rapports des taux d'incidence avec intervalles de confiance

Pour obtenir les rapports des taux d'incidence, on exponentie les coefficients et leurs intervalles de confiance. Avec le package broom et la fonction tidy(), la sortie prend la forme d'un tableau soigné :

library(broom)
# exponentiate = TRUE donne des rapports des taux d'incidence (IRR) au lieu des estimations sur l'échelle logarithmique
tidy(model.1, exponentiate = TRUE, conf.int = TRUE)

# Même chose en R de base, avec des intervalles de confiance de Wald
exp(cbind(IRR = coef(model.1), confint.default(model.1)))

Les rapports des taux d'incidence avec leurs intervalles de confiance à 95 % (Wald) sont :

Comment interpréter le résultat ?

• Ville. Le rapport des taux d'incidence pour Dallas est de 2,23. L'incidence du cancer de la peau est donc 2,23 fois plus élevée à Dallas qu'à Minneapolis (123 % de plus) à structure d'âge identique. Il s'agit d'une différence relative de taux, et non de 2,23 cas de plus par personne. Minneapolis est la catégorie de référence et n'apparaît donc pas dans le tableau. Cette différence est cohérente avec le rayonnement UV plus intense à Dallas, située nettement plus au sud.
• Âge. La classe la plus âgée (85+) est la catégorie de référence. Toutes les classes d'âge plus jeunes ont des rapports des taux d'incidence inférieurs à 1, c'est-à-dire une incidence plus faible que les plus âgés, et le rapport augmente progressivement avec l'âge. Les personnes de 15 à 24 ans ont une incidence d'environ 0,2 % de celle des personnes de 85 ans et plus. Pour les 75–84 ans, l'intervalle de confiance inclut 1 (p = 0,30) : l'incidence ne diffère pas de façon certaine de celle de la classe la plus âgée.

Les rapports des taux d'incidence peuvent être représentés dans un forest plot. L'échelle est logarithmique, de sorte que les rapports inférieurs et supérieurs à 1 se trouvent à la même distance de la ligne d'absence de différence. Avec le package sjPlot, un tel graphique s'obtient avec plot_model(model.1) ; le graphique ci-dessous présente les mêmes estimations.

{
  "$schema": "https://vega.github.io/schema/vega-lite/v5.json",
  "description": "Rapports des taux d'incidence (IRR) avec intervalles de confiance à 95 % issus du modèle de Poisson pour le cancer de la peau : ville et classe d'âge.",
  "width": "container",
  "height": 300,
  "config": {"view": {"stroke": null}},
  "data": {
    "values": [
      {"term": "Dallas (par rapport à Minneapolis)", "irr": 2.2342, "lo": 2.017, "hi": 2.4749},
      {"term": "75–84 ans", "irr": 0.8907, "lo": 0.7167, "hi": 1.107},
      {"term": "65–74 ans", "irr": 0.5893, "lo": 0.4763, "hi": 0.7291},
      {"term": "55–64 ans", "irr": 0.3373, "lo": 0.2714, "hi": 0.4191},
      {"term": "45–54 ans", "irr": 0.2062, "lo": 0.165, "hi": 0.2577},
      {"term": "35–44 ans", "irr": 0.0976, "lo": 0.076, "hi": 0.1253},
      {"term": "25–34 ans", "irr": 0.0289, "lo": 0.0208, "hi": 0.0401},
      {"term": "15–24 ans", "irr": 0.0021, "lo": 0.0008, "hi": 0.0051}
    ]
  },
  "encoding": {
    "y": {"field": "term", "type": "nominal", "sort": null, "title": null, "axis": {"labelFontSize": 13, "labelLimit": 260}},
    "x": {"type": "quantitative", "scale": {"type": "log", "domain": [0.0005, 5]}, "title": "Rapport des taux d'incidence (IRR), échelle logarithmique", "axis": {"labelFontSize": 12, "titleFontSize": 13, "values": [0.001, 0.01, 0.1, 0.5, 1, 2, 5]}}
  },
  "layer": [
    {"mark": {"type": "rule", "strokeDash": [4, 4], "color": "#8a8a96"}, "encoding": {"x": {"datum": 1}, "y": null}},
    {"mark": {"type": "rule", "strokeWidth": 2, "color": "#0e7490"}, "encoding": {"x": {"field": "lo"}, "x2": {"field": "hi"}}},
    {"mark": {"type": "point", "filled": true, "size": 80, "color": "#0e7490"}, "encoding": {"x": {"field": "irr"}, "tooltip": [{"field": "term", "title": "Comparaison"}, {"field": "irr", "title": "IRR", "format": ".3f"}, {"field": "lo", "title": "Borne inférieure", "format": ".3f"}, {"field": "hi", "title": "Borne supérieure", "format": ".3f"}]}}
  ]
}

Le modèle est-il adapté ?

Deux vérifications doivent toujours être effectuées après une régression de Poisson.

Adéquation (goodness of fit). La déviance résiduelle (8,26 pour 7 degrés de liberté) peut être comparée à une loi du khi-deux. Une valeur p élevée signifie qu'il n'y a pas d'indice d'une mauvaise adéquation du modèle :

pchisq(deviance(model.1), df.residual(model.1), lower.tail = FALSE)
#> [1] 0.3103526

Surdispersion. Le rapport entre la déviance résiduelle et les degrés de liberté est de 8,26 / 7 = 1,18, donc proche de 1. L'hypothèse de Poisson d'égalité entre moyenne et variance paraît raisonnable, et il n'est pas nécessaire de recourir à un modèle quasi-Poisson ou à une régression binomiale négative dans cet exemple.

Prédictions

Avec predict(), le modèle peut prédire l'incidence pour n'importe quelle combinaison de ville et de classe d'âge. Le modèle étant logarithmique, on exponentie la prédiction (ou on indique type = "response"). En fixant la population n à 1, on obtient le nombre attendu de cas par personne, et avec n = 1 000, le nombre pour 1 000 personnes. Un nombre pour 1 000 personnes est plus facile à interpréter qu'un petit nombre décimal.

# Nombre attendu de cas par personne, Minneapolis, 85 ans et plus
predict(model.1, newdata = data.frame(city = "Minneapolis", age.range = "85+", n = 1),
        type = "response")
#>           1
#> 0.004155093

# Pour 1 000 personnes
predict(model.1, newdata = data.frame(city = "Minneapolis", age.range = "85+", n = 1000),
        type = "response")
#>        1
#> 4.155093

Dans la classe d'âge la plus élevée à Minneapolis, on attend donc 4,2 cas pour 1 000 personnes pendant la période d'étude. Le chiffre correspondant à Dallas est de 4,2 × 2,23 = 9,3 pour 1 000.

Les prédictions pour toutes les lignes du jeu de données s'obtiennent de la même manière :

# Nouveaux jeux de données avec la même ville et la même classe d'âge que l'original
newdat1 <- nonmel
newdat2 <- nonmel

# La population est fixée à 1 pour obtenir le nombre de cas par personne
newdat1$n <- 1
nonmel$pred.cases.per.one <- predict(model.1, newdat1, type = "response")

# La population est fixée à 1 000 pour obtenir le nombre de cas pour 1 000 personnes
newdat2$n <- 1000
nonmel$pred.cases.per.thousand <- predict(model.1, newdat2, type = "response")

# Nombre attendu de cas dans la population réelle
nonmel$pred.cases <- predict(model.1, type = "response")

nonmel[, c("city", "age.range", "cases", "n", "pred.cases.per.thousand", "pred.cases")]

          city age.range cases      n pred.cases.per.thousand pred.cases
1  Minneapolis     15_24     1 172675              0.00865295   1.494148
2  Minneapolis     25_34    16 123065              0.12007322  14.776810
3  Minneapolis     35_44    30  96216              0.40555928  39.021292
4  Minneapolis     45_54    71  92051              0.85670789  78.860818
5  Minneapolis     55_64   102  72159              1.40134673 101.119779
6  Minneapolis     65_74   130  54722              2.44870028 133.997777
7  Minneapolis     75_84   133  32185              3.70128210 119.125764
8  Minneapolis       85+    40   8328              4.15509270  34.603612
9       Dallas     15_24     4 181343              0.01933271   3.505852
10      Dallas     25_34    38 146207              0.26827163  39.223190
11      Dallas     35_44   119 121374              0.90611423 109.978708
12      Dallas     45_54   221 111353              1.91408567 213.139182
13      Dallas     55_64   259  83004              3.13093612 259.880221
14      Dallas     65_74   310  55932              5.47096874 306.002223
15      Dallas     75_84   226  29007              8.26952928 239.874236
16      Dallas       85+    65   7583              9.28344824  70.396388

La colonne pred.cases peut être comparée directement au nombre de cas observés (cases). La concordance est bonne dans tous les groupes, ce qui est cohérent avec la faible déviance résiduelle.

Résumé

• La régression de Poisson est utilisée lorsque le critère de jugement est un nombre d'événements, et avec un offset lorsque le critère de jugement est un taux (événements par unité de temps, de surface ou de population).
• Le modèle est un GLM avec une fonction de lien logarithmique. Les coefficients exponentiés sont des rapports des taux d'incidence (IRR), qui s'interprètent de façon multiplicative.
• La loi de Poisson suppose que la variance est égale à la moyenne. Vérifiez toujours l'adéquation et la dispersion du modèle.
• En cas de surdispersion, on utilise un modèle quasi-Poisson, une régression binomiale négative ou des erreurs standard robustes. En cas d'excès de zéros, on utilise des modèles à inflation de zéros.