Calculer des indices standardisés de mortalité avec R et les visualiser avec ggplot2

, par Joël Girès

Je présente dans cet article une manière de calculer des indices comparatifs standardisés avec R, en prenant l’exemple de l’indice standardisé de mortalité (Standardized Mortality Ratio ou SMR en anglais). Un script directement exécutable est téléchargeable en fin d’article.

Les indices comparatifs standardisés ne sont pas d’usage courant en sociologie (ils le sont surtout en démographie ou en épidémiologie), bien que leur interprétation puisse tout à fait être intégrée à une analyse sociologique [1]. Voici par exemple une graphique qui représente les indices standardisés de mortalité par groupe socio-professionnel en Belgique. Nous allons voir comment calculer ces indices et construire ce graphique :

Indices comparatifs

Avant de nous intéresser à la standardisation, voyons rapidement ce qu’est un indice comparatif. Cet indice permet de comparer la fréquence d’une occurrence entre des groupes. Disons que nous désirons comparer la fréquence des décès entre les différentes groupes socio-professionnels. Une manière de le faire est de mesurer le nombre de décès observés pour chaque groupe et de les comparer aux décès théoriques que l’on observerait si chaque groupe avait la même probabilité de décès (vous excuserez la notation mathématique hasardeuse d’un sociologue !) :

${Indice\:comparatif_{\text{(groupe)}}}=\frac{N_{\text{observé(groupe)}}}{N_{\text{théorique(groupe)}}}$

Ces décès théoriques sont simplement calculés en multipliant les effectifs de chaque groupe par la probabilité de décès moyenne, calculée sur l’ensemble de la population, quel que soit le groupe :

${N_{\text{théorique (groupe)}}} = {effectif\:total_{\text{(groupe)}}}\times\:probabilit\acute{e}\:de\:d\acute{e}c\grave{e}s_{\text{(moyenne générale)}}$

Voilà la manière dont on peut calculer cet indice comparatif avec R. Considérons que nous disposons d’une base de donnée d avec des informations de mortalité (variable d$dead), et que l’on désire calculer des indices comparatifs entre groupes socio-professionnels (variable d$class_j_simp). Au préalable, nous dupliquons les variables que nous désirons utiliser dans des colonnes génériques d$GROUP_SMR (les groupes - ici socio-professionnels - à comparer) et d$occurence (l’occurence de comparaison - dans notre cas le décès). Cela nous permet de définir les variables à utiliser au tout début du script, ce qui nous évitera de changer tout le script par la suite quand nous réaliserons cette même analyse avec d’autres variables.

# Considérons que nous disposons d'une base de donnée d.
 
# Je crée les variables retenues dans des colonnes génériques
# * Une variable de groupe (que l'on veut comparer)
d$GROUP_SMR <- d$class_j_simp
 
# * L'occurence de comparaison (habituellement en démo : le décès)
# Note : la variable doit prendre la forme : 1 = occurence ; 0 = non-occurence
d$occurence <- d$dead

Concernant les indice standardisés, le package dplyr (faisant partie du tidyverse) permet d’en faire le calcul très facilement :

  • D’abord, nous calculons la probabilité de décès moyenne de toute la population dans une colonne quotient à l’aide de la fonction mutate. Toutes les lignes de cette colonne sont similaires et comprennent la même valeur de la probabilité moyenne de décès.
  • Ensuite, nous enclenchons l’agrégation des résultats par groupe socio-professionnel avec la fonction group_by.
  • Enfin, nous calculons avec la fonction summarise les effectifs par groupe (n_tot), le nombre de décès observés par groupe (obs_n_tot), le nombre de décès théoriques si le groupe rencontrait les probabilités de décès moyennes (th_n_tot = probabilité de décès moyenne X effectifs du groupe) et le rapport entre les deux, qui débouche sur les indices comparatifs (SMR).
library(tidyverse)
 
