DATA_MINING.Rmd

---
title: "DEVOIR_S2 DATA_MINING"
author: "SABAYE Fried-Junior, Caleb KASHALA, Mohamad EL KAWASS, Hicham AZOUD"
header-includes:
- \usepackage[Glenn]{fncychap}
- \usepackage{fancyhdr}
fontsize: 14pt
lang: 'fr'
geometry: a4paper,top=2cm,bottom=2cm,left=2cm,right=2cm
output: 
   pdf_document : 
     highlight: tango
     keep_tex: yes
     number_section: yes
     toc: yes
---
```{r, echo=FALSE}
knitr::opts_chunk$set(echo = FALSE, warning=FALSE, message=FALSE,fig.pos = "!h", sanitize=TRUE, fig.align='center',strip.white = TRUE)
```



```{r}
library(ISLR)
library(missMDA)
library(MASS)
library(ROCR)
library(dplyr)# en particulier pour la syntae %>%
library(rsample)
library(class)# package contenant la focntion knn
library(rpart) # package pour les arbres CART
library(rpart.plot)
library(randomForest)
library(e1071)
library(tidyr)
library(ggplot2)
library(parallel)
library(doParallel)
library(caret)
library(parallel)
library(FactoMineR)
library(mice)
library(plotROC)
library(ROCR)
library(ada)
library(precrec)
library(gridExtra)
library(formattable)
library(stargazer)
library(kableExtra)
set.seed(124)
```

\newpage 

# Contexte et problématique 

Cette étude fait suite à l'analyse factorielle éffectuée au **semestre précédent**. 

La base de données renseigne sur diverses variables et paramètres concernant des accidents de vélo. 

Les dix premières lignes et colonnes de la base de données figurent dans le tableau suivant : 

```{r}
data<-read.table("bike_crash.csv", header = T, sep = ";")
newdataqual<-data[,-c(1,2,3,4,12,13,14,15,16,17,18,19,20,21,22,23,24,25,27,28,30,34,38,39,40,41,44,45,46,47,50,51,54)]
newdataqual$Bike_Injur<-as.character(newdataqual$Bike_Injur)
newdataqual$Bike_Injur[newdataqual$Bike_Injur=="A: Disabling Injury"]<-"Oui"
newdataqual$Bike_Injur[newdataqual$Bike_Injur=="B: Evident Injury"]<-"Oui"
newdataqual$Bike_Injur[newdataqual$Bike_Injur=="C: Possible Injury"]<-"Non"
newdataqual$Bike_Injur[newdataqual$Bike_Injur=="Injury"]<-"Oui"
newdataqual$Bike_Injur[newdataqual$Bike_Injur=="K: Killed"]<-"Oui"
newdataqual$Bike_Injur[newdataqual$Bike_Injur=="O: No Injury"]<-"Non"
df<-newdataqual[complete.cases(newdataqual), ]
```

```{r}
kable(df[1:10,1:10],align = "c") %>%
kable_styling(latex_options =c("striped","hold_position","scale_down"),full_width = F) 
```




Le but de cette étude est de prédire la variable : *Bike_Injur*. Pour cela nous allons utiliser diverses méthodes d'apprentissage supervisé. 

```{r}
d<-data.frame(table(newdataqual$Bike_Injur))
colnames(d)<-c("Bike_Injur","Nombre")

kable(d,align = "c",caption="Effectif de la variable à prédire") %>%
kable_styling(latex_options =c("striped","hold_position"),full_width = F) 
```


&nbsp;

&nbsp;

&nbsp;


## Méthodologie  

Nous allons séparer la base de données en deux échantillons : Apprentissage et Test. 

La base de données sera séparée comme suit :  2/3 des données constitueront l'*échantillon d'Apprentissage* qui sera utilisé pour construire les modèles ; les 1/3 restant qui constitueront l'*échantillon Test* seront utilisés pour tester nos modèles.

Pour chacunes des méthodes d'apprentissage supervisé qui sera utilisée : on commencera par construire le(s) modèle(s) en fonction des différents paramètres propres à chacun d'eux et ce, sur l'échantillon d'apprentissage. Ensuite, dans le but de déterminer le meilleur paramétrage de chaque méthode : nous appliquerons ce(s) modèle(s) ainsi construit(s) sur l'échantillon test, puis nous comparerons les diverses mesures de performance. Ce qui nous permettra de choisir, si il le faut, le meilleur modèle de chaque méthode. 
Enfin, nous comparerons les résultats des meilleurs modèles de chaque méthode afin de déterminer celui permettant de mieux prédire.

```{r}
data_split <- df %>% initial_split(prop = 2/3)
test_data <- data_split %>% testing()
train_data <- data_split %>% training()
df2<-train_data
```

\newpage 

**Quelques précisions** 

$\textbf{1)}$ Après la transformation, la sélection des variables et l'affectation des effectifs réalisée $\textbf{au premier semestre}$, les variables restantes pouvant être utilisées  sont : *Bike_Age,  Bike_Alc_D, Bike_Dir, Bike_Injur, Bike_Pos, Bike_Race, Bike_Sex, Crash_Hour, Crash_Loc, Crash_Time, Crash_Type, Crash_Ty_1, Developmen, DrvrAge_Gr, Drvr_Age, Drvr_Alc_D, Drvr_Injur, Drvr_Race, Drvr_Sex, Hit_Run, Light_Cond, Num_Lanes, Num_Units,  Rd_Charact, Rd_Class, Rd_Conditi, Region, Rural_Urba et Workzone_I.*

&nbsp;


$\textbf{2)}$ En apprentissage automatique supervisé, une matrice de confusion est une matrice qui mesure la qualité d'un système de classification. Elles s'interprétent comme suit : 

```{r}
ligne<-c("Vrais négatifs","Faux négatifs")
lignee<-c("Faux positifs","Vrais positifs")
tabbb<-rbind(ligne,lignee)
rownames(tabbb)<-c("Non","Oui")
colnames(tabbb)<-c("Non","Oui")

kable(tabbb,align = "c",caption="Effectif de la variable à prédire") %>%
kable_styling(latex_options =c("striped","hold_position"),full_width = F)%>%
add_header_above(c("Classe réelle ", "Classe estimée par le classificateur" = 2))


```
Un vrai positif(VP) est un résultat où le modèle prédit correctement la classe positive. De façon analogue, un vrai négatif(VN) est un résultat où le modèle prédit correctement la classe négative.

Un faux positif(FP) est un résultat où le modèle prédit incorrectement la classe positive et un faux négatif(FN) est un résultat où le modèle prédit incorrectement la classe négative.

&nbsp;

$\textbf{3)}$ À partir de la matrice de confusion on peut dériver tout un tas de critères de performance. Parmi lesquels :

$\bullet$ **La sensibilité** : le taux de vrais positifs, c’est à dire la proportion de bléssés que l’on a correctement identifiés. C’est la capacité de notre modèle à détecter toutes les blessures.

$$Sensibilité=\frac{VP}{VP+FN}$$

$\bullet$ **La spécificité** : le taux de vrais négatifs, autrement dit la capacité à détecter toutes les  situations où il n’y a pas de bléssures. C’est une mesure complémentaire de la sensibilité. 

$$Spécificité=\frac{VN}{FP+VN}$$

$\bullet$  **La précision** : la proportion de prédictions correctes parmi les blessures que l’on a prédites positives. C’est la capacité de notre modèle à prédire qu'un individu à été bléssé que si il à été réellement bléssé.

$$Précision=\frac{VP}{VP+FP}$$






# Arbres

Nous allons commencer par réaliser un modèle d'arbre complet sur *l'échantillon d'apprentissage* et ce, pour avoir une vue d'ensemble.

## Arbre complet

&nbsp;

&nbsp;

```{r}
set.seed(124)

control.max <- rpart.control(cp = 0, max.depth = 0, minbucket = 1, minsplit = 1,xval = 10)
tree <- rpart(Bike_Injur~. , data = df2, control = control.max, parms = list(split = "information"))

```

\newpage

```{r,fig.height=9}
set.seed(124)
tree_prunedre <- prune(tree,cp=0)
prp(tree_prunedre,col="red",main="Arbre complet")
```


