From 9d218ffe669d8f5c887c72b5569d0fb982b2ad0f Mon Sep 17 00:00:00 2001 From: Thibault Verrier Date: Thu, 2 Jul 2020 10:23:30 +0200 Subject: [PATCH] =?UTF-8?q?D=C3=A9p=C3=B4t=20du=20document=20computationne?= =?UTF-8?q?l?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- module3/exo3/exercice_fr.Rmd | 299 +++++++++++++++++++++++++++++++++-- 1 file changed, 284 insertions(+), 15 deletions(-) diff --git a/module3/exo3/exercice_fr.Rmd b/module3/exo3/exercice_fr.Rmd index 7eece5e..408d9bf 100644 --- a/module3/exo3/exercice_fr.Rmd +++ b/module3/exo3/exercice_fr.Rmd @@ -1,8 +1,8 @@ --- -title: "Votre titre" -author: "Votre nom" -date: "La date du jour" -output: html_document +title: "Sujet 6 : Autour du Paradoxe de Simpson" +author: "Thibault Verrier" +date: "11 Juin 2020" +output: pdf_document --- @@ -10,24 +10,293 @@ output: html_document knitr::opts_chunk$set(echo = TRUE) ``` -## Quelques explications +## Préparation des données +Les données décrivent le tabagisme de 1314 femmes de la ville de Whickham, en Angleterre sur la période 1972-1974, et si elles sont vivantes en 1995. Les données du sujet 6 sont disponibles sur le git [**mooc-rr-ressources**](https://gitlab.inria.fr/learninglab/mooc-rr/mooc-rr-ressources/-/raw/master/module3/Practical_session/Subject6_smoking.csv?inline=false) à l'adresse : -Ceci est un document R markdown que vous pouvez aisément exporter au format HTML, PDF, et MS Word. Pour plus de détails sur R Markdown consultez . -Lorsque vous cliquerez sur le bouton **Knit** ce document sera compilé afin de ré-exécuter le code R et d'inclure les résultats dans un document final. Comme nous vous l'avons montré dans la vidéo, on inclue du code R de la façon suivante: +```{r} +data_url <- "https://gitlab.inria.fr/learninglab/mooc-rr/mooc-rr-ressources/-/ +raw/master/module3/Practical_session/Subject6_smoking.csv?inline=false" +``` + +Ces données comportent les colonnes suivantes : + +| Nom de colonne | Modalitées | +|----------------+-----------------------------------------------------------------------------------------------------------------------------------| +| `Smoker` | Yes/No - Oui/Non | +| `Status` | Alive/Dead - Vivante/Morte | +| `Age` | Age de l'individu au moment de l'étude de 1972-1974 | + | + + +### Téléchargement et création du fichier CSV local +Le fichier "SympsonParadox.csv" est téléchargé à l'adresse fournie. On créée le fichier local contenant ces données si il n'est pas déjà présent. + +```{r} +data_file <- "SimpsonParadox.csv" +if(!file.exists(data_file)) { + download.file(url = data_url, destfile = data_file, method = "auto") +} +data <- read.csv(data_file, header = TRUE) +``` +On vérifie que le jeu de données a été correctement importé : + +```{r} +head(data) +tail(data) +``` +Les données ont été importées correctement, on a le même nombre d'individus que donné dans l'énnoncé. + +## Inspection des données + +La bibliothèque dplyr est utilisée pour manipuler plus facilement les *data.frame* +```{r, echo = FALSE} +library(dplyr) +``` + +On vérifie qu'il n'y a pas de données manquantes : +```{r} +na_records <- apply(data, 1, function (x) any(is.na(x))) +data[na_records,] +``` + +On vérifie que l'on retrouve bien le nombre total d'individus en additionnant les différentes modalités, afin de vérifier que les données sont bien homogènes : +```{r} +sum(data$Smoker=="No") + sum(data$Smoker=="Yes") +sum(data$Status=="Alive") + sum(data$Status=="Dead") +``` -```{r cars} -summary(cars) +## Comparaison des taux de mortalité +On réalise ensuite un tableau de contingence afin de mieux comparer le nombre d'individues en vie ou non, selon si elles fumaient lors de la première étude. + +```{r} +tableContingeance <- table(data$Smoker, data$Status) +tableContingeance +``` +On calcule le taux de mortalité pour chaque groupe (fumeuses/non-fumeuses) : $$Mortalité_{f} = n_{Morts_{f}} /n_{Total_{f}}$$ +Avec $f=fumeuse, non-fumeuse$, $n_{Morts_{f}}$ le nombre d'individus possédant la modalité ```Dead``` et $n_{Total_{f}}$ le nombre total d'individus partageant soit la modalité ```Yes```(fumeuses), soit ```No``` (non-fumeuses) pour leur tabagisme. On calcule ainsi les taux de mortalités des individus (```dcdFumeuses``` et ```dcdNonFumeuses```). + +On calcule tout d'abord celui des femmes fumeuse : +```{r} +totFumeuses <- sum(tableContingeance[2,]) +dcdFumeuses <- tableContingeance[2,2] +mortaliteFumeuses <- dcdFumeuses/totFumeuses +mortaliteFumeuses +``` +On obtient ainsi un taux de 0.24 pour les femmes fumeuse. Puis on calcule celui des femmes non-fumeuses : +```{r} +totNonFumeuses <- sum(tableContingeance[1,]) +dcdNonFumeuses <- tableContingeance[1,2] +mortaliteNonFumeuses <- dcdNonFumeuses/totNonFumeuses +mortaliteNonFumeuses +``` +On obtient un taux de 0.32 pour les femmes non-fumeuses. +La table de contingeance est convertie en *dataframe* afin d'être utilisée pour réaliser le graphique. La colonne ```Freq``` qui a été ajoutée représente les effectifs selon les modalitées ```Smoker``` et ```Status```. + +```{r} +tableContingeance <- data.frame(tableContingeance) +tableContingeance <- rename(tableContingeance, Smoker = Var1, Status = Var2) +head(tableContingeance) +``` + +On ajoute ensuite la colonne ```PercentBySmoker``` qui représente le pourcentage d'individus morts ou vivants selon le tabagisme et le statut. Pour cela, on divise ```Freq``` par le nombre total d'individus partageant les mêmes habitudes de tabagisme. + +```{r} +tableContingeance <- mutate(tableContingeance, + PercentBySmoker = ifelse(Smoker == "No", + Freq/sum(Freq[Smoker == "No"]), + Freq/sum(Freq[Smoker == "Yes"]) )) +head(tableContingeance) +``` +Le graphique suivant présente ainsi le taux de mortalité des individus selon leur tabagisme en 1972-1974 à partir des valeurs de la variable ```PercentBySmoker``` + +```{r} +library(ggplot2) +ggplot(tableContingeance, aes(fill=Status, y=PercentBySmoker, x=Smoker)) + + geom_bar(position="stack", stat="identity")+ + ggtitle("Taux de mortalité selon le tabagisme en 1972-1974")+ + theme(plot.title = element_text(hjust = 0.5))+ + xlab("Tabagisme")+ylab("Taux de mortalité")+ + guides(fill=guide_legend(title="Status")) +``` +Etonnemment on trouve un taux de mortalité bien supérieur chez les femmes non-fumeuses que chez les femmes fumeuses. Cette constatation semble contre-intuitive dans la mesure où il est connu que fumer a de nombreux effets nocifs sur la santé. On se serait donc attendu à observer le phénomène inverse, c'est-à-dire une mortalité supérieur chez les femmes fumeuse. + +## Séparation selon les classes d'âge + +Pour pousser plus loin l'analyse, les individus sont séparés selon leur âge en quatre catégories arbitraires : 18-34 ans, 35-54 ans, 55-64 ans, et les plus de 65 ans. + +```{r} +data$Age[data$Age <= 34] <- "18-34 ans" +data$Age[data$Age > 34 & data$Age <= 54 ] <- "35-54 ans" +data$Age[data$Age > 54 & data$Age <= 64 ] <- "55-64 ans" +data$Age[data$Age > 64] <- "65+ ans" ``` -Et on peut aussi aisément inclure des figures. Par exemple: +Une nouvelle table de contingence est réalisée afin d'obtenir les effectifs dans chaque catégorie d'âge selon leur tabagisme et leur mortalité +```{r} +table2 = data.frame(table(data$Smoker, data$Status, data$Age)) +table2 <- rename(table2, Smoker = Var1, Status = Var2, Classe = Var3) +head(table2) +``` -```{r pressure, echo=FALSE} -plot(pressure) +Le graphique suivant présente la part de décès dans les effectifs des individus selon leur catégorie d'âge et leur tabagisme, et permet de mieux rendre compte de la répartition des effectifs selon ces modalités. +```{r} +ggplot(table2, aes(fill=Status, y=Freq, x=Classe)) + + geom_bar(position="stack", stat="identity")+ + ggtitle("Part de mortalité et effectifs des individus en 1995 selon leur + classe d'âge et leur tabagisme entre 1972-1974")+ + theme(plot.title = element_text(hjust = 0.5, size = 10))+ + xlab("Catégories d'âge")+ylab("Effectifs")+ + facet_wrap(~Smoker,nrow=1)+ + guides(fill=guide_legend(title="Status")) ``` +La fonction ```tauxMorta``` calcule les taux de mortalité de chaque classe d'âge selon leur tabagisme. L'argument ```data``` spécifie le jeu de données utilisé, l'argument ```fumeur``` prend soit ```"Yes"``` ou ````"No"```` selon le tabagisme étudié, et l'argument classe prend, en chaîne de caractère la classe d'âge à étudier. + +```{r} +taux_Morta <- function(data, fumeur, classe){ + nbrSmokerYesNoDead = length(data[["Smoker"]][data["Smoker"] == fumeur + & data["Age"] == classe + & data["Status"] == "Dead"]) + nbrTotal = length(data[["Smoker"]][data["Age"] == classe & data["Smoker"] == fumeur]) + tx = nbrSmokerYesNoDead/nbrTotal + return(tx) +} +``` + + +```{r} +tabac <- c("Yes", "No") +ageClasse <-c("18-34 ans", "35-54 ans", "55-64 ans", "65+ ans") +tauxM <- c() +Class <- c() +Smoker <- c() +for (t in tabac){ + for (c in ageClasse){ + tauxM <- append(tauxM, taux_Morta(data = data, fumeur = t, classe = c)) + Class <- append(Class, c) + Smoker <- append(Smoker, t) + } +} +tauxMortalite <- cbind(Smoker,Class, tauxM) +tauxMortalite <- as.data.frame(tauxMortalite) +``` +On réalise ensuite le graphique en séparant selon le tabagisme. + +```{r} +ggplot(tauxMortalite, aes(fill=tauxM, y=tauxM, x=Class, label =sprintf("%0.2f", tauxM))) + + geom_bar(position="stack", stat="identity", fill = "red")+ + facet_wrap(~Smoker,nrow=1)+ + xlab("Catégories d'âge")+ylab("taux de mortalité")+ + ggtitle("Mortalité des individus en 1995 selon leur + classe d'âge et leur tabagisme entre 1972-1974")+ + theme(plot.title = element_text(hjust = 0.5, size = 10)) +``` + + +Le premier graphique ci-dessus montre une très grande hétérogénéité de l'âge des fumeuses : la plupart des fumeuses sont relativement jeunes (entre 18 et 54 ans) et très peu ont plus 54 ans. On remarque que, naturellement, ce sont ces mêmes femmes (pour la plupart non-fumeuses mais plus agée), qui ont le taux de mortalité le plus élevé. Le second graphique montre que la répartition des taux de mortalité a un profil très similaire chez les fumeuses et les non-fumeuses. Enfin, et bien que ces profils soient très proches, on observe que les taux de mortalité des fumeuses sont légèrement plus élevés que chez les non-fumeuses (en comparant visuellement catégorie à catégorie), ce qui est le contraire de la première observation qui ne prenait en compte que la mortalité tout âge confondu. + + +## Regression logistique + +Enfin, on cherche à vérifier si le tabagisme, ou l'absence de tabagisme, joue effectivement un rôle sur le taux de mortalité. On va pour cela réaliser une régression logistique afin d'évaluer la relation entre le tabagisme et l'âge avec le taux de mortalité des individus. +Ce type de modélisation est adapté à notre cas dans la mesure ou la mortalité suit une loi binomiale, individus "décédé" ou "en vie" (le ```status```), de paramètres inconnus, et cette modélisation nous permettra de nous affranchir du regroupement des individus en classe d'âge arbitraire. On va donc modéliser la réponse ```status``` selon l'âge et le tabagisme des individus. +On créé ainsi un modèle expliquant la mortalité en fonction de l'âge et du tabagisme. On utilise pour cela notre jeu de données initial auquel on fera correspondre les status "décédé" et "vivant" aux valeurs 1 et 0 respectivement. En effet, l'utilisation de la fonction ```glm``` nécéssite que l'ensemble des données soit de type ```numeric```. + +```{r} +data <- read.csv(data_file, header = TRUE) +data$Status <- as.character(data$Status) +data$Status[data$Status == "Dead"] <- "1" +data$Status[data$Status == "Alive"] <- "0" +data$Status <- as.numeric(data$Status) +data$Age <- as.numeric(data$Age) +data$Smoker <- as.character(data$Smoker) +data$Smoker[data$Smoker == "Yes"] <- "1" +data$Smoker[data$Smoker == "No"] <- "0" +data$Smoker <- as.numeric(data$Smoker) +head(data) +``` +On réalise tout d'abord un modèle ```modLogAgeTabac``` avec la fonction ```glm``` prenant en compte l'âge et le tabagisme pour expliquer la mortalité : +```{r} +library(stats) +modLogAgeTabac <- glm(data = data, Status ~ Smoker + Age, family = binomial(logit)) +summary(modLogAgeTabac) +``` +On calcule ensuite les valeurs des status prédits par ce model ```mortalitePred```. Pour cela on génère tout d'abort une séquence d'individus d'âge ```agev``` compris entre 18 et 90 ans par pas de 0.5 an, ainsi qu'une variable ```smokerv``` qui calcule, selon une loi binomial de probabilité 0.5, si l'individu est fumeur ou non. On estime en effet qu'un individu a la même probabilité de fumer que de ne pas fumer. Ces deux variables seront utilisées pour calculer le status des individus selon le modèle. Le status représente l'état "mort" ou "en vie" 20 ans après la première étude, et correspond respectivement à 1 et 0. Ainso plus la valeur prédite sera proche de 1 et plus l'individu aura de chance d'être décédé 20 ans après, et inversement l'individu aura plus de chance d'être en vie si la valeur est proche de 0. +```{r} +agev <- (seq(from=18, to=90, by = .5)) +smokerv <- rbinom(length(agev), 1, prob= 0.5 ) +smokerv <- as.numeric(smokerv) +mortalitePred <- predict(modLogAgeTabac,list(Age=agev, Smoker=smokerv),type="response") +``` +On crée un nouveau *dataframe* contenant le status prédit et l'age des individus. + +```{r} +new.data <- data.frame(mortalitePred, agev) +``` +On réalise ensuite un graphique qui présente le status en ordonnée (1 correspond à "décédé" en 1995, 0 à "en vie" en 1995) et l'âge des individus en 1972-1974. Les points représentent les individus issus de notre jeu de données initial. La couleur bleu correspond aux fumeurs, la couleur noir aux non-fumeurs. Les individus au status prédit par le modèle sont représentés par la courbe. La bande sombre correspond à l'interval de confiance des prédictions (à 95% par défaut). +```{r} +ggplot(data = data, aes(x = Age, y = Status ))+ + geom_point(data = data, aes(x = Age, y = Status, color = Smoker))+ + geom_line(data=new.data,inherit.aes = T, aes(y = mortalitePred, x = agev))+ + guides()+ + stat_smooth(method="glm", method.args=list(family="binomial"), se=TRUE)+ + ggtitle("Mortalité selon le modèle logistique prenant en compte l'âge et le tabagisme")+ + theme(plot.title = element_text(hjust = 0.5, size = 10)) +``` +On réalise ensuite un nouveau modèle ```modLog2``` ne prenant en compte que l'âge pour expliquer la mortalité : +```{r} +modLogAge <- glm(data = data, Status ~ Age, family = binomial(logit)) +summary(modLogAge) +``` +On calcule ensuite les status prédits par ce model ```mortalitePred2```. De même, on génère une séquence d'individus d'âge ```agev2```` entre 18 et 90 ans par pas de 0.5 an. On calcule le status qui leur est associé selon le modèle. +```{r} +agev2 <- (seq(from=18, to=90, by = .5)) +mortalitePred2 <- predict(modLogAge,list(Age = agev2), type="response") +``` + +Le graphique ci-dessous représente le status en fonction de l'âge des individus. Les points représentent les individus issus de notre jeu de données initial. La couleur bleu correspond aux fumeurs, la couleur noir aux non-fumeurs. Les individus au status prédit par le modèle sont représentés par la courbe. La bande sombre correspond à l'interval de confiance des prédictions à 95% . +```{r} +new.data2 <- data.frame(cbind(mortalitePred2, agev2)) +ggplot(data = data, aes(x = Age, y = Status, ))+ + geom_point(data = data, aes(x = Age, y = Status, color = Smoker))+ + geom_line(data=new.data2, aes(y = mortalitePred2, x = agev2))+ + guides()+ + stat_smooth(method="glm", method.args=list(family="binomial"), se=TRUE)+ + ggtitle("Mortalité selon le modèle logistique prenant en compte uniquement l'âge")+ + theme(plot.title = element_text(hjust = 0.5, size = 10)) +``` +On note tout d'abord qu'il semble y avoir une forte relation entre l'age et la mortalité prédite (le status). En effet la p-value issue des 2 modèles pour le paramêtre ```Age``` est très significative tandis que le tabagisme, pour le premier modèle, ne l'est pas du tout. +Cela concorde avec l'observation de la répartition des points des individus issus de nos données. qui sont répartis entre les status 0 et 1 de manière assez uniforme selon leur tabagisme, et il semble y avoir plus d'individus ayant de fortes valeurs d'âge qui prennent le status 1. +Il est intéressant de noter que les legères oscillations de la courbe prédite par le premier model disparaissent quand on supprime le paramètre ```Smoker``` chez le deuxieme model. Ces petites différences de status semble donc être liées à l'influence du tabagisme. + +Enfin, on vérifie la solidité des modèles avec une cross-validation. Pour cela, on compare les status prédits par les modèles avec les status réels des individus. On utilise d'abord le premier modèle, ne prenant en compte que l'âge, et on calcule le nombre de prédiction juste. On fait de même avec le second modèle, prenant en compte l'âge et le tabagisme, et on compare ensuite le pourcentage de prédiction juste des deux modèles. + +```{r} +tablePredModLogAge <- table(actual = data$Status, + predictions = as.numeric(plogis(predict(modLogAge)) > 0.5)) +precisionTablePredModLogAge <-(tablePredModLogAge[1,1] + + tablePredModLogAge[2,2])/sum(tablePredModLogAge) +precisionTablePredModLogAge + +tablePredModLogAgeTabac <- table(actual = data$Status, + predictions = as.numeric(plogis(predict(modLogAgeTabac)) > 0.5)) +precisionTablePredModLogAgeTabac <- (tablePredModLogAgeTabac[1,1] + + tablePredModLogAgeTabac[2,2])/sum(tablePredModLogAgeTabac) +precisionTablePredModLogAgeTabac +``` +Les deux modèles ont chacun 85% de prédictions justes, ainsi le tabagisme ne semble pas jouer un rôle important dans la prédiction de la mortalité des individus, mais c'est bien leur âge qui est déterminant. + +Pour connaitre l'influence relative au tabagisme sur la mortalité de notre modèle, on calcule la différence relative moyenne entre les prédictions de notre modèle avec tabagisme et sans tabagisme. On obtient ainsi la proportion (entre 0 et 1) de la variabilité du status prédits dû à l'influence du tabagisme. +Cette différence est exprimée $$Diff=\frac{1}{N}\sum^{N}_{i=1}\frac {|mortalitéPred_i - mortalitéPred2_i|} {MortalitéPred1_i}$$. +Où $mortalitéPred$ est la prédiction du modèle incluant le tabagisme, $mortalitéPred2$ la prédiction du modèle sans le tabagisme, et $N$ lenombre total d'individus. + +```{r} +mean(abs(mortalitePred - mortalitePred2)/mortalitePred) +``` + +La différence absolue moyenne de mortalité entre les 2 prédictions est comprise entre 8 et 9%. Il semble donc que l'influence du tabagisme sur la mortalité (avec notre jeu de données) soit relativement faible en comparaison de l'âge des individus au moment de la première étude. +Nous avons donc montré que, contrairement à ce que la première observation suggérait, le tabagisme n'augmente pas l'espérence de vie des individus, mais semble plutôt avoir l'effet contraire. En analysant et en séparant les modalitées de notre jeu de données nous avons pu réduire les biais (ici l'âge des individus) qui nous aurait conduit à des observation fausses et des conclusions érronnées. -Vous remarquerez le paramètre `echo = FALSE` qui indique que le code ne doit pas apparaître dans la version finale du document. Nous vous recommandons dans le cadre de ce MOOC de ne pas utiliser ce paramètre car l'objectif est que vos analyses de données soient parfaitement transparentes pour être reproductibles. +On a donc ici un bel exemple du Paradoxe de Simpson, où le manque d'information peut mener à des résultats faux, comme ici quand la population est étudié tout âge confondu plutôt que par tranche d'âge. -Comme les résultats ne sont pas stockés dans les fichiers Rmd, pour faciliter la relecture de vos analyses par d'autres personnes, vous aurez donc intérêt à générer un HTML ou un PDF et à le commiter. -Maintenant, à vous de jouer! Vous pouvez effacer toutes ces informations et les remplacer par votre document computationnel. -- 2.18.1