# Calcul de l'indice standardisé de mortalité
SMR <- d %>%
  filter(!is.na(GROUP_SMR) & !is.na(occurence)) %>% # Je supprime les valeurs manquantes
  mutate(quotient = sum(occurence == 1, na.rm = TRUE) / n() ) %>% # Je calcule la probabilité de décès pour l'ensemble de la population
  group_by(GROUP_SMR) %>%
  summarise(age = mean(age_reel),
            n_tot = first(n()), # Ici artifice car si je fais summarise avec n*quotient, il ne réduit pas les lignes => A résoudre pour que ce soit plus propre.
            obs_n_tot = first (sum(occurence == 1, na.rm = TRUE)),
            th_n_tot = first(n_tot*quotient),
            SMR = (obs_n_tot/th_n_tot)*100)

On voit ci-dessous le tableau des résultats produit par le script R précédent. On y constate par exemple que les professions libérales présentent moins de décès que la moyenne (0,82 des décès théoriques) et les ouvriers qualifiés davantage que la moyenne (1,28). Chose étonnante, les grands chefs d’entreprise subissent plus de décès que la moyenne (1,21), alors qu’on aurait pu leur croire des conditions de vie plus favorables que celles du reste des autres catégories :

En réalité, interpréter ces différences en termes d’inégalités de conditions de vie n’est pas correct, car cela fait l’impasse sur un élément essentiel : les différents groupes socio-professionnels ne sont pas du même âge et du même sexe. Les chefs d’entreprise, par exemple, sont le plus âgé des groupes (52,8 ans), et ils sont en grande majorité des hommes (87,8%). Or, la chance de décéder augmente avec l’âge, et est plus importante pour les hommes (qui décèdent plus jeunes que les femmes). Ce qui nous intéresserait, en réalité, c’est l’inégalité face à la mort entre les groupe socio-professionnels à caractéristiques d’âge et de sexe égales ! C’est exactement ce à quoi sert la standardisation.

Indices comparatifs standardisés

Quel est donc le principe de la standardisation ? L’indice comparatif standardisé se calcule, comme précédemment, en divisant des décès observés par des décès théoriques pour chacun des groupes :

${Indice\:comparatif\:standardis\acute{e}_{\text{(groupe)}}}=\frac{N_{\text{observé(groupe)}}}{N_{\text{théorique(groupe)}}}$

La différence réside dans le fait que les décès théoriques ne sont pas calculés de la même manière. Si nous voulons réaliser une standardisation à deux facteurs que sont l’âge et le sexe, il faut calculer séparément la probabilité de décès pour chaque tranche d’âge et ce pour chaque sexe, et appliquer ces probabilités aux effectifs de chaque tranche d’âge, par sexe, des différents groupes, avant de les sommer. On obtient alors des décès théoriques qui tiennent comptent de la composition variable de chaque groupe socio-professionnel en termes d’âge et de sexe :

${N_{\text{théorique (groupe)}}} = {\sum\:effectifs_{\text{(groupe, age, sexe)}}}\times\:probabilit\acute{e}s\:de\:d\acute{e}c\grave{e}s_{\text{(age, sexe)}}$

Comment réalise-t-on ce calcul avec R ? Nous dupliquons au préalable les deux variables de standardisation dans des variables génériques d$STANDARD_SMR_1 et d$STANDARD_SMR_2. Nous passons ensuite au calcul :

  • D’abord, nous calculons la probabilité de décès, non pas pour toute la population, mais par âge et par sexe (avec un premier group_by). La fonction mutate crée ainsi une colonne de probabilités de décès qui varient selon l’appartenance d’âge et de sexe de la personne à la ligne considérée.
  • Ensuite, nous calculons de manière agrégée, avec la fonction summarise, par groupe, âge et sexe (suite à un nouveau group_by), les effectifs (n_tot), le nombre de décès observés (obs_n_tot) et le nombre de décès théoriques (th_n_tot).
  • Enfin, nous sommons les résultats par groupe socio-professionnel (avec un dernier group_by et un dernier summarise). Nous pouvons alors calculer les indices comparatifs standardisés en divisant les décès observé par les décès théoriques (SMR) :
# * Une variable de standardisation 1
d$STANDARD_SMR_1 <- d$age_reel
 
# * Une variable de standardisation 2
d$STANDARD_SMR_2 <- d$Sexe
 