Nous savons qu'un arbre complet n'est généralement pas utilisé pour les prédictions car il ne commet aucune erreur sur l'échantillon sur lequel il est construit, puisqu'il en épouse toutes les caractéristiques. Par conséquent il est difficilement généralisable. 

Nous devons donc élaguer notre arbre selon un certain niveau de complexité afin de palier à ce problème de surraprentissage. 

## Arbre élagué

Cherchons le niveau de complexité qui minimise l'erreur estimée. 

Nous allons réalisé un graphique qui nous montre le taux d'erreur en fonction de la complexité.

```{r,fig.height=3}
plotcp(tree)
```

```{r}
cp=tree$cptable[which.min(tree$cptable[,4]),1]
```

Ce graphique nous montre que la complexité qui permet de minimiser l'erreur estimée est de $0.001420455$ avec une erreur égale à $0.88$ environ. 

En appliquant ce niveau de complexité on obtient l'arbre élagué suivant :

&nbsp;

```{r}
set.seed(124)
tree_pruned <- prune(tree, cp =cp )
prp(tree_pruned,col="red", main="Arbre élagué")
```

\


Pour se convaincre de l'utilisation d'un arbre élagué par rapport à l'arbre complet on peut représenter les **courbes ROC**.



```{r}
pred_tree<- predict(tree,
                 df2[,-4], type='prob')
pred_prune <- predict(tree_pruned,
                 df2[,-4], type='prob')

```


```{r}
predi_tree<- prediction(pred_tree[,1], (pull(df2, Bike_Injur)=="Non"))
perf_tree <- performance(predi_tree, measure = "tpr", x.measure = "fpr")
predi_prune <- prediction(pred_prune[,1], (pull(df2, Bike_Injur)=="Non"))
perf_prune <- performance(predi_prune, measure = "tpr", x.measure = "fpr")
```



```{r}
#par(mfrow=c(1,2))
#performance(predi_tree, measure = "acc") %>% plot(col="green",main="Arbre complet")
#performance(predi_prune, measure = "acc") %>% plot(col="green",main="Arbre élagué")

```

Une courbe ROC (receiver operating characteristic) est une représentation graphique de la relation qu'il existe entre la sensibilité (qui mesure sa capacité à donner un résultat positif lorsqu'une hypothèse est vérifiée) et la spécificité (qui mesure la capacité d'un test à donner un résultat négatif lorsque l'hypothèse est vérifiée) d'un test pour chaque valeur seuil considérée. 



```{r,fig.height=3.4}
r1<-data.frame(pred_model1 = pred_prune[,1], obs = (pull(df2, Bike_Injur)=="Non")) %>%
  ggplot()+aes(d=obs,m=pred_model1)+geom_roc() + theme_light()+ggtitle("Arbre élagué") 

r2<-data.frame(pred_model1 = pred_tree[,1], obs = (pull(df2, Bike_Injur)=="Non")) %>%
  ggplot()+aes(d=obs,m=pred_model1)+geom_roc() + theme_light()+ggtitle("Arbre complet")
grid.arrange(r2,r1,ncol=2)
```


Comme prévu, la courbe ROC de l'arbre complet est parfaite à cause du surapprentissage. La courbe ROC de l'arbre élagué est près de la bissectrice mais suffisament au dessus pour conclure que l'utilisation de l'arbre élagué est préférable à celui d'un classificateur aléatoire .

&nbsp; 

Nous poursuivrons la comparaison en appliquant ces deux modèles à l'échantillon d'apprentissage

&nbsp;
## Application à l'échantillon test

En appliquant ces modèles à *l'échantillon Test* on obtient la matrice de confusion suivante : 

```{r}
test_data$Bike_Injur<-as.factor(test_data$Bike_Injur)
c1<-confusionMatrix(data = predict(tree, test_data[,-4], type='class'),reference = pull(test_data,Bike_Injur), 
                    mode = "everything", positive = "Non")

c2<-confusionMatrix(data = predict(tree_pruned, test_data[,-4], type='class'),
                reference = pull(test_data,Bike_Injur), 
                mode = "everything", positive = "Non")


vide<-c("","")


conf<-cbind(c1$table ,vide,c2$table)
colnames(conf)<-c("Non","Oui"," ","Non","Oui")

kable(conf,align = "c",caption = "Matrice de confusion Arbres : test") %>%
kable_styling(latex_options=c("striped","hold_position"),full_width = F,font_size = 16) %>%
add_header_above(c("Prédiction ", "Arbre complet" = 2," ","Arbre élagué" = 2))%>%
add_header_above(c(" ", "Référence" = 5))


```
\newpage

\begin{center}
$\textbf{Courbe ROC}$
\end{center}

```{r}
set.seed(124)
k<-predict(tree,test_data)
pl<-test_data
pl$Bike_Injur<-as.factor(pl$Bike_Injur)
pl$Bike_Injur<-as.numeric(pl$Bike_Injur)
pl$Bike_Injur[pl$Bike_Injur==2 ]<-0
pl$Bike_Injur[pl$Bike_Injur==1 ]<-1
pred<-prediction(k[,1],pl$Bike_Injur)


```

```{r}
set.seed(124)
k<-predict(tree_pruned,test_data)
pl<-test_data
pl$Bike_Injur<-as.factor(pl$Bike_Injur)
pl$Bike_Injur<-as.numeric(pl$Bike_Injur)
pl$Bike_Injur[pl$Bike_Injur==2 ]<-0
pl$Bike_Injur[pl$Bike_Injur==1 ]<-1
pred1<-prediction(k[,1],pl$Bike_Injur)
```

```{r,fig.height=3.4}
set.seed(124)
par(mfrow=c(1,2))
performance(pred,"tpr","fpr")%>%plot(col="black",main="Arbre complet")
performance(pred1,"tpr","fpr")%>%plot(col="black",main="Arbre élagué")

```

```{r}
a1<-performance(pred,"auc")@y.values[[1]]
a2<-performance(pred1,"auc")@y.values[[1]]
AUC<-c(a1,a2)
```


Les mesures de performance de nos arbres appliqués à l'échantillon test sont renseignés dans le tableau suivant : 
```{r}
T<-rbind(c1$byClass[c(1,2,5)],c2$byClass[c(1,2,5)])
Ta<-cbind(T,AUC)
rownames(Ta)<-c("Arbre complet","Arbre élagué")
kable(Ta,align = "c",caption = "mesures de performance") %>%
kable_styling(latex_options=c("striped","hold_position"),full_width = F,font_size = 16) 


```

$\textbf{Conclusion:}$ Il semblerait, dans le cas particulier de notre base de données, que l'arbre complet et l'arbre élagué aient des performances qui se valent.

&nbsp;

# Bagging et Random forest

Dans cette partie nous allons comparer la méthode *RandomForest* classique avec celle du *Bagging*. Nous allons d'abord construire ces modèles sur les données d'apprentissage, comparer les résultats et enfin nous les appliquerons sur nos données test avant de conclure.

## Random forest

L'algorithme Random Forest, est l’un des plus couramment utilisé, il s’agit d’un type spécial de bagging appliqué aux arbres de décision. Cet algorithme de forêt aléatoire de Breiman(basé sur le code Fortran original de Breiman et Cutler) combine les concepts de sous-espaces aléatoires et de bagging. L'algorithme des forêts d'arbres décisionnels effectue un apprentissage sur de multiples arbres de décision entraînés sur des sous-ensembles de données légèrement différents.

Commencons par déterminer l'importance des variables selon la méthode de *Mean Decrease Accuracy*.

\newpage

```{r}
set.seed(124)
df2$Bike_Injur<-as.factor(df2$Bike_Injur)
rf <- randomForest(Bike_Injur~., data = df2, method = "class",
                   parms = list(split = "gini"), na.action = na.roughfix,
                   keep.forest = TRUE, importance = TRUE)

varImpPlot(rf, main = "Random Forest", cex = 0.8)

```

Ce graphique montre l'importance qu'a chaque variable dans la construction de notre forêt. On peut noter que les variables : *Bike_Race,Bike_Dir, Bike_Pos* et *Rd_Class* sont les plus importantes.

En construisant la forêt aléatoire sur nos données d'apprentissage on obtient :

