[Extrait du Sujet 6: Autour du Paradoxe de Simpson](https://www.fun-mooc.fr/courses/course-v1:inria+41016+self-paced/courseware/5b932aa591d245d48d8943385cb3120a/57c96f2c7f7b42018eaac3e6b34546f4/)
"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."
Pour nous protéger contre une éventuelle disparition ou modification du fichier, nous faisons une copie locale de ce jeu 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.
```{r}
data_file = "Subject6_smoking.csv"
if (!file.exists(data_file)) {
download.file(data_url, data_file, method="auto")
}
```
Voici l'explication des colonnes données décrit dans le [sujet 6](https://www.fun-mooc.fr/courses/course-v1:inria+41016+self-paced/courseware/5b932aa591d245d48d8943385cb3120a/57c96f2c7f7b42018eaac3e6b34546f4/)
Chaque ligne représente une femme.
Nom de colonne | Decription de colonne
------------- | -------------
`Smoke` | La femme interrogée est fumeuse ou non
`Statut` | La femme interrogée est vivante ou morte au moment de la seconde étude
`Age` |Age de la femme interrogée au moment du premier sondage
# 2. Lecture des données
```{r}
data=read.csv(data_file, header=T)
```
1. Observation de l'aspect des données
Pour visualiser comment la lecture du jeu de données a été effectuée, on choisit d'observer le début et la fin de notre jeu de données *via* les fontions head et tail
2. Y a-t-il des points manquants dans nos données ?
```{r}
na_records = apply(data, 1, function (x) any(is.na(x)))
data[na_records,]
```
Aucune ligne ne présente de données manquantes, il n'y aura pas besoin de faire attention à ce point dans la suite des analyses
3. Vérification de la nature des variables
On vérifie maintenant que R classifie correctement nos variables.
On attend que:
`Smoker` et `Status`soient considérés comme des facteurs
`Age` soit considéré comme un numérique
```{r}
str(data)
```
Tout est correct.
__Conclusion__
Le jeu de données est correctement lu par R et le jeu de données ne présente pas de données manquantes.
Nous pouvons donc commencer à répondre au question de l'exercice.
# 3. Question 1
## 3.1 Enoncé
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. 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). Vous pourrez proposer une représentation graphique de ces données et calculer les intervalles de confiance. En quoi ce résultat est-il surprenant ?
## 3.2 Résolusion de la question 1
### 3.2.1 Tableau du nombre total de femmes vivantes et décédées sur la période
### en fonction de leur habitude de tabagisme
On crée une table de contigence
```{r}
statut=table(data$Status, data$Smoker)
statut
```
### 3.2.2 Calculs des taux de mortalité
### en fonction des habitudes de tabagisme
On exprime le taux de mortalité par habitude de tabagisme en appliquant pour un
groupe la formule suivante:
$mortalité = Ndécédée / (Nvivante + Ndécédée)$
N = nombre de femmes du groupe répondant au statut mentionné (décédée ou vivante)
Pour obtenir ces taux de mortalité, nous allons calculer les pourcentage globaux pour la table de contingence précédemment créée en fonction de l'habitude de tabagisme
```{r}
prop_statut= prop.table(statut,2)
prop_statut
```
L'argument 2 nous permet de specifier que nous souhaitons obtenir un calcul des fréquences par colonne.
Les colonnes correspondant à l'habitude de tabagisme, c'est ce que nous recherchons
__Conclusion__
__Le taux de mortalité chez les non fumeuses est de 31.4%__
__Le taux de motalité chez les fumeuses est de 23.9%__
### 3.2.3 Représentation graphique des données
Ici on choisit de visualiser les taux de mortalité en fonction de l'habitude de tabagisme à l'aide d'un graphique en barres
title(xlab="Habitude de Tabagisme", ylab="Taux de mortalité", main="Taux de mortalité en fonction de l'habitude de tabagisme")
```
## Quelques explications
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 <http://rmarkdown.rstudio.com>.
### 3.2.4 Calculs des intervalles de confiances (IC)
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:
On commence par récupérer le nombre d'observations pour les femmes non fumeuses pour tous les statuts
```{r}
nb.No <- as.vector(statut[,1])
```
On fait de même pour les fumeuses
```{r}
nb.Yes <- as.vector(statut[,2])
```
On récupère ensuite le nombre total d'observations en fonction des habitudes de tabagisme
```{r}
#Pour les non fumeuses
nbT.No<-sum(nb.No)
#Pour les fumeuse
nbT.Yes<-sum(nb.Yes)
```
On calcule nos IC grace à la fonction `binom.confit` associé au package `binom`. On cherche donc le package `binom`
__On retrouve bien un taux moyen de mortalité = 0.2388__
__IC pour le taux de mortalité des non fumeuses = [0.205;0.276]__
### 3.2.5 En quoi ce résultat est-il surprenant ?
Ce résultat est surprenant car on s'attendrait à voir que, dans ce contexte de travaux sur les maladies thyroïdiennes et cardiaques le taux de mortalité est plus importante chez les femmes fumeuses. C'est l'inverse qui ressort ici: le taux de mortalité est plus important chez les femmes n'ayant jamais fumé.
# 4. Question 2
## 4.1 Enoncé
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. 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.
## 4.2 Réponse à la question 2
### 4.2.1 Ajout de la variable âge
Ici on propose d'ajouter à notre jeu de données la variable âge par classes.
4 classes ont été définit:
1. 18-34 ans
2. 34-54 ans
3. 55-64 ans
4. plus de 65 ans
__Remarque__:
Certaines classes ne sont pas continues (_i.e._: 34-54 et 55-64) ceci s'explique par le fait qu'aucune femme interrogée n'a entre [54.1-54.9] ans. De même il n'y a pas de femme interrogée ayant entre 64.0 et 64.9 ans.
Néanmoins ces classes ne sont pas par défaut contenues dans notre jeu de données. Il va donc falloir les définir dans notre tableau.
Nous allons créer pour cela une nouvelle colonne dans notre tableau que nous appellerons `class_age`.
Dans cette variable nous allons "découper" la variable `Age` selon les classes définies précédemment.
Pour cela nous allons utiliser la fonction cut qui crée des classes en excluant la valeur de la borne inférieure et en incluant celle de la borne supérieure.
Mathématiquement cela donne ceci : $)a;b]$
Avant cela on remarque que les âges dans notre fichier contiennent des virgules, il nous faut donc être plus précis sur la définition des bornes de nos classes. Ainsi nous allons définir à R nos 4 classes de la façon suivante:
1. )17.9-33.9] ans -> Nous indiquons 17.9 pour que les sujets de 18.0 ans soit inclus dans la classe et 33.9 pour que la classe suivante démarre à 34 ans
2. )33.9-54.0] ans -> Nous indiquons 33.9 pour que les sujets de 34.0 ans soit inclus dans la classe et 54.0 pour inclure les sujets de 54.0 ans dans la classe
3. )54.0-64.0] ans-> Nous indiquons 54.0 pour que les sujets de 55 ans soit inclus dans la classe (NB: il n'y a pas de valeur entre 54.1 et 54.9) et 64.0 pour inclure les sujets de 64.0 ans dans la classe
4. )64.0-89.9]ans -> Nous indiquons 64.0 pour que les sujets de 65 ans soit inclus dans la classe (NB: il n'y a pas de valeur entre 64.1 et 64.9) et 89.9 car ceci correspond à l'âge des sujets les plus agés
On vérifie visuellement si les différentes classes sont bien découpées. Pour nous faciliter la tâche nous allons d'abord réorganiser notre jeu de données en fontion de l'âge.
```{r}
data[order(data$Age),]
```
La classification s'est bien passée.
Néanmoins pour nous éviter des erreurs d'interprétation par la suite, nous allons renommer nos classes pour qu'elles correspondent aux noms de celles définies dans l'énoncé de la question 2.
```{r}
rename_class_age<-class_age
levels(rename_class_age)
levels(rename_class_age)[c(4)]<-">65"
levels(rename_class_age)
levels(rename_class_age)[c(3)]<-"55-64"
levels(rename_class_age)
levels(rename_class_age)[c(2)]<-"34-54"
levels(rename_class_age)
levels(rename_class_age)[c(1)]<-"18-34"
levels(rename_class_age)
data$class_age<-rename_class_age
```
On vérifie que le renommage s'est bien passé en regardant les niveaux des facteurs de la variable `class_age`
```{r}
levels(data$class_age)
```
Tout s'est bien passé
### 4.2.2 Tableau du statuts des femmes en fonction
### de l'habitude de tabagisme et de leur age.
Comme précedemment, on commence par représenter dans une table les effectifs pour chaque classe considérée
Pour avoir un résultat plus visualisable comme précédemment nous cherchons à représenter pour chaque habitude de tabagisme et classe d'âge le taux de mortalité.
Néanmoins nous ne pouvons pas simplement appliquer la fonction `prop.table`aussi simplement que précédement, c'est à dire juste sur la table d'effectifs créée. Ceci provient du fait que les calculs de proportion avec cette fonction ne peuvent tenir compte que d'un facteur à la fois. Il n'est donc pas possible d'avoir les proportions de mortalité par classe d'âge pour les fumeuses et les non fumeuses sans modifier un peu notre code.
Pour contourner ce problème nous allons, avant de calculer nos proportions, créer 2 sous-tableaux : l'un ne contenant que les données relatives aux fumeuses et l'autre au non fumeuses. ET nous allons calculer nos taux de mortalité pour chacune des habitudes de tabagisme à partir de ces tables
__Création des sous tableaux__
Et on peut aussi aisément inclure des figures. Par exemple:
```{r}
# Création des sous tableaux
fumeuse=data[data$Smoker=="Yes",]
non_fumeuse=data[data$Smoker=="No",]
summary(fumeuse)
summary(non_fumeuse)
```{r pressure, echo=FALSE}
plot(pressure)
```
Avec les `summary` de nos table `fumeuse`et `non_fumeuse` on peut vérifier que les effectifs comprennent bien uniquement des données sur le groupe en question. On retrouve bien pour la table `fumeuse` 582 observations et 732 dans la table `non fumeuse`; La somme des deux nous donnant bien notre nombre total d'observations : 1314.
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.
__Calculs des taux de mortalité__
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 les fumeuses_
Maintenant, à vous de jouer! Vous pouvez effacer toutes ces informations et les remplacer par votre document computationnel.
On commence par réeffectuer une table des effectifs du statut en fonction de l'âge pour les fumeuses.
```{r}
fum=table(fumeuse$Status, fumeuse$class_age)
fum
```
On note que l'on retrouve bien les mêmes effectifs que ceux indiqués dans la table `table_age`précedemment créée avec l'ensemble du jeu de données. Ceci nous confirme une nouvelle fois que la création de nos sous-tables s'est bien passé.
On calcule ensuite nos taux de mortalité
```{r}
prop_fum<-prop.table(fum,2)
prop_fum
```
On obtient que les taux de mortalité pour les classes d'âge considérées sont de:
Classe d'âge | Taux de mortalité (en %) arrondi à 0.1
On note que l'on retrouve bien les mêmes effectifs que ceux indiqués dans la table `table_age`précedemment créée avec l'ensemble du jeu de données. Ceci nous confirme une nouvelle fois que la création de nos sous-tables s'est bien passé.
Calcul des taux de mortalité
```{r}
prop_nonfum<-prop.table(nonfum,2)
prop_nonfum
```
On obtient que les taux de mortalité pour les classes d'âge considérées sont de:
Classe d'âge | Taux de mortalité (en %) arrondi à 0.1
Pour avoir une meilleure visualisation de nos données nous allons les représenter comme précédemment sur un barre plot
Pour pouvoir représenter uniquement les taux de mortalité en fonction des variables `classe_age` et de `Smoke` nous allons créer un nouveau data.frame contenant uniquement ces 3 informations
On commence par créer deux vecteurs répertoriant les taux de mortalité dans chacun des groupes.
```{r}
Vprop_fum <- as.vector(prop_fum[2,])
Vprop_nonfum <- as.vector(prop_nonfum[2,])
prop_morta<-c(Vprop_fum,Vprop_nonfum)
```
Puis on crée deux nouveaux vecteurs contenant pour chacun des 8 taux de mortalité répértoriés dans `prop_morta` la classe d'âge et d'habitude de tabagisme associée
Enfin on crée un data.frame que l'on appelle `morta`avec les 3 vecteurs que nous venons de créer
```{r}
morta<-data.frame(smoke,age,prop_morta)
```
On peut maintenant faire un barplot avec ces données.
Pour faire ce graphique nous allons faire appelle à la library `GrapheR`.
GrapheR est une interface utilisateur multiplateforme (Linux, Mac OS, Windows) permettant de réaliser des graphes hautement paramétrables sous R.(plus de détails: [Introduction à GrapheR](http://www2.uaem.mx/r-mirror/web/packages/GrapheR/vignettes/manual_fr.pdf))
```{r}
library(GrapheR)
run.GrapheR()
```
GrapheR crée des graphiques en fonction des choix fait dans son interface. Nous choisissons de créer un graphique en barres pour représenter le taux de mortalité en fonction de l'habitude de tabagisme. Les lignes de Code suivante correspondent au code sorti par `GrapheR` qu'il est possible de sauvegarder.
```{r}
# Loading of the dataset
dataset <- morta
attach(dataset)
# Preliminary data creation
means <- tapply(prop_morta,list(smoke,age),function(x) mean(x,na.rm=TRUE))
Ce résultat apparait suprenant car il semble contradictoire avec le précédent. Il donne l'impression que les non fumeuses sont moins impactées que les fumeuses. En effet à partir de la deuxième classe d'âge les taux de mortalité sont moins élevés chez les non fumeuses que chez les fumeuses. Or l'analyse faite à la question précédente indique que le taux de mortalité était plus élevé chez les non fumeuses.
### 4.2.5 Arrivez-vous à expliquer ce paradoxe ?
Ce résultat est en fait le fruit d'un paradoxe de Simpson.
Ce paradoxe peut être défini comme le fait qu’une corrélation peut disparaître ou même s’inverser suivant que l’on considère les données dans leur ensemble, ou bien segmentées par groupes. Ce paradoxe se produit lorsqu'il existe un __facteur de confusion__ et __lorsqu'un échantillon que l'on étudie n'est pas distribué de manière homogène__.Si ces deux conditions sont réunies alors à cause de la distribution hétérogène de l'échantillon, faire une analyse globale peut orienter un résultat qui peut être faux, ici le tabac est bon pour la santé car réduit le taux de mortalité. Cette tendance disparait lorsque l'on effectue une analyse en utilisant le facteur de confusion comme facteur de séparation des données.
Dans notre cas:
* Le facteur de confusion est l'âge des personnes qui joue sur le taux de mortalité
* La distribution non homogène provient du fait que dans l'échantillon non fumeuse on peut remarquer qu'il y a plus de personnes agée interrogées que dans la classe fumeuse
# 5.Question 3
## 5.1 Enoncé
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).
## 5.2 Réponse à la question 3
### 5.2.1 Conversion de la variable statut en binaire
__Conversion de la variable qualitative `Statut`en variable binaire: `Death`__
Nous allons commencer par recoder la variable `Statut` comme une variable binomiale que nous allons nommée `Death` que nous réutiliserons dans la suite de cet exercice.
Nous allons coder de la façon suivante les catérogies de la variable `Statut`:
Alive = 0
Dead = 1
`Death`sera donc une variable binaire (0 ou 1) et sera considérer comme une variable numérique
1. Transformation de la variable qualitative `Status` en variable binaire `Death`
On commence par créer un objet `Death` pour qu'il soit le reflet de la variable `Status` de notre jeu de donnée
```{r}
Death<-data$Status
```
On commence alors à modifier dans ce vecteur les noms des niveaux de cet objets
Ainsi, on remplace le niveau _Death_ par _1.0_ et _Alive_ par _0.0_
```{r}
levels(Death)
levels(Death)[c(2)]<-1.0
levels(Death)
levels(Death)[c(1)]<-0.0
levels(Death)
```
On vérifie la nature de `Death
```{r}
class(Death)
```
Avant de l'implémenter dans notre jeu de donnée, nous convestissons `Death` en une variable numérique
```{r}
Death<-as.integer(Death)
Death
summary(Death)
```
On constate ici un problème qui c'est produit alors de converssion:
R a remplacé chaque 0 par 1 et chaque 1 par 2
On corrige ceci pour avoir un vecteur correctement coder en 0=Alive et 1=Death
```{r}
Death2 <- numeric(length(Death))
for (i in 1:length(Death)) if (Death[i] == 2) Death2[i] <- 1 else Death2[i] <- 0
Death2
```
Le problème semble corrigé
On intègre donc ce vecteur `Death2` bien codé en 0 et 1 à notre jeu de donnée. Cette nouvelle variable de `data` est nommée `Death`
```{r}
data$Death<-Death2
View(data)
```
Vérification de la concordance entre les variables `Status` et `Death`
```{r}
table(data$Status,data$Death)
```
Toute les données Alive sont bien codées en 0 et les données Death en 1
On peut aussi vérifier que l'on retrouve les mêmes valeurs de taux de mortalité que ceux calculer dans la partie 3.2.2
```{r}
tableDeath<-table(data$Death,data$Smoker)
prop.table(tableDeath,2)
```
On retrouve bien un taux de mortalité de 31.4% pour les non fumeuse et de 23.9% pour les fumeuses.
__Notre conversion de la variable `Status`en binaire nommé `Death`semble donc correcte__
### 5.2.2 Etude du modèle `Death ~ Age`en fonction
### de l'habitude de tabagisme
Pour créer un modèle `Death ~ Age` par catégorie d'habitude de tabagisme, nous refaire des sous tableaux par catégorie d'habitude de tabagisme. Nous refaissons ces sous tableaux car les tables précédentes créées pour la question 2 n'intègre pas la variable `Death`
On commence donc par créer ces deux tables
```{r}
fumeuse=data[data$Smoker=="Yes",]
non_fumeuse=data[data$Smoker=="No",]
summary(fumeuse)
summary(non_fumeuse)
```
On peut donc ensuite créer nos régressions logistiques
__Analyse pour les fumeuses__
_Création de la régréssion logistique_
```{r}
reg_fumeuse = glm(data=fumeuse, Death ~ Age,
family=binomial(link='logit'))
summary(reg_fumeuse)
```
L'estimateur le plus probable du paramètre age dans le groupe fumeur est égale à __0.088977__ et l'erreur standard de cet estimateur de 0.008721
_Représentation graphique_
Nous allons utiliser le package `ggplot2`pour faire ce graphique. On charge aussi les packages `lattice` et `Matrix` nécéssaires aux calculs de la droite de régression
reg_non_fumeuse = glm(data=non_fumeuse, Death ~ Age,
family=binomial(link='logit'))
summary(reg_non_fumeuse)
```
L'estimateur le plus probable du paramètre age dans le groupe fumeur est égale à __0.107275__ et l'erreur standard de cet estimateut de 0.007806
_Représentation graphique_
Nous allons utiliser le package `ggplot2`pour faire ce graphique. On charge aussi les packages `lattice` et `Matrix` nécéssaires aux calculs de la droite de régression
### 5.2.3 Conclusion: Ces régressions vous permettent-elles
### de conclure sur la nocivité du tabagisme ?
Non ces regressions ne nous permettent pas de conclure sur la novicité du tabagisme dans ce contexte. En effet, les courbes des deux régressions sont relativement semblables.
Pour s'assurer de la non influence du tabac dans ce contexte nous nous proposons de construire un nouveau modèle de régression logistique qui tiendra compte de l'ensemble des paramètres à notre disposition pour évaluer quels paramètres mesurés peuvent influencer la mortalité.
Pour cela nous allons construire le modèle suivant:
`Death`~`Age`+ `Smoker` + `Age:Smoker`
Ce modèle tient donc compte de l'age, de l'habitude du tabasgisme et de l'interaction possible entre ces deux facteurs
```{r}
reg = glm(data=data, Death ~ Age*Smoker,
family=binomial(link='logit'))
```
Maintenant pour évaluer si nos paramètres ont une influence statistiquement significative sur notre variable `Death`, nous allons faire appel à la fonction `Anova` du package `car` qui permet un calcul selon la méthode Anova de type 2. Cette méthodologie de calcul permet notament de ne pas tenir compte de la potentielle confusion entre les variables.
On commence donc par charger le package `car`puis on applique la fonction Anova sur le modèle précédemment crée.
```{r}
library(car)
Anova(reg)
```
Ici on constate alors que seul l'âge à une p value inférieure à 5%. Ainsi ici le modèle le plus juste pour expliquer la variablité dans notre jeu de données semble être un modèle ne tenant compte que de l'âge des sujets.
Ceci nous conforte donc dans notre idée que __ce jeu de données ne nous permet de conclure sur les effets nocifs du tabac__
