--- title: "Sujet 6 : Autour du Paradoxe de Simpson" author: "Louis Hognon" date: "19/04/2020" output: html_document : toc : true toc_depth: 3 --- # Préparation des données #### Les données **autour du Paradoxe de Simpson** se trouvent sont acessibles via le site de [gitlab INRIA](https://gitlab.inria.fr/learninglab/mooc-rr/mooc-rr-ressources/blob/master/module3/Practical_session/Subject6_smoking.csv) L'URL de ce lien est : ```{r} data_url= "https://gitlab.inria.fr/learninglab/mooc-rr/mooc-rr-ressources/-/raw/master/module3/Practical_session/Subject6_smoking.csv?inline=false" ``` Pour **eviter une éventuelle disparition des données du serveur gitlab**, nous faisons une **copie locale de ce jeux de données** que nous préservons Il est inutile et même risquée de télécharger les données à chaque exécution, car dans le cas d'une panne nous pourrions remplacer nos données par un fichier défectueux. Pour cette raison, nous téléchargeons les données seulement si la copie locale n'existe pas. ```{r} data_file = "simpsons.csv" if (!file.exists(data_file)) { download.file(data_url, data_file, method="auto") } ``` ## Téléchargement des données ### Nous téléchargons les données via la commande **read.csv**, appliquée à data_url et data_file, dans le cas d'une éventuelle disparition des données du serveur gitlab. ```{r} data_exo6 = read.csv(data_url) data_exo6 = read.csv(data_file) ``` # Description de la base de donnée ### Avant de commencer les analyses sur la base de données, nous allons d'abord réaliser une description de celle-ci. Puis vérifier si certains paramètres : type de données, données manquantes, données abberantes, ## Afficher la structure de la base de donnée Avec la commande **str** nous allons vérifier le type de données dans chaque colonne ainsi que les premières valeurs. ```{r} str(data_exo6) ``` On retrouve bien nos 3 colonnes: - Smoker : si la personne fume ou non - Status : si elle est vivante ou décédée au moment de la seconde étude - Age : son âge lors du premier sondage Pour rappel : *En 1972-1974, à Whickham, une ville du nord-est de l'Angleterre, située à environ 6,5 kilomètres au sud-ouest de Newcastle upon Tyne, un sondage d'un sixième des électeurs a été effectué afin 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.* ## Vérifications : données abbérantes On utilise la commande summary afin de voir s'il y a des données abbérantes dans chacune des colones. On précise la colone à inspecter grâce à la commande $. ```{r} summary(data_exo6$Smoker) summary(data_exo6$Status) summary(data_exo6$Age) ``` Il n'y a pas de données abbérantes, cad, des réponses autres que No ou Yes, Alive ou Dead et concernant l'âge pas d'âge < 0 et > 100 ans. ## Vérificaitons : données manquantes ```{r} na_records = apply(data_exo6, 1, function (x) any(is.na(x))) data_exo6[na_records,] ``` Il n'y a aucune données manquantes # Exercices ## Partie 1 ### **Consgines : ** Représentez 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. On utilise la fonction table avec les variables demandées : ici Smokers et status ```{r} tableau_mortalite= table(data_exo6$Smoker,data_exo6$Status) tableau_mortalite ``` ### **Consignes :**Calculez dans chaque groupe (fumeuses / non fumeuses) le taux de mortalité (le rapport entre le nombre de femmes décédées dans un groupe et le nombre total de femmes dans ce groupe). On utilise la fonciton prop.table pour avoir les taux de mortalité dans chaque groupe. La fonction round permet d'avoir un arrondi du résultat ```{r} round(((prop.table((tableau_mortalite)))*100),2) ``` On obtient alors un pourcentage de mortalité de : - 17.50 % : Groupe Non Fumeuse - 10.58 % : groupe Fumeuse Afin de ne pas utiliser manuellement les données écrites ci-dessus, nous allons les conserver dans des vecteurs ```{r} tableau = round(((prop.table(table(data_exo6$Smoker,data_exo6$Status)))*100),2) #creation d'un data frame à partir du tableau export_data_fram = as.data.frame(tableau) # selection des colonnes et variables qui nous interesse et creation d'un nouveau data fram new_data_frame = as.data.frame(export_data_fram[3:4,3]) # attribution des valeurs de mortalite pour les groupes à des vecteurs mortalite_non_fum = new_data_frame[1,] mortalite_fumeuse = new_data_frame[2,] ``` ###**Consignes : ** Vous pourrez proposer une représentation graphique de ces données et calculer des intervalles de confiance si vous le souhaitez. Pour la représentation graphique nous allons procéder en 2 étapes : - créer un data.frame, que l'on nomme figure de prop.table((tableau_mortalité)) avec la fonction as.data.frame - utiliser le package ggplot2 pour faire une représentation graphique de ce data frame ```{r} figure = as.data.frame((prop.table((tableau_mortalite)))*100) library(ggplot2) ggplot(figure, aes(Var1,Freq,fill = Var2)) + geom_bar(stat = "identity",position= "dodge") + labs(title="Statut vital en fonction de la consommation tabagique", x ="Smoker", y = "Fréquence",fill = "Statut vital") +theme_classic()+ theme(plot.title = element_text(color="blue", size=14, hjust= 0.7, face="bold.italic")) + scale_y_continuous( limits=c(0, 45)) ``` Pour le calcul des intervals de confiance nous allons avoir besoin d'utiliser le package binom. ```{r, echo=TRUE} library(binom) ``` Puis on utilise la fonction binom.confint en precisant la proportion de mortalité dans le groupe fumeuse et ensuite non fumeuse en fonction de la population totale ```{r} # Tout d'abord je vais attribuer à la variable N le nombre total de femme dans la population N = length(data_exo6$Age) # On utilise la fonction binom.confint, et ici on choisit la methode exact IC_mort_fume = binom.confint(mortalite_fumeuse, N, methods = "exact") IC_mort_non_fum = binom.confint(mortalite_non_fum, N, methods = "exact") IC_mort_fume IC_mort_non_fum ``` ### **Consignes : **En quoi ce résultat est-il surprenant ? Ces données parraissent surprennantes puisqu'elles sont contradictoires avec la [littérature](https://www.em-consulte.com/rmr/article/157296). Les femmes qui fumaîent semblent avoir un taux de mortalité amoindrie en comparaison avec celles qui ne fumaient pas. ## Partie 2 ### **Consignes : ** Reprenez la question 1 (effectifs et taux de mortalité) en rajoutant une nouvelle catégorie liée à la classe d'âge. On considérera par exemple les classes suivantes : 18-34 ans, 34-54 ans, 55-64 ans, plus de 65 ans. Nous allons créer une nouvelle variable que l'on nomme age_classe correspondant aux différentes catégories d'âges. ```{r} # D'abord recuperer l'âge maximal sur la variable Age max = max(data_exo6$Age) data_exo6$Age_classe = cut(data_exo6$Age, c(18,34.01,54,64,max),right = TRUE,include.lowest = TRUE,labels = c("18-34 ans", "34-54 ans", "55-64 ans", "> 65 ans")) table(data_exo6$Age_classe) #On verifie qu'on obtient le même nombre d'individu : ici 1314 codée avec la variable N length(data_exo6$Age_classe)== N ``` Puis nous allons créer un nouveau tableau en rentrant la variable Age_classe ```{r} new_tableau_mortalite= table(data_exo6$Smoker,data_exo6$Status,data_exo6$Age_classe) new_tableau_mortalite ``` Puis nous allons regarder le pourcentage de décès dans chacun des groupes ```{r} round(((prop.table((new_tableau_mortalite)))*100),2) tab_mort_classe_age = as.data.frame((prop.table(new_tableau_mortalite))*100) tab_mort_classe_age ``` Nous souhaitons disposer que du statut "Dead" pour cela nous allons créer un data frame en utilisant la fonction data.frame et utiliser la fonction filter, en utilisant avant le package dplyr pour selectionner dans la variable dead. ```{r} names(tab_mort_classe_age) names(tab_mort_classe_age)[2]= "Status" library(dplyr) filter(tab_mort_classe_age, Status == "Dead") #Création du nouveau data_frame new_tab_mort_classe_age = as.data.frame(filter(tab_mort_classe_age,Status == "Dead")) ``` Réalisation d'une figure de la mortalité selon l'âge et le statut tabagique Pour cela nous utilison le package ggplot2 ```{r} library(ggplot2) ggplot(new_tab_mort_classe_age, aes(Var3, Freq, fill = Var1)) + geom_bar(stat = "identity",position= "dodge") + labs(title="Pourcentade de décès en fonction \n de la consommation tabagique selon la classe d'âge", x ="Smoker", y = "Fréquence",fill = "Fumeuse")+ theme_classic() + theme(plot.title = element_text(color="blue", size=14, hjust= 0.5, face="bold.italic")) ``` ###En quoi ce résultat est-il surprenant ? Arrivez-vous à expliquer ce paradoxe ? De même, vous pourrez proposer une représentation graphique de ces données pour étayer vos explications. On remarque que si l'on raissone par classe d'âge on obtient des différences de mortalité entre les groupes fumeuses et non fumeuses. Cependant, on observe que les femmes non fumeuses décèdent **3 fois plus** dans la tranche d'âge 64 et plus. Cela pourrait s'expliquer par une moyenne d'âge plus élevée Nous pouvons le vérifier en : - calculant la moyenne d'âge dans le groupe fumeuse et non fumeuse avec la fonction tapply - realiser un boxplot avec la fonction plot - réaliser un test t de student en supposant que la normalité et la variance sont respectée ```{r} tapply(data_exo6$Age,data_exo6$Smoker, mean) plot(data_exo6$Smoker,data_exo6$Age) # Verfication normalité distribution et normalité de la variance hist(data_exo6$Age) by(data_exo6$Age,data_exo6$Smoker, sd, na.rm = TRUE) t.test(data_exo6$Age~data_exo6$Smoker, var.equal=TRUE) ``` On peut compléter ces résultats en calculant la moyenne d'âge dans chacun des groupes dans la tranche d'âge 65 ans et +. Pour cela nous allons : - utiliser le package dplyr pour utiliser la fonction filter - creer un nouveau data frame contenant seulement les individu > 65 ans - calculer la moyenne dans chacun des groupes ```{r} library(dplyr) verification_age = as.data.frame(filter(data_exo6, Age_classe== "> 65 ans")) by(verification_age$Age,verification_age$Smoker,mean) ``` Il y a une différence significative d'âge entre les deux groupes. Le groupe non fumeuse est significativement plus âgée. De plus il y a un pourcentage de femme non fumeuse importante et par conséquent cela renforce les chiffres de mortalité à la fin. Plus de femme âgée et non fumeuse en comparaison avec peu de femme fumeuse et moins jeunes. Toutefois on n'observe pas de différence dans la tranche d'âge 65 ans et plus. ## Partie 3 ##### ### **Conignes :** Afin d'éviter un biais induit par des regroupements en tranches d'âges arbitraires et non régulières, il est envisageable d'essayer 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. Ces régressions vous permettent-elles de conclure sur la nocivité du tabagisme ? Vous pourrez proposer une représentation graphique de ces régressions (en n'omettant pas les régions de confiance). On commence par créer une nouvelle variable Death avec la fonction recode Alive = 0 Dead = 1 et qui devra être numérique ```{r} library(questionr) library(dplyr) data_exo6$Death <- recode(data_exo6$Status, "Alive" = "0", "Dead" = "1") data_exo6$Death <- as.numeric(as.character(data_exo6$Death)) class(data_exo6$Death) ``` Realisation de la regression logistique avec la fonction glm ```{r} #on va utiliser utiliser une regression logistique multiple avec fonction glm # on precise binomial car on se trouve dans une regression logistique mod4= glm(Death~Age*Smoker, data= data_exo6, family="binomial") summary(mod4) # Age * Smokers : permet de voir l'interraction entre les variables age et Smokers ``` Dans le cadre d’un modèle logistique, habituelement on ne présente pas les coefficients du modèle mais leur valeur exponentielle, cette dernière correspondant en effet à des odds ratio Pour cela on utilise le package **questionr** avec la fonction odds.ratio() ```{r} library(questionr) odds.ratio(mod4) ``` Pour réaliser une représentation graphique des regression logistiques, nous allons d'abord utiliser le package **forestmodel** et la fonction forest_model. Cela permet d'afficher une représentation graphique visuelle et tabulaire ```{r} # la fonction forest_model va permettre de representer un tableau d'ODDS ratio du modèle de regression logistique et de présenter les p-values correspondantes # On doit d'abord regarder l'interraction, si elle n'est pas significative on n'accorde pas d'importante aux facteurs principaux. library(forestmodel) forest_model(mod4) ``` Pour réaliser une représentation graphique de la régression logistique en fonction du statut tabagique on peut utiliser la fonction ggeffects ```{r} library(effects) library(ggeffects) #Pour afficher un tableau ggeffect(mod4, "Smokers") #Pour afficher un graphique plot(ggeffect(mod4, "Smoker")) + labs(title = "Prediction de la mortalité ") + theme(plot.title = element_text(color="blue", size=20, hjust= 0.5, face="bold.italic")) ``` A partir de ce dernier graphique on peut induire que les femmes fumeuses ont plus de risque de décès que les femmes non fumeuses. # En conclusion :####### Ces résultats ne permettent pas de conclure sur la nocivité du tabac. En effet, dans un premier temps il apparaissait que les femmes fumeuses avaient un taux de mortalité amoindrie en comparaison avec les non fumeuses.Cela pouvait amener à première vue à dire que les femmes fumeuses ont un taux de mortalité amoindrie Néanmoins, cet effet sur l'ensemble de la population semble s'inverser lorsqu'on réalise des catégories sur les tranches d'âges. De plus en réalisant ces tranches d'âges, on peut vérifier la moyenne d'âge dans les groupes fumeuses et non fumeuses. On observe alors une différence significative. Les femmes non fumeuses sont plus âgées en moyennes de façon significative comparée aux femmes fumeuses. Enfin en réalisant une régression logistique, aucun résultat significatif apparaît entre la mortalité et l'interaction entre l'âge et le statut tabagique. Ce résultats sont liés à des éléments qui ne sont pas pris en compte (comme la présence de variables non indépendantes ou de différences d'effectifs entre les groupes, et la manière d'observer et d'analyser les variables) Le youtubeur et bloggeur [Science Etonnante](https://sciencetonnante.wordpress.com/2013/04/29/le-paradoxe-de-simpson/) l'explique parfaitement dans son article.