---
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)
```

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)
library(knitr)
kable(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]

```{r,, fig.height=9, fig.width=14}
library(ggplot2)
library(ggthemes)
ggplot(tab_exo1, aes(x=Smoker, y=Taux_mortalite, group=Smoker, color=Smoker)) + 
geom_pointrange(aes(ymin=ICdown, ymax=ICup),lwd=1.8) + 
  scale_color_manual(values=c("blue","red"))
```

**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]="18-34"
data$classe_age[data$Age>=35 & data$Age<54]="35-54"
data$classe_age[data$Age>=55 & data$Age<64]="55-64"
data$classe_age[data$Age>=65 ]="65 et +"
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)
```

**Vous pourrez proposer une représentation graphique de ces données et calculer des intervalles de confiance si vous le souhaitez.**
```{r,, fig.height=9, fig.width=14}
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)
library(knitr)
kable(tab_exo2)
```
**En quoi ce résultat est-il surprenant?**
Dans ce cas là, on voit que l'on a une tendance à 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.Cependant les interval de confiance sont chevauchant 

```{r,, fig.height=9, fig.width=14}
tab_exo2$age_smoker=as.factor(paste(tab_exo2$Class_Age,tab_exo2$Smoker, sep="_"))
library(ggplot2)
library(ggthemes)
ggplot(tab_exo2, aes(x=age_smoker, y=Taux_mortalite, group=Smoker, color=Smoker)) + 
geom_pointrange(aes(ymin=ICdown, ymax=ICup),lwd=1.8) + 
  scale_color_manual(values=c("blue","red"))
```
**Arrivez-vous à expliquer ce paradoxe? De même, vous pourrez proposer une représentation graphique de ces données pour étayer vos explications.**

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)
```

### 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))
``` 

- 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)
```

**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 bien qu'il y ait une tendance à ce que les femmes dans la tranche d'âge autour de 40 ans qui sont fumeuses aient un taux de mortalité légèrement supérieur aux femmes non fumeuses

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à.