---
title: "TD 1 : Modèle de Poisson (R)"
subtitle: "Années universitaires 2025–2026 et 2026–2027 • Parcours Économie de la santé & Développement durable"
description: "Régression log-linéaire : Tabac & cancer du poumon"
author: "Pierre Beaucoral"
format:
html:
theme: [cosmo, ../edu.scss]
slide-number: true
progress: true
hash: true
center: false
transition: slide
toc: true
number-sections: true
code-fold: true
code-tools: true
code-overflow: scroll
smaller: true
self-contained: true
pdf: default
editor:
markdown:
wrap: 72
---
> **Contexte**: On étudie d’anciennes données reliant **tabagisme** et
> **décès par cancer du poumon**. Variables : `age` (classes),
> `smoking status` (4 classes), `population` (centaines de milliers),
> `deaths` (décès annuels).
# Préparation
## Objectifs de ce TD
- Importer et préparer un tableau **comptages + exposition**
(population à risque).
- Ajuster un **GLM Poisson** avec *offset* (log-exposition).
- Évaluer l’ajustement : **déviance** (vs modèle saturé) &
**Pearson**.
- Comparer des modèles via **tests de rapport de vraisemblance (LR)**.
- Interpréter en **ratios de taux d'incidence (IRR)** et produire des
**comptes attendus**.
## Packages
```{r}
#| label: setup
# installer si nécessaire : install.packages(c("readxl","dplyr","tidyr","janitor","ggplot2","broom","gt","performance","DescTools"))
library(readxl)
library(dplyr)
library(tidyr)
library(janitor)
library(ggplot2)
library(broom)
library(gt)
library(performance) # check_overdispersion
library(DescTools) # PChisq for GOF (si besoin)
# Formatage des nombres dans le texte (virgule décimale, p-values)
fr <- function(x, d = 2) formatC(x, format = "f", digits = d, decimal.mark = ",")
fmt_p <- function(p) if (p < 0.001) "< 0,001" else paste("=", fr(p, 3))
```
::: {.callout-note icon="📘" title="Rappel de cours : Régression de Poisson"}
La **régression de Poisson** est un modèle linéaire généralisé (GLM)
adapté aux données de **comptage** (ex. nombre de décès, d’accidents, de
visites).
**Formulation :**
- Variable dépendante : un comptage $Y_i \in \{0,1,2,\dots\}$. - Loi
supposée : $Y_i \sim \text{Poisson}(\mu_i)$ avec
$\mathbb{E}[Y_i] = \mu_i$.
- Lien log : $$
\log(\mu_i) = \beta_0 + \beta_1 X_{1i} + \dots + \beta_p X_{pi} + \log(\text{exposition}_i)
$$ où l’offset $\log(\text{exposition}_i)$ tient compte de la taille
de la population ou du temps d’observation.
**Pourquoi utiliser ce modèle ?**
- Les données sont des comptages positifs (non négatifs).
- La variance est proportionnelle à la moyenne ($\text{Var}(Y)=\mu$).
- On cherche à modéliser des **taux d’incidence** (décès/population,
accidents/temps).
**Interprétation des coefficients :**
- Les $\beta_j$ s’interprètent via l’**Incidence Rate Ratio (IRR)** :\
$$
IRR_j = e^{\beta_j}
$$ → $IRR_j > 1$ : le taux est plus élevé que la référence.\
→ $IRR_j < 1$ : le taux est plus faible.
**Diagnostics courants :**
- Tests de qualité d’ajustement (déviance, Pearson).
- Vérification de la sur-dispersion (si $\text{Var}(Y) > \mu$,
préférer quasi-Poisson ou binomiale négative).
**Extensions :**
- Modèle de Poisson avec offset (exposition).
- Quasi-Poisson pour corriger la variance.
- Binomiale négative pour sur-dispersion forte.
:::
# Import et manipulation des données
## Importer les données `smoking_dat.xlsx`
[Télécharger `smoking_dat.xlsx`](data/smoking_dat.xlsx){download="smoking_dat.xlsx"}
> Dans l’énoncé, les données à importer sont `smoking_dat.xlsx` et les
> variables `age`, `smoking status`, `population`, `deaths`.
::: callout-note
Dictionnaire des variables :\
• age: en classes (40-44, 45-49, 50-54, 55-59, 60-64, 65-69, 70-74,
75-79, 80+).\
• smoking status: 4 classes (ne fume pas / fume le cigare ou la pipe /
fume la cigarette et le cigare ou la pipe ; fume seulement la
cigarette)\
• population: en centaine de milliers de personnes\
• deaths: comptage des décès par cancer du poumon en un an.
:::
```{r}
#| label: import
# Chemin suggéré : placez le fichier dans data/smoking_dat.xlsx
# Si vous avez un CSV, remplacez read_excel par read.csv(...)
data_path <- "data/smoking_dat.xlsx"
df <- read_excel(data_path) |>
clean_names()
# Harmonisation de noms
# On s'attend à des colonnes: age (classes), smoking_status (4 classes), population, deaths
df <- df |>
rename(
age = matches("^age$|^age_class|^agecat"),
smoking_status = matches("^smoking|^smoke"),
population = matches("^pop|^population"),
deaths = matches("^deaths|^dead")
)
glimpse(df)
```
## Coder les deux variables enregistrées en texte avec des chiffres
> En Stata on ferait `encode` + `i.variable`.\
> En R, il suffit de déclarer les variables comme `factor`.
```{r}
#| label: factors
df <- df |>
mutate(
age = factor(age, ordered = FALSE),
smoking_status = factor(smoking_status, ordered = FALSE)
)
# Vérification des niveaux
levels(df$age); levels(df$smoking_status)
```
::: callout-note
**Rappel** : Les facteurs indiquent à R qu’il s’agit de variables
qualitatives. Chaque modalité sera transformée en **variable
indicatrice** (dummy) dans la régression.
:::
### Unité d'exposition
L’énoncé précise que `population` est en **centaines de milliers**. Pour
une interprétation plus intuitive, on peut ramener l’exposition à
l’unité **personne** (facultatif) :
```{r}
#| label: exposure-scale
# Ici, on transforme 'population' en nombre de personnes si besoin.
# Exemple: si population = 2.3 signifie 2.3 * 100 000 personnes :
expo_personnes <- TRUE
scale_factor <- 1e5
df <- df |>
mutate(
exposure = if (expo_personnes) population * scale_factor else population
)
summary(df$exposure)
```
# Estimations
## Appliquer un premier modèle de régression log-linéaire
Modèle Poisson **log-linéaire** avec effets de `smoking_status` et
`age`, et *offset* `log(exposure)` : c’est l’équivalent de Stata
`poisson deaths i.smokecod i.agecod, exposure(pop)`.
::: callout-tip
**Pourquoi un modèle de Poisson ?**\
Les données sont des **comptages** (nombre de décès).\
Le modèle de Poisson relie l’**espérance** de ces comptages à des
variables explicatives par une fonction de lien log :\
$$
\log(\mathbb{E}[Y]) = X\beta
$$\
Cette structure garantit que la prédiction est **positive** et que la
variance est proportionnelle à la moyenne (hypothèse de Poisson).
:::
Nous voulons expliquer le nombre de décès par l’âge et le statut
tabagique, en tenant compte de l’exposition.
```{r}
#| label: poisson1
mod1 <- glm(deaths ~ smoking_status + age + offset(log(exposure)),
family = poisson(link = "log"), data = df)
summary(mod1)
```
::: callout-important
L’argument `offset(log(exposure))` ajoute le **log de l’exposition**
avec un coefficient fixé à 1. Cela revient à modéliser un **taux** de
décès (décès / population).
**Quelle interprétation?**
On peut utiliser $\exp(\beta_i)$ pour retrouver **l'Incidence Rate Ratio
(IRR).** On le lit comme le coefficient multiplicateur de l'occurence de
$Y$ (ici le décès) par rapport à la catégorie de référence.
:::
## Calculer la déviance du modèle ajusté. DEV1
```{r}
#| label: deviance1
dev1 <- deviance(mod1) # Deviance du modèle vs saturé
df_dev1 <- df.residual(mod1)
c(dev1 = dev1, df = df_dev1, p_value = pchisq(dev1, df_dev1, lower.tail = FALSE))
```
**Réponse.** La déviance du modèle ajusté vaut **DEV1 =
`r fr(dev1)`** pour **`r df_dev1` degrés de liberté**. La déviance
mesure l’écart entre notre modèle et le modèle **saturé** (un paramètre
par cellule, qui reproduit exactement les données) :
$DEV = 2\,[\ell(\text{saturé}) - \ell(\text{modèle})]$. Une déviance
proche de son nombre de degrés de liberté (ici `r fr(dev1)` pour
`r df_dev1`) indique déjà un bon ajustement ; on le teste formellement
plus bas.
## Quel est l’effet de l’âge sur la probabilité de décès par cancer du poumon ?
```{r}
#| label: irr-age
irr1 <- exp(coef(mod1))
ci1 <- exp(confint.default(mod1)) # IC 95 % de Wald, exponentiés
tibble(age = levels(df$age)[-1],
IRR = irr1[paste0("age", age)],
ic_bas = ci1[paste0("age", age), 1],
ic_haut = ci1[paste0("age", age), 2]) |>
gt() |>
fmt_number(columns = -age, decimals = 2) |>
tab_header(title = "IRR de l'âge (référence : 40-44 ans)")
```
**Réponse.** L’âge est le déterminant **le plus fort** de la mortalité
par cancer du poumon. Tous les coefficients d’âge sont positifs et
hautement significatifs (p \< 0,001), et les IRR **croissent de façon
monotone** avec la classe d’âge. À statut tabagique donné, et par
rapport aux 40-44 ans, le taux de décès est multiplié par :
- **`r fr(irr1["age45-59"])`** chez les 45-49 ans ;
- **`r fr(irr1["age60-64"])`** chez les 60-64 ans ;
- **`r fr(irr1["age80+"], 1)`** chez les 80 ans et plus (IC 95 % :
`r fr(ci1["age80+", 1], 1)` à `r fr(ci1["age80+", 2], 1)`).
Autrement dit, d’une classe quinquennale à la suivante, le taux augmente
d’environ 30 à 50 % (+74 % entre les deux premières classes). Ce gradient reflète l’accumulation de l’exposition
et le temps de latence du cancer du poumon.
::: callout-note
**À noter.** Le fichier nomme la 2^e^ classe `45-59` : c’est une coquille
pour **45-49**, puisque la classe suivante est 50-54. Par ailleurs,
l’ordonnée à l’origine dépend de l’unité choisie pour l’exposition (la constante
vaut le log du taux de la cellule de référence), alors que les IRR n’en
dépendent pas. Ce sont eux qu’on interprète.
:::
## Sauvegarde du modèle (utile pour LR tests)
```{r}
#| label: store1
# En R, on garde l'objet en mémoire (mod1). Pas besoin d'"estimates store".
# On peut aussi l'ajouter à une liste si on veut gérer plusieurs modèles :
models <- list(mod1 = mod1)
```
## Tester l’ajustement de ce modèle
Pour cela nous allons réaliser **deux tests** :
- **Deviance GOF** (modèle ajusté vs **saturé**).
- **Pearson GOF** (comparaison effectifs attendus vs observés). Les
**degrés de liberté** correspondent ici au nombre de cellules moins
le nombre de paramètres estimés (y compris l’ordonnée à l’origine). L’énoncé
suggère **24 df** pour chacun de ces tests (voir justification plus
bas).
```{r}
#| label: gof
# Deviance test (déjà calculé ci-dessus)
dev_stat <- deviance(mod1)
dev_df <- df.residual(mod1)
dev_p <- pchisq(dev_stat, dev_df, lower.tail = FALSE)
# Pearson GOF (somme (y - mu)^2 / mu) et chi2 approx avec mêmes df résiduels)
mu_hat <- fitted(mod1) # comptages attendus
y_obs <- df$deaths
pearson <- sum((y_obs - mu_hat)^2 / mu_hat)
pearson_p <- pchisq(pearson, dev_df, lower.tail = FALSE)
tibble(
test = c("Deviance GOF", "Pearson GOF"),
statistic = c(dev_stat, pearson),
df = c(dev_df, dev_df),
p_value = c(dev_p, pearson_p)
) |>
gt()
```
**Réponse.** Les deux tests concluent au **bon ajustement** du modèle :
- **Déviance** : $G^2 =$ `r fr(dev_stat)`, `r dev_df` df, p
`r fmt_p(dev_p)` ;
- **Pearson** : $X^2 =$ `r fr(pearson)`, `r dev_df` df, p
`r fmt_p(pearson_p)`.
Les hypothèses testées sont :
- H0 : le modèle reproduit les données aussi bien que le modèle saturé
(les écarts observés − attendus sont dus au hasard) ;
- H1 : le modèle est mal spécifié (il manque des effets, par exemple
une interaction âge × tabac).
Les deux p-values sont largement supérieures à 0,05 : on **ne rejette
pas H0**. Un modèle **additif** (effets de l’âge et du tabac qui se
multiplient sur le taux) suffit donc à décrire les 36 cellules, et il
n’est pas nécessaire d’ajouter une interaction âge × tabac.
::: callout-important
Attention au sens du test : ici, c’est une p-value **élevée** qui
rassure. Ne pas rejeter H0 veut dire que les données ne contredisent pas
le modèle ; cela ne prouve pas qu’il est « vrai ».
:::
## Justifiez pourquoi ces tests sont effectués avec le même **df**.
```{r}
#| label: df-check
c(cellules = nrow(df), parametres = length(coef(mod1)),
df = nrow(df) - length(coef(mod1)),
cellules_mu_inf_5 = sum(mu_hat < 5), mu_min = min(mu_hat))
```
**Réponse.** Les deux statistiques comparent **la même chose** : les
comptages observés $y_i$ et les comptages ajustés $\hat\mu_i$ du même
modèle, sur les mêmes cellules. La déviance le fait par un rapport de
vraisemblances,
$G^2 = 2\sum_i y_i \log(y_i/\hat\mu_i) - (y_i - \hat\mu_i)$, et Pearson
par des écarts au carré normalisés,
$X^2 = \sum_i (y_i - \hat\mu_i)^2/\hat\mu_i$. Ce sont deux
approximations asymptotiques équivalentes du même test (modèle contre
modèle saturé), qui suivent sous H0 la même loi
$\chi^2(N - p)$ :
- $N =$ `r nrow(df)` cellules (9 classes d’âge × 4 statuts
tabagiques), soit les paramètres du modèle saturé ;
- $p =$ `r length(coef(mod1))` paramètres estimés (1 constante + 3
indicatrices tabac + 8 indicatrices âge) ;
- d’où **df = `r nrow(df)` − `r length(coef(mod1))` =
`r nrow(df) - length(coef(mod1))`**, ce qui correspond aux 24 df de
l’énoncé.
Condition de validité : l’approximation χ² demande des effectifs
attendus suffisants (règle usuelle $\hat\mu_i \geq 5$). Ici, seules
`r sum(mu_hat < 5)` cellules sur `r nrow(df)` ont $\hat\mu_i < 5$
(minimum `r fr(min(mu_hat))`) : ce sont les fumeurs de cigare/pipe de
40-44 et 45-49 ans, des cellules à faible effectif. On peut
donc se fier à l’approximation, et les deux statistiques sont d’ailleurs
très proches.
## Ajuster un modèle sans la variable smoke, et effectuer un test de rapport de vraisemblance entre ce nouveau modèle et celui précédemment sauvegardé
On ajuste un modèle **sans tabac** et on compare à `mod1` par **LR
test**. En Stata : `lrtest`. En R : `anova(mod0, mod1, test="Chisq")`.
```{r}
#| label: poisson0
mod0 <- glm(deaths ~ age + offset(log(exposure)),
family = poisson, data = df)
lr01 <- anova(mod0, mod1, test = "Chisq")
lr01
```
**Réponse.** Le test du rapport de vraisemblance compare deux modèles
**emboîtés** : `mod0` (âge seul) est un cas particulier de `mod1` (âge +
tabac), obtenu en imposant que les 3 coefficients du tabac soient nuls.
La statistique est la différence des déviances,
$LR = DEV_0 - DEV_1 =$ `r fr(lr01$Deviance[2])`, qui suit sous H0 un
$\chi^2(3)$. La p-value est `r sub("^= ", "", fmt_p(lr01$"Pr(>Chi)"[2]))` : on
**rejette nettement H0**. En plus, `mod0` s’ajuste mal à lui seul
(déviance `r fr(deviance(mod0))` pour `r df.residual(mod0)` df), alors
que `mod1` s’ajuste bien.
# Conclure sur l’impact de l’usage du tabac sur la probabilité de décès par cancer du poumon.
**Réponse.** Le statut tabagique a un effet **significatif** sur la
mortalité par cancer du poumon, même à âge donné. Le test LR indique
**que** le tabac compte, mais pas **dans quel sens** : pour cela, il faut
lire les coefficients.
::: callout-warning
**Catégorie de référence.** R prend la première modalité par ordre
alphabétique, ici `cigarPipeOnly` (cigare ou pipe seulement), et non
les non-fumeurs. Les IRR ci-dessous se lisent donc **par rapport aux
fumeurs de cigare/pipe**. Pour prendre les non-fumeurs comme
référence : `df$smoking_status <- relevel(df$smoking_status, ref = "no")`.
:::
À âge égal, par rapport aux fumeurs de cigare/pipe seulement :
- **cigarette seulement** : taux multiplié par
**`r fr(irr1["smoking_statuscigarretteOnly"])`** (+
`r fr(100 * (irr1["smoking_statuscigarretteOnly"] - 1), 0)` %,
IC 95 % `r fr(ci1["smoking_statuscigarretteOnly", 1])` à
`r fr(ci1["smoking_statuscigarretteOnly", 2])`) ;
- **cigarette + cigare/pipe** : ×
**`r fr(irr1["smoking_statuscigarrettePlus"])`** (IC 95 %
`r fr(ci1["smoking_statuscigarrettePlus", 1])` à
`r fr(ci1["smoking_statuscigarrettePlus", 2])`) ;
- **non-fumeurs** : × **`r fr(irr1["smoking_statusno"])`**, différence
**non significative** (p
`r fmt_p(summary(mod1)$coefficients["smoking_statusno", 4])`).
Par rapport aux **non-fumeurs**, il suffit de faire le rapport des IRR :
les fumeurs de cigarettes seules ont un taux
`r fr(irr1["smoking_statuscigarretteOnly"] / irr1["smoking_statusno"])`
fois plus élevé (+
`r fr(100 * (irr1["smoking_statuscigarretteOnly"] / irr1["smoking_statusno"] - 1), 0)`
%). **Conclusion** : fumer des cigarettes augmente nettement le risque de
décès par cancer du poumon, alors que le cigare ou la pipe seuls ne se
distinguent pas des non-fumeurs dans ces données.
## Construire une nouvelle variable qui prend la valeur 1 si l’individu fume des cigarettes, 0 s’il n’en fume pas.
Créer une variable `cigarette_user` égale à **1 si l’individu fume des
cigarettes**, **0 sinon** : l’énoncé demande de distinguer le type de
produit et de concentrer l’attention sur la **cigarette**.
```{r}
#| label: binary-cig
# Adaptez le motif à vos libellés (ex.: "fume seulement la cigarette", "cigarette + cigare/pipe", etc.)
# On classera 1 si l'étiquette contient "cigarette", 0 sinon.
df <- df |>
mutate(
cigarette_user = as.integer(grepl("cigarrette", tolower(as.character(smoking_status)))))
table(df$cigarette_user, df$smoking_status)
```
**Réponse.** Le tableau croisé sert de vérification : `cigarette_user`
vaut 1 pour les deux modalités qui contiennent la cigarette
(`cigarretteOnly`, `cigarrettePlus`) et 0 pour `no` et `cigarPipeOnly`,
soit 18 cellules de chaque côté. Le motif `"cigarrette"` reprend
l’orthographe du fichier (avec deux *r*) : un motif `"cigarette"` ne
trouverait aucune ligne, et tout le monde serait codé 0.
## Ajuster un troisième modèle avec effet de l’âge et de cette variable d’usage de la cigarette.
```{r}
#| label: poisson2
mod2 <- glm(deaths ~ cigarette_user + age + offset(log(exposure)),
family = poisson, data = df)
summary(mod2)
anova(mod2, test = "Chisq")
irr2 <- exp(coef(mod2))["cigarette_user"]
ci2 <- exp(confint.default(mod2))["cigarette_user", ]
```
**Réponse.** À âge égal, les fumeurs de cigarettes (seules ou avec
cigare/pipe) ont un taux de décès multiplié par **`r fr(irr2)`** par
rapport aux autres (IC 95 % `r fr(ci2[1])` à `r fr(ci2[2])`,
p \< 0,001). L’effet est significatif, mais ce modèle **s’ajuste mal** :
sa déviance résiduelle vaut `r fr(deviance(mod2))` pour
`r df.residual(mod2)` df (p
`r fmt_p(pchisq(deviance(mod2), df.residual(mod2), lower.tail = FALSE))`),
contre `r fr(dev1)` pour `r df_dev1` df pour `mod1`. En regroupant les
modalités, on perd de l’information. Le test suivant le confirme.
# Comparer ce modèle avec le modèle initial par un test de rapport de vraisemblance
On compare le **modèle multiniveaux par type de tabac** (`mod1`) avec le
**modèle binaire cigarette vs non** (`mod2`).
```{r}
#| label: lr-1-2
lr21 <- anova(mod2, mod1, test = "Chisq")
lr21
```
**Réponse.** `mod2` est emboîté dans `mod1` : il revient à imposer
$\beta_{\text{cigOnly}} = \beta_{\text{cigPlus}}$ et
$\beta_{\text{no}} = \beta_{\text{cigarPipe}} = 0$, soit **2
contraintes** (`mod1` a 12 paramètres, `mod2` en a 10). La statistique
vaut $LR =$ `r fr(lr21$Deviance[2])` pour `r lr21$Df[2]` df, avec
p `r fmt_p(lr21$"Pr(>Chi)"[2])` : on **rejette** le modèle simplifié.
La variable binaire « cigarette oui/non » est **trop grossière** : les
deux groupes qu’elle fusionne n’ont pas le même risque.
# Le type de produit fumé semble-t-il influencer la probabilité de décès par cancer du poumon ?
En Poisson log-linéaire : **IRR = exp(coef)**. On reporte aussi des **IC
95%** exponentiés.
```{r}
#| label: irr
tidy(mod1, conf.int = TRUE, exponentiate = TRUE) |>
select(term, estimate, conf.low, conf.high, p.value) |>
arrange(desc(estimate)) |>
gt() |>
tab_header(title = "IRR (modèle 1) : exp(coef) avec IC 95%")
```
**Réponse.** Oui, le type de produit compte. Trois constats (référence :
cigare/pipe seulement, à âge égal) :
1. **La cigarette est le produit le plus nocif.** Les fumeurs de
cigarettes seules ont le taux le plus élevé
(IRR = `r fr(irr1["smoking_statuscigarretteOnly"])`, soit
`r fr(irr1["smoking_statuscigarretteOnly"] / irr1["smoking_statusno"])`
fois celui des non-fumeurs).
2. **Cigare et pipe seuls ne se distinguent pas des non-fumeurs.**
L’IRR des non-fumeurs par rapport au cigare/pipe
(`r fr(irr1["smoking_statusno"])`) n’est pas significativement
différent de 1, et son IC contient 1.
3. **Les fumeurs mixtes (cigarette + cigare/pipe)** ont un risque
intermédiaire (IRR = `r fr(irr1["smoking_statuscigarrettePlus"])`),
significativement plus faible que celui des fumeurs de cigarettes
seules : les IC ne se recouvrent pas. Une explication possible :
ces fumeurs remplacent une partie de leurs cigarettes par le cigare
ou la pipe, et consomment donc moins de cigarettes.
C’est pour cela que le test précédent rejette `mod2` : il met dans le
même groupe les fumeurs de cigarettes seules et les fumeurs mixtes, dont
les risques diffèrent, et dans un autre groupe les non-fumeurs et les
fumeurs de cigare/pipe, dont les risques sont proches. Le second
regroupement est acceptable, le premier ne l’est pas.
::: callout-note
Ces données sont **observationnelles** : les IRR décrivent des
associations ajustées sur l’âge. On ne contrôle ni l’intensité de la
consommation (nombre de cigarettes par jour), ni sa durée, ni d’autres
facteurs de confusion.
:::
# Extensions
## Comptes attendus & tableau Observé vs Attendu
```{r}
#| label: nij
df_preds <- df |>
mutate(
y_obs = deaths,
mu_hat = fitted(mod1)
)
df_preds |>
select(age, smoking_status, exposure, y_obs, mu_hat) |>
arrange(age, smoking_status) |>
gt() |>
fmt_number(columns = c(exposure, y_obs, mu_hat), decimals = 2)
res_pearson <- residuals(mod1, type = "pearson")
```
**Lecture.** Les comptages attendus $\hat\mu_{ij}$ sont très proches des
comptages observés dans toutes les cellules. Les résidus de Pearson,
$(y - \hat\mu)/\sqrt{\hat\mu}$, sont tous compris entre
`r fr(min(res_pearson))` et `r fr(max(res_pearson))`. Seule la cellule
des non-fumeurs de 60-64 ans atteint tout juste −2 (un peu moins de
décès que prévu), ce qui est attendu par hasard sur 36 cellules : pas
de cellule franchement mal prédite. Le modèle additif âge + tabac ne présente donc pas d’écart
systématique, ce qui est cohérent avec les tests d’ajustement.
## Vérification d’éventuelle sur-dispersion
> Le GLM Poisson suppose `Var(Y) = E[Y]`. Si `Var(Y) >> E[Y]`, la
> **sur-dispersion** peut invalider les tests usuels (SE sous-estimés).
```{r}
#| label: overdisp
check_overdispersion(mod1)
```
**Lecture.** Le rapport de dispersion vaut
$X^2/df =$ `r fr(pearson / dev_df)`, proche de 1 et même légèrement
inférieur (p `r fmt_p(pearson_p)`) : **pas de sur-dispersion**.
L’hypothèse $\text{Var}(Y) = \mathbb{E}[Y]$ est compatible avec les
données, donc les erreurs standards et les tests du modèle de Poisson
sont fiables. Il n’est pas nécessaire de passer à un modèle
quasi-Poisson ou binomial négatif.
> En cas de sur-dispersion marquée, envisagez **Quasi-Poisson**
> (`family = quasipoisson`) ou **Négative Binomiale** (`MASS::glm.nb`)
> et comparez l’ajustement.
## Graphiques (facultatifs)
```{r}
#| label: plot-fitted
shapes_9 <- c(16, 17, 15, 3, 7, 8, 0, 1, 2) # à ta convenance
ggplot(df_preds, aes(mu_hat, y_obs, color = smoking_status, shape = age)) +
geom_point(size = 2) +
geom_abline(intercept = 0, slope = 1, linetype = 2) +
scale_shape_manual(values = shapes_9) +
labs(x = "Comptes attendus (mod1)", y = "Comptes observés",
title = "Observé vs Attendu : GLM Poisson (mod1)")
```
## Bilan
- Le modèle de Poisson avec offset permet d’estimer des **taux de
décès** en fonction du tabagisme et de l’âge.
- Le modèle additif âge + tabac (`mod1`) **s’ajuste bien** (déviance
`r fr(dev1)` pour `r df_dev1` df, p `r fmt_p(dev_p)`), sans
sur-dispersion.
- L’**âge** est le facteur dominant : le taux est multiplié par
`r fr(irr1["age80+"], 1)` entre 40-44 ans et 80 ans et plus.
- Le **tabac** a un effet significatif (LR = `r fr(lr01$Deviance[2])`,
3 df), qui tient à la **cigarette** : à âge égal, les fumeurs de
cigarettes seules ont un taux
`r fr(irr1["smoking_statuscigarretteOnly"] / irr1["smoking_statusno"])`
fois plus élevé que les non-fumeurs, alors que cigare et pipe seuls
ne s’en distinguent pas.
- La variable binaire « cigarette oui/non » est **trop grossière**
(rejetée par le test LR contre `mod1`) : le type de produit compte.
- Les **IRR** offrent une interprétation intuitive : *par rapport à la
catégorie de référence, combien de fois le taux de décès est-il
multiplié*. Il faut toujours vérifier quelle est cette référence.
> **Conseil pratique** : en recherche appliquée, vérifiez toujours la
> sur-dispersion et documentez les hypothèses de variance (Poisson vs
> quasi-Poisson vs binomiale négative).