SMR_standard <- d %>%
  filter(!is.na(GROUP_SMR) & !is.na(occurence) & !is.na(STANDARD_SMR_1) & !is.na(STANDARD_SMR_2)) %>% 
  group_by(STANDARD_SMR_1, STANDARD_SMR_2) %>%
  mutate(quotient = sum(occurence == 1, na.rm = TRUE) / n() ) %>% 
  ungroup() %>% 
  group_by(STANDARD_SMR_1, STANDARD_SMR_2, GROUP_SMR) %>%
  summarise(n = first(n()),
            obs_n = first (sum(occurence == 1, na.rm = TRUE)),
            th_n = first(n*quotient)) %>% 
  ungroup() %>% 
  group_by(GROUP_SMR) %>%
  summarise(n_tot = sum(n),
            obs_n_tot = sum(obs_n),
            th_n_tot = sum(th_n),
            SMR = obs_n_tot/th_n_tot)

On peut voir le résultat dans le tableau ci-dessous. On remarque ainsi des résultats très différents : les chefs d’entreprise présentent désormais l’indice standardisé de mortalité le plus bas ! Les agents d’entretien, à l’inverse, montrent maintenant un indice standardisé de mortalité largement supérieur à la moyenne.

Mise en forme avec ggplot2

Les résultats sont déjà marquants et explicites, mais ils le seraient encore davantage dans une mise en forme graphique qui les mette en valeur. Voyons comment réaliser le "graphique divergent" (divergent bar graph) montré en début d’article avec ggplot2. Le principe du graphique est d’afficher des barres par groupes socio-professionnels qui s’écartent d’autant plus de la moyenne (standardisée) que l’indice est élevé (positivement vers la droite, négativement vers la gauche).

En premier lieu, il nous faut préparer les données du tableau SMR que nous avons créé :

  • Nous calculons les intervalles de confiance des indices comparatifs standardisés dans des bornes inférieures et supérieures (SMR$SMR_IC_inf et SMR$SMR_IC_sup) [2] ;
  • Nous créons une nouvelle variable (SMR$group) dont la valeur est 0 si l’indice est significativement négatif (au sens statistique, présumant que nous travaillons sur un échantillon) ; 1 si l’indice est non significatif ; 2 si l’indice est significativement positif. Cette variable nous permettra de faire varier la couleur sur le graphique ggplot2 en vert lorsque positif, orange lorsque non significatif, rouge lorsque négatif ;
  • Nous recodons l’indice en un écart à la moyenne en pourcentage, celui-ci étant plus facilement compréhensible (SMR$SMR2digits et SMR$SMR2digits_char) ;
  • J’utilise une astuce pour afficher les noms des groupes socio-professionnels à droite et à gauche de l’axe central qui représentera la moyenne : nous dupliquons les noms des groupes socio-professionnels dans la colonne SMR$prof0 lorsque l’indice est inférieur à la moyenne, et dans la colonne SMR$prof1 lorsque l’indice est supérieur à la moyenne. Il ne s’agit sans doute pas de la manière de procéder la plus propre, mais elle fonctionne très bien :)
# Je calcule les intervalles de confiance
conf.level <- 0.95 # Je définis le niveau de confiance
z.score <- qnorm((1-conf.level)/2,lower.tail=FALSE) # Je transforme ce niveau de confiance en valeur sur une normale centrée réduite
SMR$SMR_IC_inf <- ((sqrt(SMR$obs_n_tot)-(z.score*0.5))^2)/SMR$th_n_tot
SMR$SMR_IC_sup <- ((sqrt(SMR$obs_n_tot)+(z.score*0.5))^2)/SMR$th_n_tot
 
# Je définis la valeur a/ 0 si SMR significativement négatif b/ 1 si SMR non significatif c/ 2 si SMR significativement positif
SMR$group <- 1
SMR$group[SMR$SMR < 1 & SMR$SMR_IC_sup < 1] <- 0
SMR$group[SMR$SMR > 1 & SMR$SMR_IC_inf > 1] <- 2
SMR$group <- as.factor(SMR$group)
 