```{r}
kable(rf$confusion,align = "c",caption = "Matrice de confusion Random Forest : apprentissage") %>%
kable_styling(latex_options=c("striped","bordered","hold_position"),full_width = F,font_size = 16)%>%
add_header_above(c("Prédiction ", "Référence" = 2, "")) 
```

Cette méthode permet d'obtenir une erreur OOB ( Out Of Bag ) de 42.96%.
\newpage
## Bagging

*Bagging* signifie bootstrap aggregating. C'est un méta-algorithme d'ensemble d'apprentissage automatique conçu pour améliorer la stabilité et la précision des algorithmes d'apprentissage automatique utilisés dans la classification statistique et la régression.

```{r}
set.seed(124)
df2$Bike_Injur<-as.factor(df2$Bike_Injur)
nvar <- ncol(df2)-1
bag <- randomForest(Bike_Injur~., data = df2, method = "class",
                    parms = list(split="gini"), mtry = nvar,na.action = na.roughfix)
```



```{r}

kable(bag$confusion,align = "c",caption = "Matrice de confusion bagging : apprentissage") %>%
kable_styling(latex_options=c("striped","bordered","hold_position"),full_width = F,font_size = 16)%>%
add_header_above(c("Prédiction ", "Référence" = 2, "")) 
```

En utilisant le Bagging : on obtient une erreur OOB ( Out Of Bag ) égale à 44.6% sur nos données d'entraînement.




## Comparaison bagging  et Random forest

Comparons les niveaux d'erreurs des deux modèles. 
Pour 500 arbres, on a : 

&nbsp; 

```{r}
set.seed(124)
data.frame(tree=1:500,
           RandomForest=rf$err.rate[,1],
           Bagging=bag$err.rate[,1])%>%
  pivot_longer(cols = c(2,3),
           names_to="method",
           values_to="error")%>%
  ggplot()+aes(y=error,x=tree,
               color=method)+
  geom_line()+theme_minimal()+ coord_cartesian(ylim = c(0.35, 0.50))

```


En comparant les erreurs sur les données d'entraînement des 2 modèles pour 500 arbres, on remarque que l'erreur provenant du *Random Forest* est plus faible que l'erreur du bagging quelque soit le nombre d'arbre considéré. 


## Application des modèles à l'échantillon Test

```{r}

rftest<-confusionMatrix(data = predict(rf, test_data[,-4], type='class'),
                reference = pull(test_data,Bike_Injur), 
                mode = "everything", positive = "Non")

bagtest<-confusionMatrix(data = predict(bag, test_data[,-4], type='class'),
                reference = pull(test_data,Bike_Injur), 
                mode = "everything", positive = "Non")


```

```{r}

conf<-cbind(rftest$table ,vide,bagtest$table)
colnames(conf)<-c("Non","Oui"," ","Non","Oui")

kable(conf,align = "c",caption = "Matrice de confusion test") %>%
kable_styling(latex_options=c("striped","hold_position"),full_width = F,font_size = 16) %>%
add_header_above(c("Prédiction ", "Random Forest" = 2," ","Bagging" = 2))%>%
add_header_above(c(" ", "Référence" = 5))

```
&nbsp; 

Les mesures de performances de ce modèle sont les suivantes :

```{r}
set.seed(124)
k<-predict(rf,test_data[,-4],type='prob')
pl<-test_data
pl$Bike_Injur<-as.factor(pl$Bike_Injur)
pl$Bike_Injur<-as.numeric(pl$Bike_Injur)
pl$Bike_Injur[pl$Bike_Injur==2 ]<-0
pl$Bike_Injur[pl$Bike_Injur==1 ]<-1
predrf<-prediction(k[,1],pl$Bike_Injur)
rfauc<-performance(predrf,"auc")@y.values[[1]]

k<-predict(bag,test_data[,-4],type='prob')
pl<-test_data
pl$Bike_Injur<-as.factor(pl$Bike_Injur)
pl$Bike_Injur<-as.numeric(pl$Bike_Injur)
pl$Bike_Injur[pl$Bike_Injur==2 ]<-0
pl$Bike_Injur[pl$Bike_Injur==1 ]<-1
prebag<-prediction(k[,1],pl$Bike_Injur)
bagauc<-performance(prebag,"auc")@y.values[[1]]
```


```{r}
T1<-rbind(c(rftest$byClass[c(1,2,5)],rfauc),c(bagtest$byClass[c(1,2,5)],bagauc))
T1<-cbind(T1)
rownames(T1)<-c("Random Forest","Bagging")
colnames(T1)<-c( "Sensitivity" ,"Specificity", "Precision" ,"AUC")
kable(T1,align = "c",caption = "Mesures de performance") %>%
kable_styling(latex_options=c("striped","hold_position"),full_width = F,font_size = 16) 

```





$\textbf{Conclusion}:$ Le Random forest permet un taux de précision superieur à celui du Random Forest.De plus

&nbsp; 

# Modèle randomForest automatisé


```{r,include=F}
set.seed(124)
### modele de 50 100 200 noeux avec 1 4 8 23 variables qu'on regarde et foret de  10 50 100
#nvar <- ncol(df2)-1
system.time(rf.tune <- tune.randomForest(x = df2[complete.cases(df2),-4],y = df2[complete.cases(df2),4],nodesize = c(10,50,100,150), mtry = c(1,4,8,nvar), ntree = c(50,100,150,250)))

formattable (rf.tune$best.parameters)
formattable(rf.tune$best.model$confusion)

```

En laissant le logiciel déterminer les paramètres optimaux pour le modèle Random Forest, on obtient les paramètres suivants : 100 noeuds, une variable et 250 arbres. 

En appliquant ces paramètres et en construisant le modèle sur notre échantillon d'apprentissage, on obtient la matrice de confusion : 


```{r}
set.seed(124)
df2$Bike_Injur<-as.factor(df2$Bike_Injur)
finalrf<-randomForest(Bike_Injur~., data = df2, method = "class",
                   parms = list(split = "gini"), na.action = na.roughfix,
                   keep.forest = TRUE, importance = TRUE,mtry=1,ntree=250,nodesize=100)
test_data$Bike_Injur<-as.factor(test_data$Bike_Injur)

c4<-confusionMatrix(data = predict(finalrf,newdata=test_data[,-4],type="class"),
                reference = pull(test_data,Bike_Injur), 
                mode = "everything", positive = "Non")


kable(finalrf$confusion,align = "c",caption = "Matrice de confusion Random Forest automatisé : apprentissage") %>%
kable_styling(latex_options=c("striped","bordered","hold_position"),full_width = F,font_size = 16)%>%
add_header_above(c( "Prédiction ","Référence" = 2," "))


```

## Application à l'échantillon test 

En appliquant ce modèle à l'échantillon test on obtient :
```{r}

kable(c4$table,align = "c",caption = "Matrice de confusion Random Forest automatisé : test") %>%
kable_styling(latex_options=c("striped","bordered","hold_position"),full_width = F,font_size = 16)%>%
add_header_above(c( "Prédiction ","Référence" = 2))
```

Les mesures de performances de ce modèle sont les suivantes :
```{r}
k<-predict(finalrf,test_data[,-4],type='prob')
pl<-test_data
pl$Bike_Injur<-as.factor(pl$Bike_Injur)
pl$Bike_Injur<-as.numeric(pl$Bike_Injur)
pl$Bike_Injur[pl$Bike_Injur==2 ]<-0
pl$Bike_Injur[pl$Bike_Injur==1 ]<-1
predrfa<-prediction(k[,1],pl$Bike_Injur)
aucrfa<-performance(predrfa,"auc")@y.values[[1]]

Tf<-rbind(c(c4$byClass[c(1,2,5)],aucrfa))
Tf<-cbind(Tf)
rownames(Tf)<-"Random Forest automatisé"
colnames(Tf)<-c( "Sensitivity" ,"Specificity", "Precision" ,"AUC")
kable(Tf,align = "c",caption = "Mesures de performance") %>%
kable_styling(latex_options=c("striped","hold_position"),full_width = F,font_size = 16) 


```

$\textbf{Conclusion}:$ Les performances de ce modèle semblent bien meilleures que tous les précédents.

