---
title: "Exercice paradoxe de Simpson"
author: "Hugues de Courson"
date: "La date du jour"
output:
  pdf_document:
    df_print: paged
---


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

### Début: chargement des packages pouvant être utiles

```{r package,message=F, warning=F}
library(epiDisplay)
library(epiR)
library(prettyR)
library(knitr)
library(kableExtra)
```


### Puis l'import des données

Dans un premier temps, il s'agissait de télécharger le fichier dans mon dossier local, ensuite on commit cette action et avec un push il est envoyé sur gitlab.

Ensuite il s'agit de le charger dans R :

```{r data import}
data <- read.csv("Subject6_smoking.csv")
View(data)
```

Vérification du format des données :

```{r data type}
str(data)
```

Tout semble ok.

### Mission 1

## Creation du tableau :

L'objectif est de créer un tableau avec le nombre total de femme vivantes ou décédées sur la période en fonction des habitudes de tabagisme. 

```{r tableau 1}
vivantes_fum <- sum(data$Status[data$Smoker=="Yes"]=="Alive")
vivantes_nonfum <- sum(data$Status[data$Smoker=="No"]=="Alive")

decedees_fum <- sum(data$Status[data$Smoker=="Yes"]=="Dead")
decedees_nonfum  <- sum(data$Status[data$Smoker=="No"]=="Dead")
tab_1 <- data.frame(Vivantes = c(vivantes_fum,vivantes_nonfum), Decedees = c(decedees_fum,decedees_nonfum), row.names = c("Fumeuses","Non fumeuses"))
kable(tab_1,align = 'l')
```


## Calcul des taux de mortalité :
# Creation d'un tableau de contigence
```{r taux mortalite}

mort_fum <- sum(data$Status[data$Smoker=="Yes"]=="Dead")
mort_nonfum <- sum(data$Status[data$Smoker=="No"]=="Dead")
vivant_fum <-sum(data$Status[data$Smoker=="Yes"]=="Alive")
vivant_nonfum <- sum(data$Status[data$Smoker=="No"]=="Alive")
nb_fum <- sum(data$Smoker=="Yes")
nb_nonfum <- sum(data$Smoker=="No")
nb_viv <- sum(data$Status=="Alive")
nb_deces <-sum(data$Status=="Dead")
tot <- sum(data$Status=="Alive"|data$Status=="Dead")

tab_2 <- data.frame(cbind(rbind(mort_fum, vivant_fum, nb_fum)), rbind(mort_nonfum,vivant_nonfum, nb_nonfum),rbind(nb_deces,nb_viv,tot))
row.names(tab_2) <- c("Decedees", "Vivantes","Total")
names(tab_2) <- c("Fumeuses", "Non Fumeuses", "Total")

kable(tab_2,align = "c")
```

# Representation graphique
Ceci passe par la création d'un tableau avec les pourcentages de mortalité
```{r representation graphique taux de mortalite}
tab_3 <- data.frame(cbind(rbind((sum(data$Status[data$Smoker=="Yes"]=="Dead")/sum(data$Smoker=="Yes"))*100,(sum(data$Status[data$Smoker=="Yes"]=="Alive")/sum(data$Smoker=="Yes"))*100)), 
                    rbind((sum(data$Status[data$Smoker=="No"]=="Dead")/sum(data$Smoker=="No"))*100,(sum(data$Status[data$Smoker=="No"]=="Alive")/sum(data$Smoker=="No"))*100))
                                
names(tab_3) <- c("Fumeuses", "Non fumeuses")
row.names(tab_3) <- c("Decedees", "Vivantes")

kable(tab_3, align = "c")

barplot(as.matrix(tab_3), col = c("black","gray"), space = 1.5,width = 0.5) 
legend("center",xpd=NA, legend = c("Decedees", "Vivantes"), fill = c("black","gray"))
```

# Calcul des intervalles de confiance des proportions
```{r intervalle de confiance}

#  creation de la fonction :

IC = function(x) {
  y <- x/100
  inf = y-(1.96*sqrt((y*(1-y))/n))
  sup = y+(1.96*sqrt((y*(1-y))/n))
  print(c(x,inf*100,sup*100))
}

```
Chez les fumeuses : 
```{r IC fumeuses}

x <- tab_3[1,1]
n <- as.numeric(tab_2[3,1])

IC(x)
```
Chez les non fumeuses :
```{r IC non fumeuses}
x <- tab_3[1,2]
n <- tab_2[3,2]

IC(x)

```

