--- title: "Autour du Paradoxe de Simpson" author: "Brahima DIARRA" date: "22 juin 2020" output: pdf_document: default html_document: default --- ```{r setup, message=FALSE, warning=FALSE, include=FALSE} knitr::opts_chunk$set(echo = TRUE, message=FALSE, warning=FALSE) ``` l'objectif de ce document est de reproduire les résultats d'un sondage sur un sixième des électeurs de la ville de Wickham, une ville au nord-est de l'Angleterre dans les années 1972-1974. Cette étude se voulait d'éclairer des travaux sur les maladies thyroïdiennes et cardiaques (Tunbridge et al. 1977). Une suite de cette étude a été menée vingt ans plus tard (Vanderpump et al. 1995). Certains des résultats avaient trait au tabagisme et cherchaient à savoir si les individus étaient toujours en vie lors de la seconde étude. Par simplicité, nous nous restreindrons aux femmes et parmi celles-ci aux 1314 qui ont été catégorisées comme "fumant actuellement" ou "n'ayant jamais fumé". Il y avait relativement peu de femmes dans le sondage initial ayant fumé et ayant arrêté depuis (162) et très peu pour lesquelles l'information n'était pas disponible (18). La survie à 20 ans a été déterminée pour l'ensemble des femmes du premier sondage. Les données utilisées sont téléchargeables en format CSV à cette [adresse](https://gitlab.inria.fr/learninglab/mooc-rr/mooc-rr-ressources/blob/master/module3/Practical_session/Subject6_smoking.csv). chaque ligne renseigne si la personne fume ou non, si elle est vivante ou décédée au moment de la seconde étude, et son âge lors du premier sondage. # Tâches préliminaires Cette section va consister à faire les travaux préliminaires avant d'aborder les tâches mentionnées dans l'étude de cas. Ces tâches sont : 1. Représenter dans un tableau le nombre total de femmes vivantes et décédées sur la période en fonction de leur habitude de tabagisme ; 2. Reprendre la question 1 (effectifs et taux de mortalité) en rajoutant une nouvelle catégorie liée à la classe d'âge ; 3. Envisager d'essayer de réaliser une régression logistique afin d'éviter un biais induit par des regroupements en tranches d'âges arbitraires et non régulières. ## Téléchargement des données Nous allons télécharger les données, en utilisant l'adresse du lien indiqué ci-dessus, pour les stocker sur notre machine pour ne pas avoir à les charger à chaque fois qu'on veut relancer l'analyse. ```{r} #Url de téléchargement data_url <- "https://gitlab.inria.fr/learninglab/mooc-rr/mooc-rr-ressources/-/raw/master/module3/Practical_session/Subject6_smoking.csv?inline=false" #Test d'existence du fichiers de données data_file <- "subject6_smoking.csv" if (!file.exists(data_file)){ download.file(data_url, data_file, method = "auto") } ``` ## Lecture des données Après avoir télécharger les données et les stocker sur le disque local, nous pouvons les charger dans R et procéder à l'inspection. ```{r} data <- read.csv(data_file) head(data) tail(data) ``` On peut aussi vérifier la présence d'observations manquantes ```{r} lignes_na <- apply(data, 1, function(x)any(is.na(x))) data[lignes_na,] ``` La sortie nous indique qu'il n'y a pas d'observations manquantes dans les données. On peut procéder aux analyses ## Installation et chargement des packages utilisées Nous utilisons une série de package pour pouvoir faire cette analyse, nous allons d'abord tester si un package parmi la liste n'est pas installé, on l'installe puis on charge l'ensemble des packages nécessaires. ```{r packg_loading, echo = F} #Liste des packages utilisés dans l'analyse pkgs <- c( "dplyr", ###Manipuler les données "ggplot2", ###Faire ds graphiques "janitor", ###Tableaux croisés avec totaux lignes et colonnes "knitr", ###Formater les tableaux "kableExtra" ###Mise en forme des tableaux ) #Test et intallation des packages lapply(pkgs[!(pkgs %in% installed.packages())], install.packages) ## Charger tous les packages lapply(pkgs, library, character.only = TRUE) ``` # Mission dévolues par le MOOC ## 1. Représentation dans un tableau du nombre total de femmes vivantes et décédées sur la période en fonction de leur habitude de tabagisme. ### Nombre total Pour avoir ce tableau, nous utilisons la fonction *tabyl* du package [janitor](https://cran.r-project.org/web/packages/janitor/vignettes/tabyls.html) qui lui-même se sert de [dplyr](https://dplyr.tidyverse.org/). ```{r} data %>% tabyl(Status, Smoker) %>% adorn_totals(c("row", "col")) %>% adorn_title(placement = "combined") %>% kable(align = c("l",rep("c",3)),caption = "Nombre total de femmes vivantes et décédées sur la période en fonction de leur habitude de tabagisme") %>% kable_styling(bootstrap_options = c("striped", "hover","responsive"), fixed_thead = T, font_size = 13, full_width = T, position = "center") %>% column_spec(1, bold = T) %>% #row_spec(seq(1,nlevels(data$Status)+2, by=1), background = "#D5E4EB") footnote(general = "Sondage auprès d'un sixième des électeurs afin d'éclairer des travaux sur les maladies thyroïdiennes et cardiaques en 1972-1974 à Whickham.", general_title = "Source: ", footnote_as_chunk = T) ``` Il ressort du tableau précédent que sur les 1314 femmes, 369 sont décédées sur la période et 945 étaient toujours en vie. Parmi les 369 femmes décédées, 230 ne fumaient pas alors que 130 fumaient. ### Calcul du taux de mortalité suivant le statut de tabagisme Le taux de mortalité est le rapport entre le nombre de femmes décédées dans un groupe et le nombre total de femmes dans ce groupe. ```{r} data %>% tabyl(Status, Smoker) %>% adorn_totals(c("row", "col")) %>% adorn_title(placement = "combined") %>% adorn_percentages("col") %>% adorn_pct_formatting(digits = 1, affix_sign = F) %>% kable(align = c("l",rep("c",3)),caption = "Taux de mortalité en fonction de l'habitude de tabagisme") %>% kable_styling(bootstrap_options = c("striped", "hover","responsive"), fixed_thead = T, font_size = 13, full_width = T, position = "center") %>% column_spec(1, bold = T) %>% #row_spec(seq(1,nlevels(data$Status)+2, by=1), background = "#D5E4EB") footnote(general = "Sondage auprès d'un sixième des électeurs afin d'éclairer des travaux sur les maladies thyroïdiennes et cardiaques en 1972-1974 à Whickham.", general_title = "Source: ", footnote_as_chunk = T) ``` Ce tableau montre que le taux de mortalité est de 31,4% pour les femmes qui ne fumaient pas et 23,9% pour celles qui fumaient. Le graphique suivant permet de mieux cerner les différences entre les deux groupes. ```{r} data_tx <- data %>% group_by(Smoker) %>% count(Status) %>% mutate(ratio=scales::percent(n/sum(n), decimal.mark = ",", accuracy = 0.1)) ggplot(data,aes(x=Smoker,fill=Status))+ geom_bar(position="fill")+ geom_text(data=data_tx, aes(y=n,label=ratio), position=position_fill(vjust=0.5))+ scale_fill_manual(values = c("green3", "red1"))+ labs(title = "Un taux de mortalité plus élevé dans le groupe des \n non-fumeuses que dans celui des fumeuses", x = "Habitude de tabagisme", fill = "Statut de décès", caption = "Source: Sondage auprès d'un sixième des électeurs afin d'éclairer \n des travaux sur les maladies thyroïdiennes et cardiaques \n en 1972-1974 à Whickham.")+ theme_bw() ``` Le test de Khi2 sur ces deux variables permet de conclure que ces différences sont significatives comme le montre la sortie suivante: ```{r} chisq.test(data$Smoker, data$Status) ``` De l'analyse des tableaux qui précède, il est constaté que le taux de mortalité est plus importante dans le groupe des femmes qui ne fumaient pas (31,4%) que dans celui des femmes qui fumaient (23,9%). Ce résultat est étonnant dans la mesure où le tabac a, de manière générale, une incidence négative sur la santé des fumeurs, en conséquence, le taux de mortalité devrait être plus important dans le groupe des femmes qui fumaient. ## 2. Prise en compte de la classe d'âge dans l'analyse ### Création des classes d'âges On considérera les classes suivantes dans l'analyse: 18-34 ans, 34-54 ans, 55-64 ans, plus de 65 ans. Pour ce faire, on peut utiliser la fonction *cut* de base ou la fonction *icut* du [package *questionner*](https://cran.r-project.org/web/packages/questionr/index.html) ```{r discretisation, echo=FALSE} summary(data$Age) data$agecl <- cut(data$Age, breaks = c(18, 34, 54, 64, 90), labels = c("18-34 ans", "35-54 ans", "55-64 ans","plus de 65 ans"), include.lowest = TRUE) tabyl(data$agecl) %>% adorn_totals() %>% adorn_pct_formatting(digits = 1, affix_sign = T) ``` La répartition en classe d'âge montre que 30,4% des femmes de l'étude sont dans la tranche d'âge 18-34 ans, le tiers d'entre elles ont entre 35 et 54 ans et 36,4% ont 55 ans ou plus. ### Croisement de la mortalité avec les classes d'âges Le tableau suivant permet de croiser la mortalité en fonction de la tranche d'âge. ```{r} data %>% tabyl(Status, agecl) %>% adorn_totals(c("row", "col")) %>% adorn_title(placement = "combined") %>% kable(align = c("l",rep("c",3)),caption = "Nombre de décès en fonction de la tranche d'âge") %>% kable_styling(bootstrap_options = c("striped", "hover","responsive"), fixed_thead = T, font_size = 13, full_width = T, position = "center") %>% column_spec(1, bold = T) %>% #row_spec(seq(1,nlevels(data$Status)+2, by=1), background = "#D5E4EB") footnote(general = "Sondage auprès d'un sixième des électeurs afin d'éclairer des travaux sur les maladies thyroïdiennes et cardiaques en 1972-1974 à Whickham.", general_title = "Source: ", footnote_as_chunk = T) ``` Il ressort de ce tableau que l'effectif du nombre de mort a tendance en augmenté en fonction de l'âge des enquêtés. ```{r} data %>% tabyl(Status, agecl) %>% adorn_totals(c("row", "col")) %>% adorn_title(placement = "combined") %>% adorn_percentages("col") %>% adorn_pct_formatting(digits = 1, affix_sign = F) %>% kable(align = c("l",rep("c",3)),caption = "Taux de mortalité en fonction de la tranche d'âge") %>% kable_styling(bootstrap_options = c("striped", "hover","responsive"), fixed_thead = T, font_size = 13, full_width = T, position = "center") %>% column_spec(1, bold = T) %>% #row_spec(seq(1,nlevels(data$Status)+2, by=1), background = "#D5E4EB") footnote(general = "Sondage auprès d'un sixième des électeurs afin d'éclairer des travaux sur les maladies thyroïdiennes et cardiaques en 1972-1974 à Whickham.", general_title = "Source: ", footnote_as_chunk = T) ``` Les fréquences relatives confirment le constat précédent. En effet, le taux de mortalité global est 28,1%; ce taux est de 2,8% pour les femmes ayant moins de 35 ans, 13,8% pour celles ayant entre 35 et 54 ans et plus de 85% pour les femmes qui ont 65 ans ou plus. Voyons à présent l'association eventuelle entre l'âge et le tabagisme. ```{r} data %>% tabyl(Smoker, agecl) %>% adorn_totals(c("row", "col")) %>% adorn_title(placement = "combined") %>% adorn_percentages("col") %>% adorn_pct_formatting(digits = 1, affix_sign = F) %>% kable(align = c("l",rep("c",3)),caption = "Habitude de tabagisme en fonction de la tranche d'âge") %>% kable_styling(bootstrap_options = c("striped", "hover","responsive"), fixed_thead = T, font_size = 13, full_width = T, position = "center") %>% column_spec(1, bold = T) %>% #row_spec(seq(1,nlevels(data$Status)+2, by=1), background = "#D5E4EB") footnote(general = "Sondage auprès d'un sixième des électeurs afin d'éclairer des travaux sur les maladies thyroïdiennes et cardiaques en 1972-1974 à Whickham.", general_title = "Source: ", footnote_as_chunk = T) ``` Aucun lien ne semble se dégager entre ces deux variables au vu du tableau, par contre le test du khi2 conclue au rejet d'absence de liaisons entre elles. ```{r} chisq.test(data$Smoker, data$agecl) ``` Le graphique suivant donne la réprésentation graphique de la mortalité en fonction de l'âge et du tabagisme. ```{r} data_tx <- data %>% group_by(Smoker, agecl) %>% count(Status) %>% group_by(Smoker, agecl) %>% mutate(n_smok = sum(n)) %>% mutate(ratio=scales::percent(n/n_smok, decimal.mark = ",", accuracy = 0.1)) ggplot(data,aes(x=Smoker,fill=Status))+ geom_bar(position="fill")+ geom_text(data=data_tx, aes(y=n,label=ratio), position=position_fill(vjust=0.5))+ scale_fill_manual(values = c("green3", "red1"))+ labs(title = "Un taux de mortalité qui croit en fonction de l'âge", x = "Habitude de tabagisme", y = "Taux de décès", fill = "Statut de décès", caption = "Source: Sondage auprès d'un sixième des électeurs afin d'éclairer \n des travaux sur les maladies thyroïdiennes et cardiaques \n en 1972-1974 à Whickham.")+ facet_grid(~agecl)+ theme_bw() ``` Ce graphique montre que le taux de mortalité est affecté par l'âge pour les femmes se sitaunt dans les tranches d'âges 35-54 et 55-64 ans. En effet, le taux de mortalité pour les femmes qui fumaient de la tranche d'âge 35-54 ans est plus important (17,3%) que pour celles qui ne fumaient pas (9,5%). Ce constat est le même pour la tranche d'âge 55-64 ans avec 44,3% de taux de mortalité pour les femmes qui fumaient contre 33,1% pour celles qui ne fumaient pas. Il ressort ainsi que ces deux tranches d'âges contredisent la conclusion générale tiré à la première question selon laquelle, le taux de mortalité était plus important dans le groupe des femmes qui ne fumaient pas que dans celui des femmes qui fumaient. Cela pourrait s'expliquer par l'importance de la mortalité dans la tranche d'âge *plus de 65 ans* dans laquelle plus de 56% de la mortalité est concentrée. Le tabagisme n'ayant pas d'effet sur le taux de mortalité dans cette tranche (le taux de mortalité sont sensiblement égaux pour les fumeuses ainsi que pour les non-fumeuses), l'effet du tabagisme sur la mortalité se trouve attenuer par ce taux important de mortalité de cette tranche d'âge. Ce paradoxe est appelé [*paradoxe de Simpson*](https://fr.wikipedia.org/wiki/Paradoxe_de_Simpson) ## 3. Utilisation de la régression logistique Afin d'éviter un biais induit par des regroupements en tranches d'âges arbitraires et non régulières, nous envisageasons de réaliser une régression logistique. Si on introduit une variable Death valant 1 ou 0 pour indiquer si l'individu est décédé durant la période de 20 ans, on peut étudier le modèle Death ~ Age pour étudier la probabilité de décès en fonction de l'âge selon que l'on considère le groupe des fumeuses ou des non fumeuses. ### Création de la variable Death Pour créer cette variable, on utilise la fonction *if_else* du package de base. ```{r creation_Death} data$Death <- ifelse(data$Status == "Dead", 1, 0) #Vérification du codage with(data, table(Status, Death)) ``` ### Estimation du modèle et représentation graphique ```{r estimation_model} logistic_reg = glm(data=data, Death ~ Age, family=binomial(link='logit')) summary(logistic_reg) ``` L'estimation du modèle montre que l'âge a effectivement un impact sur la probabilité de décès. L'augmentation de l'âge d'une unité influence positivement la probabilité de décès des femmes enquêtés indépendamment de l'habitude de tabagisme. Le graphique suivant est relatif à la probabilité de décès en fonction de l'âge selon que l'on considère le groupe des fumeuses ou des non fumeuses. ```{r graph_model} ggplot(data, aes(x=Age, y=Death)) + geom_point() + stat_smooth(method="glm", method.args=list(family="binomial"), se=TRUE)+ labs(title = "Représentation graphique de la probabilité de décès\n en fonction de l'âge suivant l'habitude de tabagisme", x = "Age du répondant", y = "Probabilité de décès")+ facet_grid(~Smoker)+ theme_bw() ``` Ce graphique montre que la probabilité de décès augmente en fonction de l'âge quel que soit l'habitude de tabagisme du répondant. La region de confiance est plus étroite pour les non-fumeuses que pour les fumeuses.