# Je calcule un écart à la moyenne en pourcentage et le transforme en caractère en ajoutant un signe et un pourcentage
SMR$SMR2digits <- (SMR$SMR - 1)
SMR$SMR2digits_char <- ifelse(SMR$SMR > 1,
                              paste0("+",
                                     round((SMR$SMR - 1)*100, digits = 0),
                                     "%"),
                              paste0(round((SMR$SMR - 1)*100, digits = 0),
                                     "%")
                              )
 
# Je place les noms des catégories socio-professionnelles dans deux colonnes distinctes selon que le SMR soit inférieur ou supérieur à la moyenne
SMR$prof0 <- ""
SMR$prof0[SMR$SMR < 1] <- as.character(SMR$GROUP_SMR[SMR$SMR < 1])
SMR$prof1 <- ""
SMR$prof1[SMR$SMR > 1] <- as.character(SMR$GROUP_SMR[SMR$SMR > 1])
 
# J'ordonne les données selon les valeurs croissantes du SMR
SMR <- SMR[order(SMR$SMR),]

On peut ensuite construire un graphique ggplot2. J’utilise pour ce faire les fonction geom_bar pour dessiner les barres, geom_hline pour dessiner la moyenne (standardisée) et geom_text pour afficher les noms des groupes socio-professionnels. J’emploie les package hrbrthemes (pour le thème ggplot2 theme_ipsum) et showtext (pour les polices) afin d’enjoliver le graphique.

library(hrbrthemes) # Des thèmes pour ggplot2
library(showtext) # Un package pour importer des polices de caractère
font_add(family = "Arial Narrow",
         regular = "jo_scripts_arial narrow plain.ttf",
         bold = "jo_scripts_arial narrow bold.ttf",
         italic = "jo_scripts_arial narrow italic.ttf",
         bolditalic = "jo_scripts_arial narrow bold italic.ttf")
showtext_auto()
 
# Titres et labels
label_main_SMR_manual <- "Indices comparatifs par catégories socio-professionnelles"
label_sub_SMR_manual <- "Standardisés par âge & sexe"
Label_occurence_SMR <- "Mortalité"
 
# Je construis le graphique
SMR_ggplot <- ggplot(SMR) + 
  geom_bar(aes(x = reorder(GROUP_SMR, -SMR2digits), y = SMR2digits, fill=group), 
           stat = "identity", 
           width = 0.8) + 
  geom_hline(yintercept = 0, # Une ligne verticale noire représentant la moyenne (standardisée)
             linetype="solid", 
             size=1) +
  geom_text(aes(x = reorder(GROUP_SMR, -SMR2digits), y = SMR2digits, label=SMR2digits_char),
            color="white",
            size=4, 
            position = position_stack(vjust = 0.5), 
            vjust = 0.35) + 
  geom_text(aes(x = reorder(GROUP_SMR, -SMR2digits), y = 0.01, label=prof0),
            hjust = "left",
            vjust = 0.35,
            size=4) +
  geom_text(aes(x = reorder(GROUP_SMR, -SMR2digits), y = -0.01, label=prof1),
            hjust = "right",
            vjust = 0.35,
            size=4)  +
  coord_flip() + 
  scale_fill_manual(  # Changement de couleur des barres en fonction de la significativité
    values=c("0" = "#35ac97", "1" = "#e2c688", "2" = "#e94e66"),
    name=Label_occurence_SMR,
    labels=c("0" = "Sous-représentation", "1" = "Ecart non significatif", "2" = "Sur-représentation")
  ) + 
  theme_ipsum() + # thème issu du package "hrbrthemes"
  theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    panel.grid.major.x = element_line(color = "grey", linetype = "dotted"),
    panel.grid.minor.x = element_line(color = "grey", linetype = "dotted"),
    panel.grid.major.y = element_blank(),
    panel.grid.minor.y = element_blank()
  ) +
  ggtitle (label_main_SMR_manual) +
  labs (subtitle = print(label_sub_SMR_manual)) +
  scale_y_continuous(labels = scales::label_percent(accuracy = 1L)) +
  ylab("") +
  xlab("")
 
SMR_ggplot

Voici le résultat final. Joli, non ?