&nbsp;

&nbsp;

# Boosting

Le *Boosting* génère une séquence de modèles de classification, chaque modèle de classification successif dans la séquence permettant de mieux prévoir la classification des observations qui était mal classée par les modèles de classification précédents. Lors du déploiement, les prévisions issues des différents modèles de classification pourront alors être combinées afin d'obtenir la meilleure prévision ou classification. Le Boosting, est similaire à la méthode bagging. En revanche, les étapes se produisent de manière séquentielles et non pas simultanées.


```{r}
set.seed(124)
df2$Bike_Injur<-as.factor(df2$Bike_Injur)
test_data$Bike_Injur<-as.factor(test_data$Bike_Injur)

boost<- ada(Bike_Injur~.,data=df2, type = "discrete", loss = "exponential", 
            control = rpart.control(cp = 0), iter = 200, nu = 1)

boostpen1 <- ada(Bike_Injur~.,data=df2, type = "discrete", loss = "exponential", 
            control = rpart.control(maxdepth = 1, cp = -1, minsplit = 0, xval = 0),
            iter = 200, nu = 0.1)

boostpen01 <- ada(Bike_Injur~.,data=df2, type = "discrete", loss = "exponential", 
            control = rpart.control(maxdepth = 1, cp = -1, minsplit = 0, xval = 0),
            iter = 200, nu = 0.01)

boostpen001 <- ada(Bike_Injur~.,data=df2, type = "discrete", loss = "exponential", 
            control = rpart.control(maxdepth = 1, cp = -1, minsplit = 0, xval = 0),
            iter = 200, nu = 0.001)


```

En construisant le boosting, on obtient la matrice de confusion suivante :

```{r}
kable(boost$confusion,align = "c",caption = "Matrice de confusion apprentissage : Boosting") %>%
kable_styling(latex_options=c("striped","bordered","hold_position"),full_width = F,font_size = 12)%>%
add_header_above(c( "Référence", "Prédiction " = 2))
```
L'erreur Out of bag est de 0.019. 

On constate, comme on pouvait s'y attendre, un problème de surapprentissage.
On décide donc d'appliquer différentes pénalisations à notre modèle et ce afin de trouver le modèle optimal.

Les pénalisations généralement utilisées sur le Boosting sont : $0.1$, $0.01$ et $0.001$.

$\textbf{Pénalisation = 0.1}$

```{r}
kable(boostpen1$confusion,align = "c",caption = "Matrice de confusion apprentissage : Boosting p=0.1") %>%
kable_styling(latex_options=c("striped","bordered","hold_position"),full_width = F,font_size = 12)%>%
add_header_above(c( "Référence", "Prédiction " = 2))
```
L'erreur Out of bag est de 0.398. 

$\textbf{Pénalisation = 0.01}$

```{r}
kable(boostpen01$confusion,align = "c",caption = "Matrice de confusion apprentissage : Boosting p=0.01") %>%
kable_styling(latex_options=c("striped","bordered","hold_position"),full_width = F,font_size = 12)%>%
add_header_above(c( "Référence", "Prédiction " = 2))
```
L'erreur Out of bag est de 0.408.

$\textbf{Pénalisation = 0.001}$

```{r}
kable(boostpen001$confusion,align = "c",caption = "Matrice de confusion apprentissage : Boosting p=0.001") %>%
kable_styling(latex_options=c("striped","bordered","hold_position"),full_width = F,font_size = 12)%>%
add_header_above(c( "Référence", "Prédiction " = 2))
```

L'erreur Out of bag est de 0.409.

Afin de déterminer le modèle boosting optimal, on va appliquer chacun d'entre eux à l'échantillon test. Le meilleur modèle sera celui permettant la meilleure prévision. 


## Application à l'échantillon test

Pour déterminer le meilleur modèle Boosting on va comparer les courbes ROC ainsi que les mesures de performance.


```{r,include=F}

modeleb<- list(boosting_normal = boost,
               bossting_pen1 = boostpen1,
               bossting_pen01 = boostpen01,
               bossting_pen001= boostpen001
               
               )
library(MLmetrics)
mPred <- as.data.frame(lapply(modeleb, FUN = function(m){ predict(m, test_data[,-4], type='prob')[,2]}))

conftest<-lapply(as.data.frame((mPred>0.5)*1), FUN = ConfusionMatrix, y_true = pull(test_data,Bike_Injur)=="Non")


```


```{r}
modelsPredictions <- lapply(modeleb, FUN = function(m){predict(m, test_data[,-4], type='prob')[,1]})

mPredb <- as.data.frame(modelsPredictions)
```

```{r}
mPredb %>% mutate( obs = (pull(test_data, Bike_Injur)=="Non")*1) %>% 
  pivot_longer(cols = 1:4, names_to ="method", values_to = "pred") %>%
  ggplot()+aes(d=obs,m=pred,color=method)+geom_roc() + style_roc() +ggtitle("Courbe de ROC boosting")
```

```{r}
set.seed(124)
k<-predict(boost,test_data[,-4],type='prob')
pl<-test_data
pl$Bike_Injur<-as.factor(pl$Bike_Injur)
pl$Bike_Injur<-as.numeric(pl$Bike_Injur)
pl$Bike_Injur[pl$Bike_Injur==2 ]<-0
pl$Bike_Injur[pl$Bike_Injur==1 ]<-1
pred1<-prediction(k[,1],pl$Bike_Injur)

k<-predict(boostpen1,test_data[,-4],type='prob')
pl<-test_data
pl$Bike_Injur<-as.factor(pl$Bike_Injur)
pl$Bike_Injur<-as.numeric(pl$Bike_Injur)
pl$Bike_Injur[pl$Bike_Injur==2 ]<-0
pl$Bike_Injur[pl$Bike_Injur==1 ]<-1
pred2<-prediction(k[,1],pl$Bike_Injur)

k<-predict(boostpen01,test_data[,-4],type='prob')
pl<-test_data
pl$Bike_Injur<-as.factor(pl$Bike_Injur)
pl$Bike_Injur<-as.numeric(pl$Bike_Injur)
pl$Bike_Injur[pl$Bike_Injur==2 ]<-0
pl$Bike_Injur[pl$Bike_Injur==1 ]<-1
pred3<-prediction(k[,1],pl$Bike_Injur)


k<-predict(boostpen001,test_data[,-4],type='prob')
pl<-test_data
pl$Bike_Injur<-as.factor(pl$Bike_Injur)
pl$Bike_Injur<-as.numeric(pl$Bike_Injur)
pl$Bike_Injur[pl$Bike_Injur==2 ]<-0
pl$Bike_Injur[pl$Bike_Injur==1 ]<-1
pred4<-prediction(k[,1],pl$Bike_Injur)

```

Le modèle ayant l'AUC la plus grande est le modèle boosting pénalisé à 0.01. 

On pourrait avant de conclure, comparer également les mesures de performance des tests : 


             
```{r}
test_data$Bike_Injur<-as.factor(test_data$Bike_Injur)

boost<-confusionMatrix(data = predict(boost,newdata=test_data[,-4],type="class"),
                reference = pull(test_data,Bike_Injur), 
                mode = "everything", positive = "Non")

bossting_pen1<-confusionMatrix(data = predict(boostpen1,newdata=test_data[,-4],type="class"),
                reference = pull(test_data,Bike_Injur), 
                mode = "everything", positive = "Non")

bossting_pen01<-confusionMatrix(data = predict(boostpen01,newdata=test_data[,-4],type="class"),
                reference = pull(test_data,Bike_Injur), 
                mode = "everything", positive = "Non")

bossting_pen001<-confusionMatrix(data = predict(boostpen001,newdata=test_data[,-4], type="class"),reference = pull(test_data,Bike_Injur) , mode = "everything", positive = "Non")

```

