You need to sign in or sign up before continuing.
Commit 3cb94dcc authored by Clair Ch's avatar Clair Ch

sujet 6 : partie régression linéaire

parent c6ed171d
---
title: "Sujet 6 : Autour du Paradoxe de Simpson"
title: 'Sujet 6 : Autour du Paradoxe de Simpson'
author: "Clair Ch"
date: "20 avril 2020"
output: html_document
output:
pdf_document: default
html_document: default
---
......@@ -10,7 +12,7 @@ output: html_document
knitr::opts_chunk$set(echo = TRUE)
```
# Instructions
## Instructions
Voici les instructions du Sujet 6 : Autour du Paradoxe de Simpson
......@@ -28,16 +30,19 @@ Voici les instructions du Sujet 6 : Autour du Paradoxe de Simpson
4. Déposez votre étude dans FUN
# Préparation des données
## Préparation des données
## Téléchargement
#### Téléchargement
Les données sont disponibles sur le [Gitlab](https://gitlab.inria.fr/learninglab/mooc-rr/mooc-rr-ressources/blob/master/module3/Practical_session/Subject6_smoking.csv) du MOOC ["Recherche Reproductible : principes méthodologiques pour une science transparente"](https://www.fun-mooc.fr/courses/course-v1:inria+41016+session01bis/about) de l'Inria. On peut les récupérer au format .csv à cette adresse :
```{r}
data_url<-"https://gitlab.inria.fr/learninglab/mooc-rr/mooc-rr-ressources/-/raw/master/module3/Practical_session/Subject6_smoking.csv?inline=false"
data_url<-"https://gitlab.inria.fr/learninglab/mooc-rr/mooc-rr-ressources/-/raw/master/
module3/Practical_session/Subject6_smoking.csv?inline=false"
```
Elles regroupent, par femme, les habitudes de tabagisme (fumeuse/non fumeuse), le statut lors de la deuxième étude de 1995 (vivante ou morte) et l'âge lors de la première étude (effectuée entre 1972 et 1974).
Pour nous protéger contre une éventuelle disparition ou modification du serveur du Gitlab du MOOC, nous faisons une copie locale de ce jeux de données que nous préservons avec notre analyse. Il est inutile et même risqué 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.
......@@ -48,7 +53,8 @@ if (!file.exists(data_file)) {
}
```
## Lecture
#### Lecture
```{r}
data<-read.csv(data_file)
......@@ -63,6 +69,7 @@ tail(data)
Nous avons les données pour 1314 femmes, avec dans l'ordre : si la personne fume ou non (colonne `Smoker` : `Yes` ou `No`), si elle est morte ou vivante au moment de la deuxième étude en 1995 (colonne `Status` : `Dead` ou `Alive`) et son âge lors de la première étude faite entre 1972 et 1974 (colonne `Age`).
Y a-t'il des points manquants dans nos données ?
```{r}
......@@ -72,6 +79,7 @@ data[na_records,]
Aucune donnée manquante, parfait !
Vérifions la classe de chaque colonne
```{r}
......@@ -79,15 +87,15 @@ class(data$Smoker)
class(data$Status)
class(data$Age)
```
La données de la colonne `Age` sont bien de classe numeric
Les données des colonnes `Smoker`et `Status`sont de classe factor, ça posera peut-être problème par la suite ?
Les données de la colonne `Age` sont bien de classe numeric, les autres colonnes de type factor.
# Décès et tabagisme
## Décès et tabagisme
Nous voulons étudier le taux de mortalité pour chaque groupe (fumeuses/non fumeuses)
## Inspection des données
#### Inspection des données
Faisons un tableau représentant le nombre de femmes par habitude de tabagisme (fumeuses/non fumeuses) et statut (vivantes/mortes)
......@@ -96,12 +104,14 @@ t<-table(data$Smoker,data$Status)
t
```
Vérifions que la somme de chaque catégorie correspond bien au nombre total de femmes de l'étude
```{r}
sum(t)==nrow(data)
```
Nous pouvons voir qu'il y a plus de femmes vivantes que de femmes mortes et plus de femmes qui ne fument pas :
```{r}
......@@ -110,6 +120,7 @@ alive
dead<-sum(t[,2]) #femmes mortes
dead
```
Il y a également plus de femmes non fumeuses que fumeuses :
```{r}
......@@ -119,56 +130,388 @@ smoke<-sum(t[2,]) #femmes fumeuses
smoke
```
## Calcul du taux de mortalité
#### Calcul du taux de mortalité
Calculons dans chaque groupe (fumeuses/non fumeuses) le taux de mortalité (rapport entre le nombre de femmes décédées dans un groupe et le nombre total de femmes dans ce groupe)
taux de mortalité pour les femmes non fumeuses :
```{r}
nsmoke_morta<-t[1,2]/nsmoke
nsmoke_morta
```
taux de mortalité pour les femmes fumeuses :
```{r}
smoke_morta<-t[2,2]/smoke
smoke_morta
```
Les femmes non fumeuses présentent un taux de mortalité plus élevé !
De combien est-il plus élevé ?
```{r}
nsmoke_morta/smoke_morta
```
## Calcul de l'intervalle de confiance à 95% du taux de mortalité
Arrondi au centième, les femmes non fumeuses présent donc un taux de mortalité `r round(nsmoke_morta/smoke_morta,digit=2)` fois plus élevé que les femmes fumeuses dans cette étude. C'est surprenant car contre-intruitif !
On peut faire un test statistique pour avoir une idée si cette différence n'est pas due au hasard (test Chi-square de Pearson sur les proportions) :
```{r}
chisq.test(t)
```
La p-value est en dessous de 0.05, il y a donc moins de 5% de chance que cette différence du taux de mortalité soit due au hasard. Mais est-elle vraiment due aux habitudes de tabagisme ?
#### Calcul de l'intervalle de confiance à 95% du taux de mortalité
Pour calculer les intervalles de confiance à 95% du taux de mortalité par catégorie (fumeuses/non fumeuses), nous allons suivre les instructions de l'article [How to Determine the Confidence Interval for a Population Proportion](https://www.dummies.com/education/math/statistics/how-to-determine-the-confidence-interval-for-a-population-proportion/) du site [dummies.com](https://www.dummmies.com)
mettre formule mathématique
La formule donnée pour calculé l'intervalle de confiance à 95% est la suivante :
$$1,96\sqrt{\frac{\beta(1-\beta)}{n}}$$
$\beta$ représente ici le taux de mortalité des femmes fumeuses ou non fumeuses et $n$ le nombre total de femmes fumeuses ou non fumeuses
Créons la fonction correspondante :
```{r}
CI95<-function(morta,n){
1.96*sqrt((morta*(1-morta))/n)
}
```
Calculons les intervalles de confiance pour les femmes non fumeuses et les femmes fumeuses :
```{r}
nsmoke_CI<-CI95(nsmoke_morta,nsmoke)
smoke_CI<-CI95(smoke_morta,smoke)
```
Arrondi au centième, les femmes fumeuses ont un taux de mortalité de `r round(smoke_morta,digit=2)`$\pm$ `r round(smoke_CI,digit=2)`
Arrondi au centième, les femmes non fumeuses ont un taux de mortalité de `r round(nsmoke_morta,digit=2)`$\pm$ `r round(nsmoke_CI,digit=2)`
#### Représentation graphique
Pour faire les graphiques, nous utiliserons la librairie `ggplot2`.
Ce code permet d'installer la libraire si nécessaire et de la loader.
```{r}
if(!require(ggplot2)){
install.packages("ggplot2")
library(ggplot2)
}
```
Pour faire le graphe, il faut d'abord mettre les données de taux de mortalité et d'intervalle de confiance dans une data frame, j'en profite pour les passer en pourcentage :
```{r}
df_morta<-data.frame(Smoker=c("No","Yes"), Mortality=c(nsmoke_morta*100,smoke_morta*100), CI95=c(nsmoke_CI*100,smoke_CI*100))
df_morta
```
Avec l'aide d'un [tutoriel pour faire
des barplots](http://www.sthda.com/french/wiki/ggplot2-barplots-guide-de-demarrage-rapide-logiciel-r-et-visualisation-de-donnees#barplot-du-comptage-des-observations), faisons un graphe représentant le taux de mortalité en fonction des habitudes de tabagisme :
```{r}
plot_morta<-ggplot(df_morta,aes(x=Smoker,y=Mortality))+
geom_bar(stat="identity",fill="steelblue",width=0.5)+
theme_minimal()+
ylim(c(0,50))+ #limite de l'axe y
geom_errorbar(aes(ymin=Mortality-CI95,ymax=Mortality+CI95),width=.2,position=position_dodge(.9))+ #barres d'erreurs
labs(title="Taux de Mortalité en fonction des habitudes de tabagisme", subtitle="Taux de Mortalité en pourcentage +/- intervalle de confiance à 95%")
plot_morta
```
Les femmes non fumeuses semblent avoir un taux de mortalité plus élevé que les femmes fumeuses. Si on considère ce résultat comme tel, on serait tenté de conclure que la cigarette est bonne pour la santé !
## Décés et Tabagisme en fonction de la classe d'âge
#### Ajout de la tranche d'âge
Rajoutons une colonne `Tranche` dans le fichier `data`donnant la tranche d'âge (`18-34`: de 18 ans inclus à 35 ans non inclus ; `35-54`: de 35 ans inclus à 55 ans non inclus ; `55-64`: de 55 ans inclus à 65 non inclus ; `>65`: plus de 65 ans inclus), mettons le message `ERROR`si une tranche d'âge n'est pas trouvée
```{r}
for(i in 1:nrow(data)){
if(data$Age[i]>=65){data$Tranche[i]<-">65"
}else if(data$Age[i]>=55){data$Tranche[i]<-"55-64"
}else if(data$Age[i]>=35){data$Tranche[i]<-"35-54"
}else if(data$Age[i]>=18){data$Tranche[i]<-"18-34"
}else(data$Tranche[i]<-"ERROR")
}
```
Vérifions que ça ait bien marché :
```{r}
head(data)
tail(data)
```
Est-ce qu'il y a eu des erreurs ?
```{r}
sum(data$Tranche=="ERROR")
```
Aucune erreur, tout va bien !
#### Calcul du taux de mortalité par tranche d'âge et tabagisme
Faisons une fonction permettant de calculer le nombre de femmes par tranche d'âge selon leur habitude de tabagisme (fumeuse/non fumeuse) et éventuellement leur statut (dead/alive)
```{r}
comptage<-function(smoker,statut,tranche){
if(statut=="NULL")
{sum(data$Smoker==smoker&data$Tranche==tranche)
}else{sum(data$Smoker==smoker&data$Status==statut&data$Tranche==tranche)}
}
```
Calculons le nombre de femmes totales et mortes par tranche d'âge et catégories
```{r}
tranches<-c("18-34","35-54","55-64",">65")
nsmoke_age<-mapply(comptage,smoker="No",statut="NULL",tranche=tranches) #femmes non fumeuses
smoke_age<-mapply(comptage,smoker="Yes",statut="NULL",tranche=tranches) #femmes fumeuses
nsmoke_age_dead<-mapply(comptage,smoker="No",statut="Dead",tranche=tranches) #femmes non fumeuses mortes
smoke_age_dead<-mapply(comptage,smoker="Yes",statut="Dead",tranche=tranches) #femmes fumeuses mortes
#nommer les tranches d'âge dans les vecteurs :
names(nsmoke_age)<-tranches
names(smoke_age)<-tranches
names(nsmoke_age_dead)<-tranches
names(smoke_age_dead)<-tranches
```
Vérifions que les nombres de femmes fumeuses et non fumeuses correspondent bien au nombre de femmes totales
```{r}
sum(smoke_age,nsmoke_age)==nrow(data)
```
Vérifions également que le nombre de femmes mortes fumeuses et non fumeuses correspondent bien au nombre de femmes totales mortes :
```{r}
sum(smoke_age_dead,nsmoke_age_dead)==nrow(data$Status=="Dead")
```
Nous pouvons maintenant calculer le taux de mortalité pour les femmes non fumeuses par tranche d'âge
```{r}
nsmoke_age_morta<-nsmoke_age_dead/nsmoke_age
nsmoke_age_morta
```
Et pour les femmes fumeuses par tranche d'âge :
```{r}
smoke_age_morta<-smoke_age_dead/smoke_age
smoke_age_morta
```
Mettons les taux de mortalité en pourcentage au centième près selon les tranches d'âges et habitude de tabagisme dans une data frame pour pouvoir les comparer :
```{r}
df_morta_age<-data.frame(Age=tranches,"Mortality non Smoker"=round(nsmoke_age_morta*100,digit=2),"Mortality Smoker"=round(smoke_age_morta*100,digit=2))
df_morta_age
```
On peut voir tout d'abord que de façon générale, le taux de mortalité est plus important pour les catégories d'âges les plus "avancées", ce qui est un résultat attendu. Ici la mortalité des femmes qui fument semble plus importante que celle des femmes qui ne fument pas dans toutes les tranches d'âge (à part celle des femmes de plus de 65 ans). On a donc une différence quand on regarde le taux de mortalité total et le taux de mortalité par tranche d'âge. Le taux de mortalité plus élevé chez les femmes non fumeuses n'était peut-être pas du au tabagisme...
#### Calcul de l'intervalle de confiance à 95%
De la même manière que pour le taux de mortalité global par tranche d'âge, nous calculons l'intervalle de confiance à 95% en prenant en compte la tranche d'âge :
```{r}
nsmoke_age_IC<-mapply(CI95, morta=nsmoke_age_morta,n=nsmoke_age)
smoke_age_IC<-mapply(CI95,morta=smoke_age_morta,n=nsmoke_age)
```
#### Graphe
Mettons au propres les données pour faire le graphe: la 1ère colonne repésente l'habitude de tabagisme, la deuxième la tranche d'âge, la troisième le taux de mortalité (en pourcentage au centième près), la quatrième l'intervalle de confiance à 95% (en pourcentage au centième près)
```{r}
df_morta_age2<-data.frame(Smoker=c(rep("No",4),rep("Yes",4)),Age=rep(c("18-34","35-54","55-64",">65"),2), Mortality=c(nsmoke_age_morta,smoke_age_morta),CI95=c(nsmoke_age_IC,smoke_age_IC))
df_morta_age2$CI95<-round(df_morta_age2$CI95*100,digits=2)
df_morta_age2$Mortality<-round(df_morta_age2$Mortality*100,digits=2)
df_morta_age2
```
Pour que les tranches d'âges soient dans le bon ordre dans le graphe, nous ordonnons les niveaux de la colonne `Age`
```{r}
df_morta_age2$Age<-ordered(df_morta_age2$Age, levels=c("18-34","35-54","55-64",">65"))
```
Faisons maintenant le graphe, toujours grâce à la libraire `ggplot2` :
```{r}
if(!require(ggplot2)){
install.packages("ggplot2")
library(ggplot2)
}
```
```{r}
plot_morta_age<-ggplot(df_morta_age2,aes(x=Age,y=Mortality,fill=Smoker))+
geom_bar(stat="identity",position="dodge")+
scale_fill_brewer(palette="Paired")+
theme_minimal()+
geom_errorbar(aes(ymin=Mortality-CI95, ymax=Mortality+CI95), width=.2, position=position_dodge(.9))+ #barres d'erreur
labs(title="Taux de Mortalité par catégorie d'âge en fonction des habitudes de tabagisme", subtitle="Taux de Mortalité en pourcentage +/- intervalle de confiance à 95%")
plot_morta_age
```
Ici il apparait que la mortalité est plus importante pour les femmes fumeuses que non fumeuses, la cigarette ne semble au final pas bon pour la santé !
#### Paradoxe
Le paradoxe entre le premier résultat (si on regarde le taux de mortalité de façon globale, il est plus élevé pour les femmes non fumeuses) et ce dernier (si on classe les femmes par catégories d'âges, on a un taux de mortalité plus élevé pour les femmes fumeuses) pourrait provenir du fait que les catégories d'âges ne sont pas représentées de manière égales entre les femmes fumeuses et non fumeuses : les femmes "jeunes" auraient plus tendance à fumer, et au contraire les femmes "âgées" auraient tendance à moins fumer et seraient plus représentées dans le groupe des non fumeuses. Ainsi le taux de mortalité "naturelle" (de vieillesse) est globalement plus important pour le groupe des non fumeuses.
Nous pouvons vérifier cette hypothèse en analysant la distribution des âges par habitude de tabagisme :
```{r}
summary(data$Age[data$Smoker=="No"])
summary(data$Age[data$Smoker=="Yes"])
```
On voit que l'étendu des âges est sensiblement la même entre les femmes non fumeuses et fumeuses (de 18 à environ 89-90 ans). Par contre l'âge médian est plus élevé chez les femmes non fumeuses que non fumeuses (48.40 ans contre 43.10 ans). C'est en concordance avec notre théorie.
Comparons maintenant les histogrammes (distribution des âges) pour avoir une meilleure idée de la distribution des âges :
```{r}
par(mfrow=c(1,2)) #pour avoir les histogrammes côte à côte
hist(data$Age[data$Smoker=="No"],main="Femmes non fumeuses",xlim=c(18,90),ylim=c(0,80),xlab="Age")
hist(data$Age[data$Smoker=="Yes"], main="Femmes fumeuses",xlim=c(18,90),ylim=c(0,80),xlab="Age")
```
En comparant ces histogrammes, on voit bien que la distribution des âges n'est pas du tout la même entre les femmes fumeuses et non fumeuses : il y a très peu de femmes "âgées" dans les femmes fumeuses. Ainsi beaucoup moins de femmes meurent de mort "naturelle" (vieillesse) dans le groupe des femmes fumeuses, ce qui abaisse le taux de mortalité global. Pour savoir si la cigarette a en effet un effet sur le taux de mortalité, il faut prendre en compte cette variable âge pour s'affranchir du taux de mortalité du à la vieillesse.
## Régression logistique
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.
La régression logistique permet d'évaluer et caractériser les relations entre deux variables, ici l'âge et la mortalité. On peut trouver une explication dans ce [document](https://gitlab.inria.fr/learninglab/mooc-rr/mooc-rr-ressources/blob/master/module2/exo5/challenger_fr.pdf) du MOOC "Recherche reproductible" et une mise en pratique dans [celui-ci](https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/blob/master/challenger.pdf)
#### Création de la variable Death
Ajoutons une variable `Death`prenant la valeur de 0 si la personne est décédée et 1 si la personne est vivante
```{r}
for(i in 1:nrow(data)){
if(data$Status[i]=="Dead")
{data$Death[i]<-0}
else(data$Death[i]<-1)
}
```
On vérifie que ça a bien fonctionné
```{r}
head(data)
```
#### Calcul de la régression logistique
Faisons une régression logistique de la mortalité en fonction de l'âge selon que la personne fume ou non :
```{r}
logistic_reg<-glm(Death~Age+Smoker,data=data,family=binomial)
summary(logistic_reg)
```
On voit qu'il y a un effet négatif de l'âge sur la mortalité (estimate de -0.099+/-0.005) qui est signicatif (p<0.05).
Il y a aussi un effet négatif de l'habitude de tabagisme (estimate de -0.27+/-0.14) mais qui n'est pas significatif (p>0.05).
Cette régression logistique ne montre aucune évidence d'un effet du tabagisme sur la mortalité dans ce jeu de données.
#### Graphes
Pour faire le graphe, nous utiliserons la librairie `ggplot2`
```{r}
nsmoke_CI<-1.96*sqrt((nsmoke_morta*(1-nsmoke_morta))/nsmoke)
smoke_CI<-1.96*sqrt((smoke_morta*(1-smoke_morta))/smoke)
if(!require(ggplot2)){
install.packages("ggplot2")
library(ggplot2)
}
```
Faisons le graphe de la régression logistique de la mortalité (Death vaut 0 is la personne est décédée et 1 si la personne est morte) en fonction des habitudes de tabagisme et de l'âge, avec en grisé l'intervalle de confiance à 95%
```{r}
ggplot(data,aes(x=Age,y=Death,color=Smoker)) + geom_point(alpha=.3,size=1) +
theme_bw() +
geom_smooth(method = "glm",
method.args = list(family = "binomial"))+
labs(title="Logistic regression")
```
Sur ce graphe, il ne semble pas globalement que le fait de fumer augmente ou non la probabilité d'être décedée 20 ans plus tard. Cependant la courbe des femmes non fumeuses semble faire un coude et être légèrement au dessus de celles des fumeuses pour un âge entre 25 et 65 ans. Re-faire une analyse sur cette tranche d'âge serait intéressant. Le tabagisme n'a peut-être aucun effet dans nos données sur le femmes les plus jeunes et les plus âgées (les plus âgées vont mourir de vieillesse, et peut-être que les femmes les plus jeunes ne fument pas depuis assez longtemps pour que ça ait un effet évident ou qu'elles sont plus protégées d'un effet néfaste du tabac).
## Représentation graphique
#### Régression logistique sur les femmes entre 25 et 65 ans
Pour faire les graphiques, nous utiliserons la librairie `ggplot2`
Pour cette partie nous utiliserons la librarie `dplyr
```{r}
library(ggplot2)
if(!require(dplyr)){
install.packages("dplyr")
library(dplyr)
}
```
Créons un deuxième jeu de données ne comprenant que les femmes entre 25 et 65 ans.
```{r}
data2<-filter(data,Age>25&Age<65)
```
# Truc
Faisons une régression logistique sur ces données :
```{r}
hist(data$Age)
summary(data$Age)
logistic_reg2<-glm(Death~Age+Smoker,data=data2,family=binomial)
summary(logistic_reg2)
```
Ici à la fois l'âge et l'habitude de tabagisme ont un effet négatif significatif. Il semble donc que le fait de fumer affecte la probabilité de décéder pour les femmes qui avaient entre 25 et 65 ans.
Faisons le nouveau graphe :
```{r}
ggplot(data2,aes(x=Age,y=Death,color=Smoker)) + geom_point(alpha=.3,size=1) +
theme_bw() +
geom_smooth(method = "glm",
method.args = list(family = "binomial"))+
labs(title="Logistic regression 2")
```
Sur ce nouveau graphe il apparait que le fait de fumer augmente légèrement la probabilité d'être décédée 20 ans plus tard pour les personnes qui avaient entre 25 et 65 ans.
## Conclusion
Les femmes fumeuses qui avaient entre 25 et 65 ans au moment de la première étude semblent avoir une probabilité plus forte de mourir que les femmes non fumeuses. Le fait de fumer augmenterait donc la probabilité de mourir dans les 20 ans si on est d'âge "jeune-moyen". Avec ce jeu de données, il semble donc que le fait de fumer pour les femmes "âgées" ou "très jeunes" n'ait pas d'incidence sur la probabilité de mourir 20 ans plus tard. On peut faire l'interprétation que si la femme est âgée, 20 ans plus tard elle va très probalement mourir qu'elle fume ou non (pour les fumeuses à cause du tabagisme et/ou vieillesse et pour les non fumeuses de vieillesse), et pour les plus jeunes soit elles ne fument pas depuis assez longtemps pour que ça ait une incidence sur leur probabilité de mourir, soit elles sont "protégées"" par leur "jeunesse". Il serait intéressant de faire une étude à intervalles de temps plus rapproché (par exemple tous les 5 ans), pour limiter l'effet de la mort naturelle sur les personnes les plus âgées et en prenant en compte depuis combien de temps la personne fume pour voir si cela a un effet sur la probabilité de mourir.
Le tabagisme semble donc avoir un effet néfaste sur la mortalité qui est dépendant de l'âge.
\ No newline at end of file
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment