---
title: 'Exercice 3: Simson paradox'
author: "A.Changenet"
date: "02/04/2021"
output:
  pdf_document:
    pandoc_args: --listings
    toc: true
  html_document:
    toc: true
    theme: journal
documentclass: article
classoption: a4paper
header-includes:
- \hypersetup{colorlinks=true,pagebackref=true}
- \lstset{breaklines=true,basicstyle=\ttfamily}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


# Question 1

On commence par importer les données depuis le repertoire local: 
```{r}
data <- 
read.csv("~/Google Drive/Cours/mooc-rr/module3/exo3/module3_Practical_session_Subject6_smoking.csv")

```

Puis 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: 

```{r}
data.table <- table(data$Smoker, data$Status)
data.table

```

On calcule 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). L'argument margin permets de calculer les proportions relativement aux sommes des lignes (donc par statut de fumeuse ou non). Multiplier par 100 permets d'avoir un pourcentage. 
```{r}
data.table.prop <- prop.table(data.table, margin=1)*100
data.table.prop
```

Pour vérifier que notre calcul est bon, on ajoute la, porportion de vivant à la proportion de personne morte. On s'attends à avoir 100% dans les deux cas. 

```{r}
apply(data.table.prop,1,sum)
```
On a bien 100% pour les deux catégories.


Maintenant on fait un graphique de la proportion de mort en fonction du statut de smoker:

```{r}
dataf <- as.data.frame(data.table.prop)
Mort_rate <- dataf[dataf$Var2=="Dead",]
barplot(Mort_rate$Freq~Mort_rate$Var1,main = "Mortality rate VS smoking statuts",
        xlab = "Smoker",ylab = "proportion of dead")
```

On constate avec surprise que le taux de mortalité est plus elevé chez les non fumeuse que chez les fumeuse alors que l'on pourrait s'attendre à l'inverse. 


# Question 2

On commence par importer les données depuis le repertoire local: 
```{r}
data <- 
read.csv("~/Google Drive/Cours/mooc-rr/module3/exo3/module3_Practical_session_Subject6_smoking.csv")

```

Dans un premier temps, on s'assure que l'âge minimum n'est pas en dessous de 18 ans. 
```{r}
summary(data$Age)
```

C'est bien le cas (min = 18) donc on peut ajouter une variable liée à la classe d'âge. On considérera ici les classes suivantes : 18-34 ans, 34-54 ans, 55-64 ans, plus de 65 ans

```{r}
data$Ageclass <- cut(data$Age,breaks = c(18,34,54,64,Inf))
```


Puis 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}
data.table <- table(data$Smoker, data$Status,data$Ageclass)
data.table
```


On calcule dans chaque groupe (fumeuses / non fumeuses) de chaque classe d'âge, 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). L'argument margin permets de calculer les proportions relativement aux sommes des lignes (donc par statut de fumeuse ou non = 1), au sein de chaque catégorie d'âge (3). 

```{r}
data.table.prop <- prop.table(data.table, c(1,3))*100
data.table.prop
```


Maintenant on fait un graphique de la proportion de mort en fonction du statut de smoker:

```{r}
dataf <- as.data.frame(data.table.prop)
Mort_rate <- dataf[dataf$Var2=="Dead",]
barplot(Mort_rate$Freq~Mort_rate$Var1+Mort_rate$Var3,beside=T,
        main = "Mortality rate VS smoking statuts in different age classes",
        xlab = "Age",ylab = "proportion of dead",legend=T,args.legend=list(title="Smoker"))

```

On constate que cette fois, la proportion de mort est plus élevé chez les fumeurs quelque soit la tranche d'âge. 
C'est paradoxal avec le résultat précédent, ou on observait un taux de mortalité plus faible chez les fumeurs. C'est ce que l'on appelle le [paradoxe de Simpson](https://fr.wikipedia.org/wiki/Paradoxe_de_Simpson). 

Deux explications possibles: 
* L'âge que nous n'avons pas pris en compte dans la première analyse est un paramètre important (et à priori il est important car on s'attends globalement a ce que les plus jeune meurt moins (fumeurs ou pas), et que les plus vieux meurent plus (fumeurs ou pas). Cet effet de sénescence est pris en compte dans le second temps. 

Seconde explication: 
* Le déséquilibre dans les effectifs (globaux, et au sein de chaque classes) pourrait être responsable de ces conclusions trompeuses. 


## Visualisation graphique: 


```{r}
summary(as.factor(data$Smoker))
summary(as.factor(data$Status))
table(data$Smoker, data$Status)

```
* On remarque dans un premier temps ci dessus, le déséquilibre au sein des variables qualitatives "Smoker" et "Status". En revanche: 

```{r}
summary(as.factor(data$Ageclass))

```
Les classes d'âge ont l'air homogènes. Poursuivons en regardant les effectifs en fonction de plusieurs variables: 

```{r}
data.table <- table(data$Smoker, data$Status,data$Ageclass)
data.table
```
Quand on regarde les statuts et conditions en fonction de ces tranches d'âge en revanche, on remarque un total déséquilibre. 

En particulier, dans la tranche d'âge la plus élevée, les fumeurs sont compltement sous représentés en nombre par rapport aux non fumeurs. 
Comme c'est dans cette tranche dâge que l'on trouve le taux de mortalité le plus élevé (sénéscence) quelque soit le statut de fumeur ou non, c'est cette tranche d'âge qui biaise totalement nos résultats puisque même si la proportion de mort chez les fumeurs est aussi élevé que chez les non fumeurs, le **nombre** de mort lui est bien plus faible (42 contre 165), impactant fortement nos conclusions générale sur la dynamique globale qui ne reposent que sur les proportions totale !
Il est donc important de prendre en compte la variable âge. 


* Pour s'en rendre compte on peut aussi visualiser les données: 

```{r,fig.height=7}
dataf <- as.data.frame(data.table)
Mort_rate <- dataf[dataf$Var2=="Dead",]
Live_rate <- dataf[dataf$Var2=="Alive",]

