---
title: "Autour du paradoxe de Simpson"
author: "Arnaud Legrand"
date: "3/12/2020"
output: pdf_document
---

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

# Récupération des données
```{r}
data_url = "https://gitlab.inria.fr/learninglab/mooc-rr/mooc-rr-ressources/-/raw/master/module3/Practical_session/Subject6_smoking.csv"
data_filename = "smoking.csv";

if (!file.exists(data_filename)) {
    download.file(data_url, data_filename)
}
df = read.csv(data_filename);
```

Vérifions si tout à l'air en ordre. A priori, pas de valeurs manquantes, la lecture et le codage des données ont l'air s'être bien passé. On a pour chaque personne, l'indication de si elle fume ou pas, de si elle est décédé durant la période et son age au début de la période. Les visualisations futures me permettront peut-être de repérer des valeurs étranges mais en attendant, ça ira bien.
```{r}
summary(df)
str(df)
head(df)
tail(df)
```

# Calcul des effectifs et de la mortalité globale
Pour celà, nous utiliserons les paquets du tidyverse qui sont bien pratiques et très souples. Je vais aussi fixer une palette de couleur pour les différents niveaux des facteurs.
```{r}
# I need Hmisc this to easily compute binary confidence intervals. 
# It has to be loaded before dplyr otherwise this results in some name conflicts.
library(Hmisc) 
library(dplyr)
library(ggplot2)
mystyle = list(theme_classic(),scale_fill_manual(
    values = c("Alive"="#377EB8", "Dead"= "#E41A1C", "No"="#4DAF4A", "Yes"= "#E41A1C")))
```
Commençons, comme demandé, par calculer les effectifs et la mortalité selon le tabagisme des individus.
```{r}
df %>% group_by(Smoker,Status) %>% summarize(Number=n()) -> df_overview
df_overview %>% tidyr::spread(Status,Number) %>% 
    mutate(Mortality = 100*Dead/(Alive+Dead)) -> df_mortality_overview
df_mortality_overview
```

Tiens, oui, surprenant en effet que la mortalité soit supérieure chez les non fumeurs. Ça va un peu à l'encontre de ce à quoi on pourrait s'attendre. Calculons la confiance naivement.
```{r}
df_mortality_overview %>% group_by(Smoker) %>%
    do(data.frame(., Conf = Hmisc::binconf(.$Dead, .$Alive+.$Dead, 
                                           alpha = 0.05, method = "wilson"))) ->
  df_mortality_overview
df_mortality_overview
```

Une petite représentation graphique avec uniquement la mortalité.
```{r}
ggplot(df_mortality_overview, aes(x=Smoker,y=Conf.PointEst*100)) + geom_bar(aes(fill=Smoker), stat = "identity") + 
    geom_errorbar(aes(ymin=Conf.Lower*100,ymax=Conf.Upper*100),width=.2) +
    ylab("Mortality (%)") + mystyle
```

Une petite représentation graphique avec les effectifs:
```{r}
ggplot(df_mortality_overview %>% select(Smoker, Alive, Dead) %>% tidyr::gather(key=Status,value=Number,-Smoker),
       aes(x=Status,y=Number)) + geom_bar(aes(fill=Status), stat = "identity", position=position_dodge()) + 
    geom_text(data=df_mortality_overview, aes(x="Dead",y=Dead+50, label=paste0("Mortality:\n",round(Mortality,1),"%"))) +
    facet_wrap(~Smoker, labeller = label_both)  + 
    mystyle
```

# Utilisation des catégories d'âges
On nous propose de définir des catégories d'âges. Allons-y:
```{r}
df %>% mutate(Age_Cat = case_when(
  Age >=18 & Age <34 ~ "18-34",
  Age >=34 & Age <55 ~ "34-54",
  Age >=55 & Age <65 ~ "55-64",
  Age >=65 ~ "65+")) %>% 
  mutate(Age_Cat = recode(as.factor(Age_Cat), 
                "0" = "18-34", "1" = "34-54", 
                "2" = "55-64", "3" = "65+")) -> df
```