```{r}
Sensibilité<-c(sensitivity(boost$table),sensitivity(bossting_pen1$table),sensitivity(bossting_pen01$table),sensitivity(bossting_pen001$table))
Spécificité<-c(specificity(boost$table),specificity(bossting_pen1$table),specificity(bossting_pen01$table),specificity(bossting_pen001$table))
Précision<-c(precision(boost$table),precision(bossting_pen1$table),precision(bossting_pen01$table),precision(bossting_pen001$table))
AUC<-c(performance(pred1,"auc")@y.values[[1]],performance(pred2,"auc")@y.values[[1]],performance(pred3,"auc")@y.values[[1]],performance(pred4,"auc")@y.values[[1]])

test1<-cbind(Sensibilité,Spécificité,Précision,AUC)
rownames(test1)<-c("Boosting","Boosting p=0.1","Boosting p=0.01","Boosting p=0.001")

kable(test1,align = "c",caption = "Comparaions des boosting") %>%
kable_styling(latex_options=c("striped","bordered","hold_position"),full_width = F,font_size = 12)
```

$\textbf{Conclusion :}$ Les AUC étant très similaires, pour le choix du meilleur modèle Boosting on se base sur la Précision : c'est donc le Boosting pénalisé à 0.01 qui est choisi.

\newpage

# Scoring


```{r,echo=F,include=F}
data$Crash_Hour[data$Crash_Hour>= 8 & data$Crash_Hour<=19]<-"Jour"
data$Crash_Hour[data$Crash_Hour< 8 & data$Crash_Hour>19 ]<-"Nuit"
data$Crash_Hour[data$Crash_Hour<=1 ]<-"Nuit"



data$Bike_Age[data$Bike_Age<14 ]<-"Enfant"
data$Bike_Age[data$Bike_Age>= 14 & data$Bike_Age<= 30]<-"Jeune"
data$Bike_Age[data$Bike_Age> 30 & data$Bike_Age<= 70]<-"Adulte"


```


```{r,echo=F}
newdataqual<-data[,-c(1,2,4,12,13,15,16,17,20,21,22,23,24,25,27,28,30,34,35,38,40,41,44,45,46,47,50,51,52,54)]

```

&nbsp;

```{r,echo=F,include=F}
table(data$Bike_Race)
newdataqual$Bike_Race<-as.character(newdataqual$Bike_Race)
newdataqual$Bike_Race[newdataqual$Bike_Race=="Native American"]<-"Other"
newdataqual$Bike_Race[newdataqual$Bike_Race=="Asian"]<-"Other"
newdataqual$Bike_Race[newdataqual$Bike_Race=="Hispanic"]<-"Other"
newdataqual$Bike_Race[newdataqual$Bike_Race=="/Missing"]<-"Other"
table(newdataqual$Bike_Race)

table(data$Num_Lanes)
newdataqual$Num_Lanes<-as.character(newdataqual$Num_Lanes)
newdataqual$Num_Lanes[newdataqual$Num_Lanes=="1 lane"]<-"Autre"
newdataqual$Num_Lanes[newdataqual$Num_Lanes=="2 lanes"]<-"2 voie"
newdataqual$Num_Lanes[newdataqual$Num_Lanes=="3 lanes"]<-"3_4_5 Voie"
newdataqual$Num_Lanes[newdataqual$Num_Lanes=="4 lanes"]<-"3_4_5 Voie"
newdataqual$Num_Lanes[newdataqual$Num_Lanes=="5 lanes"]<-"3_4_5 Voie"
newdataqual$Num_Lanes[newdataqual$Num_Lanes=="6 lanes"]<-"Autre"
newdataqual$Num_Lanes[newdataqual$Num_Lanes=="7 lanes"]<-"Autre"
newdataqual$Num_Lanes[newdataqual$Num_Lanes=="8 lanes"]<-"Autre"
newdataqual$Num_Lanes[newdataqual$Num_Lanes=="9 or more lanes"]<-"Autre"

table(newdataqual$Num_Lanes)
```


```{r,echo=F,include=F}
levels(newdataqual$Bike_Injur)

newdataqual$Bike_Injur<-as.character(newdataqual$Bike_Injur)

newdataqual$Bike_Injur[newdataqual$Bike_Injur=="A: Disabling Injury"]<-1
newdataqual$Bike_Injur[newdataqual$Bike_Injur=="B: Evident Injury"]<-1
newdataqual$Bike_Injur[newdataqual$Bike_Injur=="C: Possible Injury"]<-0
newdataqual$Bike_Injur[newdataqual$Bike_Injur=="Injury"]<-1
newdataqual$Bike_Injur[newdataqual$Bike_Injur=="K: Killed"]<-1
newdataqual$Bike_Injur[newdataqual$Bike_Injur=="O: No Injury"]<-0
table(newdataqual$Bike_Injur)

newdataqual$Drvr_Injur<-as.character(newdataqual$Drvr_Injur)
newdataqual$Drvr_Injur[newdataqual$Drvr_Injur=="A: Disabling Injury"]<-"1"
newdataqual$Drvr_Injur[newdataqual$Drvr_Injur=="B: Evident Injury"]<-"1"
newdataqual$Drvr_Injur[newdataqual$Drvr_Injur=="C: Possible Injury"]<-"0"
newdataqual$Drvr_Injur[newdataqual$Drvr_Injur=="Injury"]<-"1"
newdataqual$Drvr_Injur[newdataqual$Drvr_Injur=="K: Killed"]<-"1"
newdataqual$Drvr_Injur[newdataqual$Drvr_Injur=="O: No Injury"]<-"0"
table(data$Drvr_Injur)



table(newdataqual$Bike_Pos)
newdataqual$Bike_Pos<-as.character(newdataqual$Bike_Pos)
newdataqual$Bike_Pos[newdataqual$Bike_Pos=="Bike Lane / Paved Shoulder"]<-"Other"
newdataqual$Bike_Pos[newdataqual$Bike_Pos=="Driveway / Alley"]<-"Other"
newdataqual$Bike_Pos[newdataqual$Bike_Pos=="Multi-use Path"]<-"Other"





levels(newdataqual$Light_Cond) 
newdataqual$Light_Cond<-as.character(newdataqual$Light_Cond)

newdataqual$Light_Cond[newdataqual$Light_Cond=="Dark -  Lighting"]<-"Tres_Sombre"
newdataqual$Light_Cond[newdataqual$Light_Cond=="Dark - Lighted Roadway"]<-"Sombre"
newdataqual$Light_Cond[newdataqual$Light_Cond=="Dark - Roadway Not Lighted"]<-"Tres_Sombre"
newdataqual$Light_Cond[newdataqual$Light_Cond=="Dawn"]<-"Sombre"
newdataqual$Light_Cond[newdataqual$Light_Cond=="Daylight"]<-"Clair"
newdataqual$Light_Cond[newdataqual$Light_Cond=="Dusk"]<-"Sombre"
newdataqual$Light_Cond[newdataqual$Light_Cond=="Other"]<-"Sombre"
table(newdataqual$Light_Cond)


levels(newdataqual$Rd_Conditi)
newdataqual$Rd_Conditi<-as.character(newdataqual$Rd_Conditi)
newdataqual$Rd_Conditi[newdataqual$Rd_Conditi=="Water (Standing, Moving)"]<-"Wet"
newdataqual$Rd_Conditi[newdataqual$Rd_Conditi=="Sand, Mud, Dirt, Gravel"]<-"Other"
newdataqual$Rd_Conditi[newdataqual$Rd_Conditi=="Ice"]<-"Other"
newdataqual$Rd_Conditi[newdataqual$Rd_Conditi=="Snow"]<-"Other"
newdataqual$Rd_Conditi[newdataqual$Rd_Conditi=="Other"]<-"Wet"
table(newdataqual$Rd_Conditi)



levels(newdataqual$Drvr_Race)
newdataqual$Drvr_Race<-as.character(newdataqual$Drvr_Race)
newdataqual$Drvr_Race[newdataqual$Drvr_Race=="Native American"]<-"Other"
newdataqual$Drvr_Race[newdataqual$Drvr_Race=="Asian"]<-"Other"
newdataqual$Drvr_Race[newdataqual$Drvr_Race=="Hispanic"]<-"Other"

table(newdataqual$Drvr_Race)

colnames(newdataqual)
newdataqual$Workzone_I<-as.character(newdataqual$Workzone_I)
newdataqual$Workzone_I[newdataqual$Workzone_I=="No0"]<-"No"
table(newdataqual$Workzone_I)

levels(newdataqual$Developmen)
table(newdataqual$Developmen)
newdataqual$Developmen<-as.character(newdataqual$Developmen)
newdataqual$Developmen[newdataqual$Developmen=="Industrial"]<-"Farms_Pastures"
newdataqual$Developmen[newdataqual$Developmen=="Farms, Woods, Pastures"]<-"Farms_Pastures"


levels(newdataqual$Rd_Class)
table(newdataqual$Rd_Class)
newdataqual$Rd_Class<-as.character(newdataqual$Rd_Class)
newdataqual$Rd_Class[newdataqual$Rd_Class=="Interstate"]<-"US Route"
newdataqual$Rd_Class[newdataqual$Rd_Class=="Private Road, Driveway"]<-"Local Street"
newdataqual$Rd_Class[newdataqual$Rd_Class=="NC Route"]<-"US Route"
newdataqual$Rd_Class[newdataqual$Rd_Class=="State Secondary Route"]<-"US Route"
newdataqual$Rd_Class[newdataqual$Rd_Class=="Public Vehicular Area"]<-"US Route"




newdataqual$Bike_Pos<-as.character(newdataqual$Bike_Pos)
newdataqual$Bike_Pos[newdataqual$Bike_Pos=="Non-Roadway"]<-"Other"
newdataqual$Bike_Pos[newdataqual$Bike_Pos=="Sidewalk / Crosswalk / Driveway Crossing"]<-"Other"

newdataqual$Bike_Alc_D<-as.character(newdataqual$Bike_Alc_D)
newdataqual$Bike_Alc_D[newdataqual$Bike_Alc_D=="Missing"]<-"Yes"


newdataqual$Crash_Hour[newdataqual$Crash_Hour>= 8 & newdataqual$Crash_Hour<=19]<-"Jour"
newdataqual$Crash_Hour[newdataqual$Crash_Hour<8 & newdataqual$Crash_Hour>19 ]<-"Nuit"

table(newdataqual$Crash_Loc)
newdataqual$Crash_Loc<-as.character(newdataqual$Crash_Loc)
newdataqual$Crash_Loc[newdataqual$Crash_Loc=="Intersection-Related"]<-"Intersection"
newdataqual$Crash_Loc[newdataqual$Crash_Loc=="Location"]<-"Non-Intersection"
newdataqual$Crash_Loc[newdataqual$Crash_Loc=="Non-Roadway"]<-"Non-Intersection"

table(newdataqual$Crash_Loc)
```