Au vu de l'ensemble des résultats, on fait le constat suivant :
  Il semblerait que le taux de mortalité est plus important chez les non-fumeuses...
  Ceci va à l'encontre des connaissances sur le tabac, on pourrait s'attendre à ce que les non fumeuses décèdent moins. 
  
### Mission 2

Il s'agit de reprendre les données, cette fois-ci en fonction de la classe d'âge

## Première étape : création de la variable classe d'âge

```{r var classe dage}

data$cl_age <- as.factor(ifelse(data$Age>=18&data$Age<35,"[18;34]",
                                ifelse(data$Age>=35&data$Age<54,"[35;54]",
                                       ifelse(data$Age>=55&data$Age<65,"[55;64]","[65;["))))

tab1(data$cl_age)

```

## Deuxième étape on refait les tableaux de contingence mais cette fois-ci de manière plus optimisée

Pour cela on utilise [ce lien](http://olivier.godechot.free.fr/hoparticle.php?id_art=465) pour connaître la méthode. 
Et [ce lien](https://haozhu233.github.io/kableExtra/awesome_table_in_pdf.pdf) pour la réalisation de beaux tableaux.

```{r table contingence}

tab_5 <- ftable(data[,c(1,2,4)])
tab_5 <- round(prop.table(tab_5,2)*100,1)
tab_5 <- rbind(tab_5, tab1(data$cl_age)$output.table[,1])
tab_5 <- as.table(tab_5)
row.names(tab_5) <- c("Non fumeuses vivantes ", "Non fumeuses decedees", "Fumeuses vivantes", "Fumeuses decedees" ,"Effectifs")
colnames(tab_5) <- c("[18;34]", "[35;54]", "[55;64]", "[65;[")

kable(tab_5,"latex", align = "c") %>% kable_styling(latex_options = "striped", stripe_index = c(1,2))
```

## Et maintenant on trace les barplots

Super [lien](http://zoonek2.free.fr/UNIX/48_R/04.html) pour les paramètres graphiques.
Et autre [lien](https://sites.google.com/site/rgraphiques/realiser-des-graphiques-avec-le-logiciel-r/diagrammes-en-barres) pour les diagrammes en barre.

```{r deuxieme barplot}

par(mar = c(5,5,5,13))
barplot(tab_5[c(1:4),], col = c("gray","black","gray","black"), density = c(30,30,40,100),angle = c(70,70,0,0))
xmin <- par("usr")[1]
xmax <- par("usr")[2]
ymin <- par("usr")[3]
ymax <- par("usr")[4]
par(xpd=TRUE)
lambda <- 0.025
legend(((1 + lambda) * par("usr")[2] - lambda * par("usr")[1]),50, legend = c("Non fumeuses vivantes", "Non fumeuses decedees", "Fumeuses vivantes", "Fumeuses decedees"), fill = c("gray","black","gray","black"),density = c(30,30,40,100),angle = (c(70,70,0,0)))


```

**NB** : dans le code précédent :

```{r explication par mar}
par(mar = c(5,5,5,15))
```
On définit les paramètres de marge, dans l'ordre c(bas,gauche,haut,droite)

```{r explication lsite usr}
xmin <- par("usr")[1]
xmax <- par("usr")[2]
ymin <- par("usr")[3]
ymax <- par("usr")[4]
```
On récupère les valeurs min et max du plot précédent


## Et maintenant on peut obtenir les intervalles de confiance
<!--# Pour le premier groupe d'age non fumeuses decedees-->
```{r IC premier groupe dage non fum, include=FALSE}

x <- tab_5[2,1]
n <- tab_5[5,1]

IC1<- round(IC(x),1)
IC1 <- paste(c(IC1[1],"[",IC1[2],"-",IC1[3],"]"),collapse = "")

```
<!--# Pour le premier groupe d'age fumeuses decedees-->
```{r IC premier groupe dage fum, include=FALSE}

x <- tab_5[4,1]
n <- tab_5[5,1]

IC2 <- round(IC(x),1)
IC2 <- paste(c(IC2[1],"[",IC2[2],"-",IC2[3],"]"), collapse = "")
```

<!--# Pour le deuxième groupe d'age non fum-->
```{r IC deuxieme groupe dage non fum, include=FALSE}

x <- tab_5[2,2]
n <- tab_5[5,2]

IC3 <- round(IC(x),1)
IC3 <- paste(c(IC3[1],"[",IC3[2],"-",IC3[3],"]"),collapse = "")

```
<!--# Pour le deuxième groupe d'age fum-->
```{r IC deuxieme groupe dage fum, include=FALSE}

x <- tab_5[4,2]
n <- tab_5[5,2]

IC4 <- round(IC(x),1)
IC4 <- paste(c(IC4[1],"[",IC4[2],"-",IC4[3],"]"),collapse = "")

```
<!--# Pour le troisieme groupe d'age non fum-->
```{r IC troisieme groupe dage non fum, include=FALSE}

x <- tab_5[2,3]
n <- tab_5[5,3]

IC5 <- round(IC(x),1)
IC5 <- paste(c(IC5[1],"[",IC5[2],"-",IC5[3],"]"),collapse = "")

```
<!--# Pour le troisieme groupe d'age  fum-->
```{r IC troisieme groupe dage  fum, include=FALSE}

x <- tab_5[4,3]
n <- tab_5[5,3]

IC6 <- round(IC(x),1)
IC6 <- paste(c(IC6[1],"[",IC6[2],"-",IC6[3],"]"),collapse = "")
```
<!--# Pour le quatrieme groupe d'age non fum-->
```{r IC quatrieme groupe dage non fum, include=FALSE}

x <- tab_5[2,4]
n <- tab_5[5,4]

IC7<- round(IC(x),1)
IC7 <- paste(c(IC7[1],"[",IC7[2],"-",IC7[3],"]"),collapse = "")

```

<!--# Pour le quatrieme groupe d'age fum-->
```{r IC quatrieme groupe dage fum, include=FALSE}

x <- tab_5[4,4]
n <- tab_5[5,4]

IC8<- round(IC(x),1)
IC8 <- paste(c(IC8[1],"[",IC8[2],"-",IC8[3],"]"), collapse = "")

```

Sous forme d'un tableau :
```{r, message=FALSE}

tab_6 <- cbind(rbind(IC1,IC2), rbind(IC3,IC4), rbind(IC5,IC6), rbind(IC7,IC8))
colnames(tab_6) <- colnames(tab_5)
row.names(tab_6) <- c("Non Fumeuses", "Fumeuses")
kable(tab_6,align = "c")

```

Une fois de plus à la vue de ces résultats on est étonnés ! En effet, autant c'est cohérent chez les non fumeurs mais chez les fumeurs il y a une discordance car ce ne sont pas les personnes les plus agées qui décèdent le plus. Il y a probablement une interaction entre l'âge et le fait d'être fumeur, où alors les catégories sont mal faites...

### Mission 3
Il s'agit maintenant de faire une régression logistique pour utiliser l'âge en continu et ainsi ne pas perdre d'information.

## On commence par bien coder la variable décès

```{r recodage deces}

data$deces <- as.factor(ifelse(data$Status=="Alive",0,1))

```

## Ensuite on fait la regression logistique, d'abord univariable.

#Pour le tabac :

```{r univar tabac}

regunita <- glm(deces~Smoker, family = binomial, data = data)

# pour obtenir les OR

logistic.display(regunita)

```

#Pour l'âge :

```{r uni age}

reguniage <- glm(deces~Age, family = binomial(), data = data)

logistic.display(reguniage)

```

La subtilité pour l'âge est qu'il faut vérifier l'hypothèse de linéarité du Logit car variable quantitative et hypothèse sous-jacente d'un modèle de régression logistique.
Pour cela on introduit un polynome fractionnaire dans le modèle de régression.

```{r linearite logit}

library(mfp)

reguniage2 <- mfp(deces~fp(Age, df =4, select = 1, scale = T), family = binomial, data = data)

summary(reguniage2)
```

En lisant le summary, le meilleur polynome est de degrès 1 donc l'hypothèse de linéarité est respectée.

## Modèle multivariable

```{r modele multi}

regtot <- glm(deces~Age+Smoker, family = binomial, data = data)

logistic.display(regtot)
```


Et là on constate une inversion de l'odds ratio associé au tabac lorsque l'on ajuste sur l'âge...