par(mfrow=c(1,2))
barplot(Mort_rate$Freq~Mort_rate$Var1+Mort_rate$Var3,beside=T,
        main = "Mortality count\n VS smoking statuts\n in different age classes",
        xlab = "Age",ylab = "Number of dead",legend=T,
        args.legend=list(x="topleft",title="Smoker"))
barplot(Live_rate$Freq~Mort_rate$Var1+Mort_rate$Var3,beside=T,
        main = "Living count\n VS smoking statuts\n in different age classes",
        xlab = "Age",ylab = "Number of living",legend=T,
        args.legend=list(title="Smoker","topright"))
```
On voit ici le déséquilibre flagrant entre tranche d'âge pour les vivants comme pour les morts. Dans la tranche d'âge la plus élevée on voit clairement que les non fumeurs morts sont surreprésentés en nombre (alors qu'ils représentent une proportion semblables chez les fumeurs).


# Question 3: modèle binomial:

Dans un premier temps on recode la variable "Status" en "Death" avec 0 indiquant la survie, et 1 un évènement de mortalité
```{r}
data$Death <- "0"
data[data$Status=="Dead","Death"] <- 1
data$Death <- as.numeric(data$Death)
```


Dans un second temps, 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. 

Le plus simple est d'inclure un effet âge et un effet "fumeur" dans notre modèle: 

```{r}
M1 <- glm(Death ~ Smoker+Age, data = data, family = binomial(logit))
```


On regarde les sorties de notre modèle:
```{r}
summary(M1)
```
On constate que l'effet $age est significatif mais que l'effet "smoker" ne l'est pas. Voyons de plus prêt à l'aide des packages 'cowplot' et 'ggeffects' que l'on pourra installer grace aux commandes: `install.packages(c("cowplot","ggeffects"))`


```{r}
cowplot::plot_grid(plotlist = plot(ggeffects::ggeffect(M1)))
```
On remarque que les fumeurs ont une probbailité plus grande de mort que les non fumeurs. On remarque aussi que la probabilité de mort augmente drastiquement avec l'âge (c'est la sénéscence).


On peut égalament regarder la même courbe pour les deux groupes.
Pour se faire, on définit une plage d'âge pour laquelle on veut des prédictions, disons 18 à 100 ans, pour un groupe de fumeur et un groupe de non fumeurs. On transforme ces vecteurs en un data frame pour pouvoir faire nos prédictions: 

```{r}
Age <- seq(18,100,1)
SmokerYes <- rep("Yes",length(Age))
SmokerNo <- rep("No",length(Age))
newdata <- as.data.frame(cbind(Age,SmokerYes,SmokerNo))

```

Ici on vérifie qu'on a bien la structure attendue: 
```{r}
newdata
```

On prédit successivement les réponses (en terme de mortalité), pour le groupe des fumeurs de 18 à 100 ans et celui des non fumeur de 18 à 100 ans. On précise que l'on veut aussi récupérer l'intervalle de confiance grace a l'argument `se.fit = T` et on récupère le tout dans deux nouveaaux objects.
```{r}
Pred.Yes <- predict(M1,list(Smoker=SmokerYes,Age=Age),se.fit = T,type="response")
Pred.No <- predict(M1,list(Smoker=SmokerNo,Age=Age),se.fit = T,type="response")
```


On représente nos points observés d'abord. Puis on ajoute successivement le groupe des fumeurs avec l'intervalle de confiance, et celui des non fumeurs, avec la légende adequat.
La commande `lines`permet d'ajouter la ligne que l'on souhaite (ici la réponse en terme de mortalité des fumeurs ou des non fumeurs). Dans un second temps on ajoute les intervalles de confiance grace à la commande `polygon`. 

```{r}
plot(data$Death~data$Age,pch = 16, xlab = "Age", ylab = "Death probability")
polygon(x=c(newdata$Age,rev(newdata$Age)),
        y=c(Pred.Yes$fit-Pred.Yes$se.fit,
            rev(Pred.Yes$fit+Pred.Yes$se.fit)),
        col=adjustcolor("blue",alpha.f=0.1),
        border=NA)
lines(Pred.Yes$fit~newdata$Age,col="blue")

polygon(x=c(newdata$Age,rev(newdata$Age)),
        y=c(Pred.No$fit-Pred.No$se.fit,
            rev(Pred.No$fit+Pred.No$se.fit)),
        col=adjustcolor("red",alpha.f=0.1),
        border=NA)
lines(Pred.No$fit~newdata$Age,col="red")

legend("left",title="Smoker",c("Yes","No"),col=c("blue","red"),fill=c("blue","red"))
```

On constate bien qu'on prédit, à âge égale, moins de mortalité chez les non fumeurs que chez les fumeurs, quelquesoit l'âge. Le tabac semble donc nocif, au regard de ces données.


Toutes les analyses ont été réalisé sur R 
```{r}
sessionInfo()
```