```{r,echo=F}
newdataqual$Bike_Injur<-as.factor(newdataqual$Bike_Injur)
newdataqual$Bike_Pos<-as.factor(newdataqual$Bike_Pos)
newdataqual$Bike_Race<-as.factor(newdataqual$Bike_Race)
newdataqual$Bike_Pos<-as.factor(newdataqual$Bike_Pos)
newdataqual$Developmen<-as.factor(newdataqual$Developmen)
newdataqual$Drvr_Injur<-as.factor(newdataqual$Drvr_Injur)
newdataqual$Drvr_Race<-as.factor(newdataqual$Drvr_Race)
newdataqual$Rd_Class<-as.factor(newdataqual$Rd_Class)
newdataqual$Rd_Conditi<-as.factor(newdataqual$Rd_Conditi)
newdataqual$Workzone_I<-as.factor(newdataqual$Workzone_I)
newdataqual$Light_Cond<-as.factor(newdataqual$Light_Cond)
newdataqual$Rural_Urba<-as.factor(newdataqual$Rural_Urba)
newdataqual$Bike_Age<-as.factor(newdataqual$Bike_Age)
newdataqual$Crash_Hour<-as.factor(newdataqual$Crash_Hour)
```

```{r}
data2<-newdataqual[complete.cases(newdataqual),-c(1,9) ]

```

Comme toute bonne démarche de modélisation, la construction d'un bon **score** se fait par une succession d'étapes : nous commencerons par vérifier la liaison entre les descripteurs, ensuite nous construirons les modèles sur l'*échantillon d'apprentissage* et enfin nous les appliquerons sur l'*échantillon test*. Nous pourrons ensuite comparer les mesures de performance afin de déterminer le meilleur modèle.

Nous comparerons plusieurs modèles et retiendrons le modèle le plus adéquat selon l'objectif de l'étude.

\

## Variable à expliquer

\
Revenons sur nos données initiales. La variable **Bike_Injur** est séparée en plusieurs modalités, comme suit : 

```{r}
bi<-data.frame(table(data$Bike_Injur))
bi1<-rbind(bi$Freq)
colnames(bi1)<-bi$Var1
```

```{r}
kable(bi1,"latex",align = "c") %>%
kable_styling(latex_options = c(" striped","hold_position"),
              position="center",full_width = F)

```

Nous allons affecter à ces modalités les valeurs 0 et 1. 

La valeur 1 : à celles qui concernent les blessures avérées et graves (Killed, Disabling Injury, Evident Injury et Injury)  et la valeur 0 à celles qui concernent l'absence, évidente ou non, de blessure (No Injury et Possible Injury). 

On obtient donc : 

```{r}
data2<-newdataqual[complete.cases(newdataqual),-c(1,9) ]

```

```{r}
re<-data.frame(table(data$Bike_Injur))
aaa<-c(0,1)
aaa1<-c(526+2199,291+2405+172+123)
C<-data.frame(rbind(aaa1))

row.names(C)<-c("Effectif")
colnames(C)<-c("0","1")

kable(C,align = "c") %>%
kable_styling(latex_options =c("striped","hold_position"),full_width = F) %>%
  add_header_above(c(" ", "Bike_Injur" = 2))
```

On souhaite déterminer la probabilité qu'un cycliste soit bléssé ou non compte tenu des paramètres de son accident.

Pour cela, nous utiliserons des **_Régression logistique_**.

Nous allons effectuer une sélection de variables et ce pour plusieurs raisons.

D'abord parce que certaines des variables de notre base de données sont inutilisables en l'état et d'autre part parce qu'un modèle avec peu de variables sera plus facilement généralisable en terme de robustesse : *Principe du rasoir d'Occam*.


\

## Analyse exploratoire

Construisons une ACM afin d'essayer d'identifier les corrélations éventuelles entre les variables explicatives.

```{r}
library("FactoMineR")
library("factoextra")
res.pca <- MCA(data2, graph = FALSE)
```


```{r,fig.height=6}
library(gridExtra)
re<-fviz_pca_var(res.pca, col.var = "black",select.var = list(cos2 =0.1),gradient.cols = c("#3399FF","#663399"),labelsize = 3)

ree<-plot(res.pca, choix = "var",title="Rapport de correlation au carre", col.var = "black",xlim=c(0,0.75),ylim=c(0,0.6),label="none")

grid.arrange(re,ree,ncol=2)
```

Plusieurs variables semblent très corrélées, ce qui pourrait entraîner des problèmes de colinéarité lors de la régression.


\


```{r}
set.seed(111)

d = sort(sample(nrow(data2), nrow(data2) * 2/3))
# Echantillon d'apprentissage
appren <- data2[d, ]
# Echantillon de test
test <- data2[-d, ]

```

\

## Construction des modèles


Le **modèle général** que nous allons construire est un modèle naïf. Il prend en compte toutes les variables, sans aucune spécification particulière : 
\begin{center}
$Bike\_Injur_i$ = $\beta_0$ + $\beta_1$ $Bike\_Age_i$ + $\beta_2$ $Bike\_Alc\_D_i$  + $\beta_3$ $Bike\_Dir_i$+ $\beta_4$ $Bike\_Pos_i$ + $\beta_5$ $Bike\_Race_i$  + $\beta_6$ $Bike\_Sex_i$ + $\beta_7$ $Crash\_Hour_i$ + $\beta_8$ $Crash\_Loc_i$ + $\beta_9$ $Developmen_i$  + $\beta_{10}$ $Drvr\_Alc\_D_i$ + $\beta_{11}$ $Drvr\_Injur_i$ + $\beta_{12}$ $Drvr\_Race_i$ + $\beta_{13}$ $Drvr\_Sex_i$ + $\beta_{14}$ $Hit\_Run_i$ + $\beta_{15}$ $Light\_Cond_i$ + $\beta_{16}$ $Num\_Lanes_i$ + $\beta_{17}$ $Rd\_Class_i$  + $\beta_{18}$ $Rd\_Conditi_i$ + $\beta_{19}$ $Region_i$ + $\beta_{20}$ $Rural\_Urba_i$ + $\beta_{21}$ $Workzone\_I_i$ + $\varepsilon_i$ 
\end{center}

