From 29d6e57bac6d35673e1e8dc9518a4a207e590594 Mon Sep 17 00:00:00 2001 From: 74a384058fe452f679f08e5f0fe15e8b <74a384058fe452f679f08e5f0fe15e8b@app-learninglab.inria.fr> Date: Wed, 24 Feb 2021 07:49:53 +0000 Subject: [PATCH] =?UTF-8?q?R=C3=A9ponse=20=C3=A0=20l'exercice=20"Autour=20?= =?UTF-8?q?du=20paradoxe=20de=20Simpson"?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- module3/exo3/exercice_fr.Rmd | 232 ++++++++++++++++++++++++++++++++--- 1 file changed, 216 insertions(+), 16 deletions(-) diff --git a/module3/exo3/exercice_fr.Rmd b/module3/exo3/exercice_fr.Rmd index 7eece5e..7fce89b 100644 --- a/module3/exo3/exercice_fr.Rmd +++ b/module3/exo3/exercice_fr.Rmd @@ -1,33 +1,233 @@ --- -title: "Votre titre" -author: "Votre nom" -date: "La date du jour" -output: html_document +title: "Autour du Paradoxe de Simpson" +author: "Laura" +date: "18 janvier 2021" +output: + pdf_document: default + html_document: default --- - ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` -## Quelques explications +Informations technique de l'ordinateur réalisant le code +```{r} +sessionInfo() +devtools::session_info() +``` + +## Importer et vérifier les données + +Aller récupérer les données sur le dépot Gitlab +```{r} +data_url="https://gitlab.inria.fr/learninglab/mooc-rr/mooc-rr-ressources/-/raw/master/ +module3/Practical_session/Subject6_smoking.csv" +data=read.csv(data_url) +library(knitr) +kable(head(data)) +``` +On cherche si il y a des données manquantes dans le fichier +```{r} +lignes_na=apply(data,1,function(x) any(is.na(x))) +data[lignes_na,] +``` + + +## Réponses aux exercices +### 1.Exercice 1 +**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** +```{r} +tab_exo1=as.data.frame(table(data$Smoker,data$Status)) +colnames(tab_exo1)=c("Smoker","Status","Number_smoker_status") +library(knitr) +kable(tab_exo1) +``` + +**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 récupère les valeurs du nombres total de fumeuses et non fumeuses +```{r} +tab_smoker_exo1=as.data.frame(table(data$Smoker)) +colnames(tab_smoker_exo1)=c("Smoker","Number_smoker") +library(knitr) +kable(tab_smoker_exo1) +``` +- On ajoute ces valeurs au tableau précédent +```{r} +tab_exo1=merge(tab_exo1,tab_smoker_exo1, by="Smoker", all.x=T) +library(knitr) +kable(tab_exo1) +``` +- On calcul 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). +```{r} +tab_exo1=subset(tab_exo1, Status=="Dead") +tab_exo1$Number_smoker_status=as.numeric(as.character(tab_exo1$Number_smoker_status)) +tab_exo1$Number_smoker=as.numeric(as.character(tab_exo1$Number_smoker)) +tab_exo1$Taux_mortalite= (tab_exo1$Number_smoker_status /tab_exo1$Number_smoker) +library(knitr) +kable(tab_exo1) +``` + +**Vous pourrez proposer une représentation graphique de ces données et calculer des intervalles de confiance si vous le souhaitez.** +```{r} +plot(Taux_mortalite ~ Smoker, data=tab_exo1) +``` + +```{r} +tab_exo1$ICdown=tab_exo1$Taux_mortalite-1.96*sqrt(tab_exo1$Taux_mortalite* + (1-tab_exo1$Taux_mortalite)/tab_exo1$Number_smoker) +tab_exo1$ICup=tab_exo1$Taux_mortalite+1.96*sqrt(tab_exo1$Taux_mortalite* + (1-tab_exo1$Taux_mortalite)/tab_exo1$Number_smoker) +tab_exo1 +``` +On a donc un taux de mortalité pour les fumeuses de 0.24 avec un IC [0.20;0.27] +et pour les non fumeuse 0.31 avec un IC [0.28;0.35] + +**En quoi ce résultat est-il surprenant?** +On a un taux de mortalité plus important pour les femmes non fumeuses. + +### 2.Exercice 2 +**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.** +Dans la consigne le 34 est pris 2 fois, j'ai fais le choix de ne le prendre que dans la première classe + +- On définit les différentes classes d'âge +```{r} +data$classe_age=NA +data$classe_age[data$Age>=18 & data$Age<35]="1" +data$classe_age[data$Age>=35 & data$Age<64]="2" +data$classe_age[data$Age>=65 ]="3" +data$classe_age=as.factor(data$classe_age) +levels(data$classe_age) +``` + +- On représente 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 et de leur classe d'âge +```{r} +tab_exo2=as.data.frame(table(data$Smoker,data$Status,data$classe_age)) +colnames(tab_exo2)=c("Smoker","Status","Class_Age","Number_smoker_status_classage") +library(knitr) +kable(tab_exo2) +``` + +- On récupère les valeurs du nombre total de fumeuses et non fumeuses par classe d'âge +```{r} +tab_smoker_classage_exo2=as.data.frame(table(data$Smoker,data$classe_age)) +colnames(tab_smoker_classage_exo2)=c("Smoker","Class_Age","Number_smoker_classage") +library(knitr) +kable(tab_smoker_classage_exo2) +``` + +- On ajoute ces valeurs au tableau précédent +```{r} +tab_exo2=merge(tab_exo2,tab_smoker_classage_exo2, by=c("Smoker", "Class_Age"), all.x=T) +library(knitr) +kable(tab_exo2) +``` +- On calcul 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). +```{r} +tab_exo2=subset(tab_exo2, Status=="Dead") +tab_exo2$Number_smoker_status_classage=as.numeric(as.character + (tab_exo2$Number_smoker_status_classage)) +tab_exo2$Number_smoker_classage=as.numeric(as.character + (tab_exo2$Number_smoker_classage)) +tab_exo2$Taux_mortalite= tab_exo2$Number_smoker_status_classage / + tab_exo2$Number_smoker_classage +library(knitr) +kable(tab_exo2) +``` -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 . +**Vous pourrez proposer une représentation graphique de ces données et calculer des intervalles de confiance si vous le souhaitez.** +```{r} +boxplot(Taux_mortalite ~ Smoker*Class_Age, data=tab_exo2) +``` + +```{r} +tab_exo2$ICdown=tab_exo2$Taux_mortalite-1.96*sqrt(tab_exo2$Taux_mortalite* + (1-tab_exo2$Taux_mortalite)/tab_exo2$Number_smoker_classage) +tab_exo2$ICup=tab_exo2$Taux_mortalite+1.96*sqrt(tab_exo2$Taux_mortalite* + (1-tab_exo2$Taux_mortalite)/tab_exo2$Number_smoker_classage) +tab_exo2 +``` +**En quoi ce résultat est-il surprenant?** +Dans ce cas là, on voit que l'on a un plus fort taux de mortalité pour les femmes fumeuses dans la classe d'âge 2 et une tendance similaire pour les femmes de la classe d'âge 1. -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: +**Arrivez-vous à expliquer ce paradoxe? De même, vous pourrez proposer une représentation graphique de ces données pour étayer vos explications.** -```{r cars} -summary(cars) +On peut voir qu'il y a beaucoup plus de femmes de la classe d'âge 3 dans les non fumeuses, donc une mortalité qui semble plus importante pour les non fumeuse car cette classe d'âge et celle ou il y a un plus fort taux de mortalité +```{r} +tab_smoker_classage_exo2=as.data.frame(table(data$Smoker,data$classe_age)) +colnames(tab_smoker_classage_exo2)=c("Smoker","Class_Age","Number_smoker_classage") +library(knitr) +kable(tab_smoker_classage_exo2) ``` -Et on peut aussi aisément inclure des figures. Par exemple: +### 3.Exercice 3 +**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 introduit une variable Death valant 1 ou 0 pour indiquer si l'individu est décédé durant la période de 20 ans + +```{r} +data$Death=NA +data$Death[data$Status=="Alive"]=0 +data$Death[data$Status=="Dead"]=1 +library(knitr) +kable(head(data)) +``` -```{r pressure, echo=FALSE} -plot(pressure) +- On étudie 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. + +Pour cela on utilise un modèle linéaire généralisé avec une distribution binomiale, comprenant l'intéraction des effets Age et Smoker +```{r} +hist(data$Death) +m1=glm(Death ~ Age * Smoker, data=data, family=binomial(logit)) +summary(m1) +``` + +- On plot les prédictions de ce modèle en indiquant les intervals de confiance à 95% +```{r} +data$Smoker=as.factor(data$Smoker) +grid = expand.grid( Smoker=levels(data$Smoker),Age=seq(min(data$Age), max(data$Age),1)) +library(AICcmodavg) +pred <- predict(m1,grid, se.fit=T,type="response") +grid$y=pred$fit +grid$low=pred$fit-1.96*pred$se.fit +grid$up=pred$fit+1.96*pred$se.fit +head(grid) +Nonfumeuse=subset(grid, Smoker=="No") +Fumeuse=subset(grid, Smoker=="Yes") +par(mar=c(5,5.5,1,2)) +plot(y~Age,xlab="Age",ylab="Death" , + data=Nonfumeuse,axes=F, cex.lab=2, col="blue",type="l",lwd=2) +lines(low~Age, data=Nonfumeuse, col="blue",lwd=1,lty=2) +lines(up~Age, data=Nonfumeuse, col="blue",lwd=1,lty=2) + +lines(y~Age, data=Fumeuse, col="red",lwd=2) +lines(low~Age, data=Fumeuse, col="red",lwd=1,lty=2) +lines(up~Age, data=Fumeuse, col="red",lwd=1,lty=2) +box(bty="l") +axis(1,cex.axis=1.5) +axis(2,cex.axis=1.5) +legend("bottomright",bty="n", legend=c("Fumeuse","Non fumeuse"), + col=c("red","blue"),lwd=2,cex=1) ``` -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. +**Ces régressions vous permettent-elles de conclure sur la nocivité du tabagisme ?** + +Non, on voit que d'une manière générale la probabilité de mourir augmente avec l'âge, mais on trouve de faibles différences entre les courbes des fumeuses et des non fumeuses, avec des intervals de confiance trés chevauchant hormis dans la tranche d'âge autour de 40 ans. -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. +Pour confirmer cela on peut faire une sélection de modèle sur le critère AICc, en comparant le modèle avec l'effet smoker et celui sans: +```{r} +m0=glm(Death ~ Age , data=data, family=binomial) +m1=glm(Death ~ Age + Smoker, data=data, family=binomial) +m2=glm(Death ~ Age * Smoker, data=data, family=binomial) +library(AICcmodavg) +Cands <- list(m0,m1,m2) +Model.names <- c("Age","Age+smoker","Age*smoker") +tab_aic2=aictab(cand.set = Cands, modnames = Model.names, sort = TRUE, second.ord = T) +library(knitr) +kable(tab_aic2) +``` +On voit bien que malgré que le modèle avec les effet Age x smoker soit celui avec un plus faible AIC, il y a un faible soutien au modèle comprenant l'interaction Age*smoker, avec un delta AICc <2 par rapport au modèle uniquement avec Age, il est donc difficle de conclure sur la nocivité du tabagisme. +Si l'on utilise le principe de parcimonie on retiendrait le modèle le plus simple avec un delta AICc<2, donc le modèle uniquement avec âge dans ce cas là. -Maintenant, à vous de jouer! Vous pouvez effacer toutes ces informations et les remplacer par votre document computationnel. -- 2.18.1