Puis recalculons les effectifs et la mortalité  comme précédemment:
```{r}
df %>% group_by(Smoker,Status,Age_Cat) %>% summarize(Number=n()) %>%
    tidyr::spread(Status,Number) %>% mutate(Mortality = 100*Dead/(Alive+Dead)) -> df_grouped
df_grouped %>% arrange(Age_Cat)
```

Une petite représentation graphique de la mortalité. Effectivement il est paradoxal que la mortalité soit globalement supérieure pour les non fumeurs mais que pour **chaque** catégorie d'âge, elle soit inférieure!
```{r}
    ggplot(df_grouped, aes(x=Age_Cat,y=Mortality,fill=Smoker)) +
        geom_bar(stat="identity", position=position_dodge()) + 
        ylab("Mortality (%)") + ylim(0,100) + mystyle

```

Faisons une autre visualisation, plus explicite sur les effectifs et permettant de commencer à comprendre d'où peut venir le problème. En gros, pas de différence de mortalité significative pour les très jeunes et les très vieux mais une différence très importante pour les gens d'âge moyen. Par contre très peu de vieux fumeurs (de vieilles fumeuses en l'occurence). Fumer tue effectivement (on meurt plus tôt) mais les effectifs trompent les ratios.
```{r}
ggplot(df %>% group_by(Smoker,Status,Age_Cat) %>% summarize(Number=n()), 
    aes(x=Age_Cat,y=Number)) +
    geom_bar(aes(fill=Status), stat = "identity",position=position_dodge()) +     
    geom_text(data=df_grouped, aes(x=Age_Cat,y=Dead+10, 
                                   label=paste0("       ",round(Mortality,1),"%"))) +
    facet_wrap(~Smoker, labeller = label_both) + 
    mystyle
```


# Régression logistique
Tentons une régression logistique qui ne devrait pas être impactée par ces catégorisations arbitraires d'âge.
```{r}
summary(glm(data=df, Status ~ Age * Smoker, family=binomial(link='logit')))
```

Ah, on voit un petit effet mais ce n'est pas vraiment significatif (pour Smoker... évidemment que Age est significatif, plus on vieillit plus on a de chances de mourrir!). 

Essayons à tout hasard une régression pour chaque catégorie on pourra voir où les courbes de régression s'intersectent en prenant en compte la confiance
```{r}
summary(glm(data=df %>% filter(Smoker=="Yes"), 
            Status ~ Age, family=binomial(link='logit')))
summary(glm(data=df %>% filter(Smoker=="No"), 
            Status ~ Age, family=binomial(link='logit')))
```

Bon, pas facile à lire mais graphiquement, ça devrait être plus explicite. Pour éviter l'*overplotting*, je met un peu de jitter vertical sur mes points (manuellement plutôt qu'avec geom_jitter afin de séparer les *Smoker*).
```{r}
df %>% mutate(Status_num=ifelse(Status=="Dead",1,0),
              y = Status_num + (as.numeric(as.factor(Smoker))-2)*0.05 +
                  0.05*runif(n())) -> df_raw
ggplot(df_raw, aes(x=Age, y=Status_num, color=Smoker)) + 
    geom_point(alpha=.4, aes(y=y)) +
    geom_smooth(method="glm", method.args = list(family = "binomial"),fullrange = TRUE) +
    theme_bw() +  xlim(15,95) + ylab("Death occurence / probability") +
    scale_color_manual(values = c("No"="#4DAF4A", "Yes"= "#E41A1C"))
```

On voit nettement le décalage entre les deux courbes même si les régions de confiances s'interceptent. Avec des effectifs plus importants, la séparation serait nette. Le point important n'est pas la mortalité après 20 ans mais l'âge auquel on meurt...