&nbsp; 

Nous savons que ce modèle sera très mauvais étant donné le nombre de variable utilisée. Pour obtenir un modèle pertinent nous devons  effectuer une sélection des variables et ce pour plusieurs raisons, la principale étant que : un modèle avec peu de variables sera plus facilement généralisable en terme de robustesse *Principe du rasoir d'Occam*.

&nbsp;

Le second modèle, le **modèle AIC** : sera obtenu en faisant une sélection automatique de variables, sur le critère d'Akaïke (AIC). Ce dernier s'écrit comme suit: $AIC= 2k-2\ln(L)$ ; où k est le nombre de paramètres à estimer du modèle et L est le maximum de la fonction de vraisemblance du modèle.
Si l'on considère un ensemble de modèles candidats, le modèle choisi sera celui qui aura la plus faible valeur d'AIC.


&nbsp;

Le dernier modèle, le **modèle AIC modifié** sera issue du second. On supprimera toutes les variables non significatives du second modèle.



Les différents résultats des régréssions obtenus sont renseignés dans le tableau çi-dessous : 

```{r}
lm1<-glm(Bike_Injur~Bike_Age+ Bike_Alc_D+Bike_Dir+ Bike_Pos+ Bike_Race+ Bike_Sex+Crash_Hour+ Crash_Loc+Developmen+Drvr_Alc_D+ Drvr_Injur+ Drvr_Race+Drvr_Sex+ Rd_Conditi+Region  +Rural_Urba+Workzone_I, data=appren, family=binomial(link=logit)) 
```

```{r,include=FALSE}
require(MASS)
slm<-step(lm1, direction = "both")
```



```{r}
logit = function(formula, lien = "logit", data = NULL) {
    glm(formula, family = binomial(link = lien), data2)
}
```

```{r}
m.logit <- logit(Bike_Injur ~ Bike_Dir + Bike_Pos + Bike_Race + 
    Developmen + Drvr_Alc_D + Rural_Urba + Workzone_I, data = appren)

m.logit2 <- logit(Bike_Injur ~ Bike_Dir + Bike_Pos + Bike_Race + Drvr_Alc_D + Rural_Urba , data = appren)

```

\newpage 

```{r,warning=F,message=F,results='asis',header=F}
library(stargazer)
stargazer(m.logit,m.logit2,lm1,type="latex",header=F,font.size = "tiny",title="Résultats")
```

**Modèle AIC**:
\begin{center}
$Bike\_Injur_i$ = $\beta_0$ +  $\beta_1$ $Bike\_Alc\_D_i$  + $\beta_2$ $Bike\_Dir_i$+ $\beta_3$ $Bike\_Pos_i$ + $\beta_4$ $Bike\_Race_i$ + $\beta_9$ $Developmen_i$  + $\beta_{5}$ $Drvr\_Alc\_D_i$ + $\beta_{6}$ $Rural\_Urba_i$ + $\beta_{7}$ $Workzone\_I_i$ + $\varepsilon_i$ 
\end{center}


**Modèle AIC modifié**:
\begin{center}
$Bike\_Injur_i$ = $\beta_0$ +  $\beta_1$ $Bike\_Alc\_D_i$  + $\beta_2$ $Bike\_Dir_i$+ $\beta_3$ $Bike\_Pos_i$ + $\beta_4$ $Bike\_Race_i$ + $\beta_{5}$ $Drvr\_Alc\_D_i$ + $\beta_{6}$ $Rural\_Urba_i$ + $\varepsilon_i$ 
\end{center}

```{r}
library(GGally)
library(gridExtra)
```
&nbsp;


## Validation des modèles : Indicateurs de qualité et de robustesse


On exclu le modèle général car la plupart des coefficients ne sont pas significatifs.  

On s'intéressera donc exclusivement aux modèles AIC (Régréssion 1) et AIC modifié (Régréssion 2).

```{r,warning=F,message=F,results='asis',header=F}
stargazer(m.logit,m.logit2,type="latex",header=F,font.size = "tiny",title="Résultats des modèles 2 et 3")
```


$\textbf{Résidus de déviances}$

Pour les régressions logistiques, on s'intéresse la plupart du temps aux résidus de déviance. Ils prennent généralement des valeurs qui oscillent entre -2 et 2. 

```{r,fig.height=8.2}
par(mfrow = c(2, 1))
plot(rstudent(m.logit),main="Modèle 1", type = "p", cex = 0.5, ylab = "Résidus studentisés ", col = "black", ylim = c(-3, 3),xlab="")
abline(h = c(-2, 2), col = "red")

plot(rstudent(m.logit2),main="Modèle 2", type = "p", cex = 0.5, ylab = "Résidus studentisés ", col = "black", ylim = c(-3, 3),xlab="")
abline(h = c(-2, 2), col = "red")
```

Il semblerait qu'il n'y ait pas de valeurs aberrantes. 

Les deux modèles sont donc utilisables dans l'état.

&nbsp;

&nbsp;



```{r}
appren.p <- cbind(appren, predict(m.logit, newdata = appren, type = "link", 
    se = TRUE))
```


```{r}
appren.p <- within(appren.p, {
    PredictedProb <- plogis(fit, lower.tail = TRUE)
    LL <- plogis(fit - (1.96 * se.fit))
    UL <- plogis(fit + (1.96 * se.fit))
})
```

```{r}
appren.p <- cbind(appren.p, pred.bike_injur = factor(ifelse(appren.p$PredictedProb > 
    0.5, 1, 0)))
```

```{r}
appren.p1 <- cbind(appren, predict(m.logit2, newdata = appren, type = "link", 
    se = TRUE))

```



```{r}
appren.p1 <- within(appren.p1, {
    PredictedProb <- plogis(fit, lower.tail = TRUE)
    LL <- plogis(fit - (1.96 * se.fit))
    UL <- plogis(fit + (1.96 * se.fit))
})

```

```{r}
appren.p1 <- cbind(appren.p1, pred.bike_injur = factor(ifelse(appren.p1$PredictedProb > 
    0.5, 1, 0)))
```

```{r}
# Matrice de confusion2
m.confusion1 <- as.matrix(table(appren.p1$pred.bike_injur, appren.p1$Bike_Injur))
```

```{r}
# Matrice de confusion
m.confusion <- as.matrix(table(appren.p$pred.bike_injur, appren.p$Bike_Injur))
vide<-c("","")

conf<-cbind(m.confusion ,vide,m.confusion1)
colnames(conf)<-c("0","1","","0","1")


```





```{r}
m.confusion <- unclass(m.confusion)
# Taux d'erreur
Tx_err <- function(y, ypred) {
    mc <- table(y, ypred)
    error <- (mc[1, 2] + mc[2, 1])/sum(mc)
    
}

m.confusion1 <- unclass(m.confusion1)
# Taux d'erreur
Tx_err1 <- function(y, ypred) {
    mc <- table(y, ypred)
    error <- (mc[1, 2] + mc[2, 1])/sum(mc)
    
}

```

```{r}
tx_a1<-Tx_err(appren.p$pred.bike_injur, appren.p$Bike_Injur)

tx_a2<-Tx_err1(appren.p1$pred.bike_injur, appren.p1$Bike_Injur)

```

```{r}
tx_a3<-cbind(tx_a1,tx_a2)

colnames(tx_a3)<-c("Modèle 1","Modèle 2")


```

 
## Application  à l'échantillon test

\newpage
En appliquant les deux modèles à l'*échantillon Test* on a :

```{r}
test.p <- cbind(test, predict(m.logit, newdata = test, type = "response", se = TRUE))
test.p <- cbind(test.p, pred.bike_injur = factor(ifelse(test.p$fit > 0.5, 1, 0)))


m.confusiontest <- as.matrix(table(test.p$pred.bike_injur, test.p$Bike_Injur))


test.p1 <- cbind(test, predict(m.logit2, newdata = test, type = "response", se = TRUE))
test.p1 <- cbind(test.p1, pred.bike_injur = factor(ifelse(test.p1$fit > 0.5, 1, 0)))


m.confusiontest1 <- as.matrix(table(test.p1$pred.bike_injur, test.p$Bike_Injur))


conf1<-cbind(m.confusiontest ,vide,m.confusiontest1)
colnames(conf1)<-c("0","1","","0","1")

kable(conf1,align = "c",caption = "Matrice de confusion Scoring : test") %>%
kable_styling(latex_options=c("striped","hold_position"),full_width = F,font_size = 14) %>%
add_header_above(c("Prédiction", "Modèle AIC" = 2," ","Modèle AIC modifié" = 2))%>%
add_header_above(c(" ", "Référence" = 5))

```



\begin{center}
Qu'en est-il du taux d'erreur ?
\end{center}

```{r}
# calcul du taux d'erreur sur l'échantillon test

tx3t<-Tx_err(test.p$pred.bike_injur, test.p$Bike_Injur)

tx3t2<-Tx_err(test.p1$pred.bike_injur, test.p1$Bike_Injur)

ttx3<-cbind(tx3t,tx3t2)

colnames(ttx3)<-c("Modèle AIC","Modèle AIC modifié")


kable(ttx3,align = "c",caption = "Taux d'erreur test", escape = TRUE) %>%
kable_styling(latex_options=c("striped","hold_position"),full_width = F,font_size = 14) 
```

Le taux d'erreur est relativement le même pour les deux modèles. Il est inférieur à 0.5, ce qui est suffisant pour conclure que les modèles sont pertinents.

&nbsp;

\newpage

\begin{center}
$\textbf{Courbes ROC}$
\end{center}

```{r}
library(ROCR)
```



```{r}
Pred = prediction(appren.p$PredictedProb, appren.p$Bike_Injur)
Perf = performance(Pred, "tpr", "fpr")
perf <- performance(Pred, "auc")

Predtest = prediction(test.p$fit, test.p$Bike_Injur)
Perftest = performance(Predtest, "tpr", "fpr")
perftest <- performance(Predtest, "auc")

```




```{r}
Pred1 = prediction(appren.p1$PredictedProb, appren.p1$Bike_Injur)
Perf1 = performance(Pred1, "tpr", "fpr")
perf1 <- performance(Pred1, "auc")

Predtest1 = prediction(test.p1$fit, test.p1$Bike_Injur)
Perftest1 = performance(Predtest1, "tpr", "fpr")
perftest1 <- performance(Predtest1, "auc")

par(mfrow = c(1, 2))
plot(Perftest1, colorize = TRUE, main = "Modèle AIC ")
plot(Perftest, colorize = TRUE, main = "Modèle AIC modifié ")

```




\begin{center}
\textbf{Aire sous les courbes}
\end{center}

```{r}

airr1<-data.frame(cbind(perftest@y.values[[1]],perftest1@y.values[[1]]))

colnames(airr1)<-c("Modèle AIC","Modèle AIC modifié")

kable(airr1,align = "c",caption = "Aire sous les courbes", escape = TRUE) %>%
kable_styling(latex_options=c("striped","hold_position"),full_width = F,font_size = 16) 
```

L'aire sous les courbes est relativement le même lui aussi. 

\newpage 

$\textbf{Mesure de performance :}$

```{r}
library(MLmetrics)

k<-cbind(sensitivity(m.confusiontest),specificity(m.confusiontest),precision(m.confusiontest))
rownames(k)<-"Modèle AIC"

k1<-cbind(sensitivity(m.confusiontest1),specificity(m.confusiontest1),precision(m.confusiontest1))
rownames(k1)<-"Modèle AIC modifié"

klm<-rbind(k,k1)
colnames(klm)<-c("Sensibilité","Spécificité","Précision")
kable(klm,align = "c",caption = "Mesure de performance : Scoring", escape = TRUE) %>%
kable_styling(latex_options=c("striped","hold_position"),full_width = F,font_size = 16) 
```


$\textbf{Conclusion :}$ Les deux modèles sont très similaires et ont des mesures de performances qui se valent mais on décide de privilégié celui étant le plus précis : le modèle AUC modifié.


# Choix du meilleur modèle


```{r}
rr<-cbind(sensitivity(bossting_pen01$table),specificity(bossting_pen01$table),precision(bossting_pen01$table),performance(pred3,"auc")@y.values[[1]])



L1<-cbind(sensitivity(m.confusiontest1),specificity(m.confusiontest1),precision(m.confusiontest1),perftest1@y.values[[1]])
rownames(L1)<-"Scoring"

Llm<-rbind(L1)


rownames(rr)<-"Boosting p=0.01"

comp<-rbind(Ta,T1,Tf,rr,Llm)
kable(comp,align = "c",caption = "Comparaison des modèles", escape = TRUE) %>%
kable_styling(latex_options=c("striped","hold_position"),full_width = F,font_size = 16)
```


# Discussion

La sélection des mesures de performance les plus pertinentes se fait en fonction de la problématique à traiter. La notre pourrait être soit de déterminer la probabilité que l’accident ait causé une blessure soit de déterminé la probabilité que l’accident n’ait pas causé de blessure. 
Il n’y a pas grand intérêt à déterminer la probabilité que l’accident ait causé une blessure pour les hôpitaux par exemple car ces derniers envoient systématiquement une ambulance. En revanche il serait très intéressant pour eux de déterminer la probabilité que l’accident n’ait pas causé de blessures et ce afin de gérer le flux des ambulances ou simplement de faire un choix prioritaire parmi deux situations par exemple. En ce sens la sensibilité ne donne pas une mesure très pertinente. 
La spécificité en revanche, qui mesure le taux de vrais négatifs (dans notre cas qu’il n’y ait pas de blessure) semble plus pertinente.

Aussi, La précision est de fait une mesure très intéressante. Dans le cadre de notre étude : c’est la capacité de nos modèles à ne prédire non à une blessure si l’accident n’a effectivement pas entraîner une blessure. Pour évaluer un compromis entre sensibilité et précision, on peut calculer la "F-mesure", qui est leur moyenne harmonique.


Le calcul du F-mesure ou encore F-score est le suivant:

$$F−mesure=\frac{2\times Précision\times sensibilité}{précision+sensibilité}$$

```{r}
FarbreC=(2 *Ta[1,1] * Ta[1,3]) / (Ta[1,1] + Ta[1,3])
FarbreE=(2 *Ta[2,1] * Ta[2,3]) / (Ta[2,1] + Ta[2,3])
Frandom=(2 *T1[1,1] * T1[1,3]) / (T1[1,1] + T1[1,3])
Fbagging=(2 *T1[2,1] * T1[2,3]) / (T1[2,1] + T1[2,3])
FrandomFA=(2 *Tf[1] * Tf[3]) / (Tf[1] + Tf[3])
Fboosting=(2 *rr[1] * rr[3]) / (rr[1] + rr[3])
Fscoring=(2 * sensitivity(m.confusiontest1) * precision(m.confusiontest1) ) / (sensitivity(m.confusiontest1)  +  precision(m.confusiontest1) )

tableau= rbind(FarbreC,FarbreE,Frandom,Fbagging,FrandomFA,Fboosting,Fscoring)
rownames(tableau)=c("Arbre complet","Arbre élagué","Random Forest","Bagging","Random Forest automatisé","Boosting p=0.01","Scoring")
colnames(tableau)="F-mesure"
kable(tableau,align = "c",caption = "Comparaison des F-mesure des modèles", escape = TRUE) %>%
kable_styling(latex_options=c("striped","hold_position"),full_width = F,font_size = 16)

```


Ces considérations ainsi prises en compte nous permettent de conclure que le modèle RandomForest automatisé (qui pour rappel est un modèle Random Forest ayant les paramètres suivant : 250 arbres, 50 nœud et 1 ) et le Boosting pénalisé à 0.01 sont les modèles permettant d’obtenir les meilleures résultats pour la problématique : Quelle est la probabilité que l’accident n’est pas engendré de blessures.