Labo_exemple_Pima.Rmd

---
title: "Labo 23 mars"
author: "Anne-Sophie Charest/Stéphane Caron"
date: "22 mars 2018"
output: html_document
---

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

## Le problème

On utilise le jeu de données PimaIndiansDiabetes2 de la librairie mlbench. On souhaite prédire si un individu souffre de diabète à partir de 8 variables explicatives. Votre travail est d'ajuster quelques classifieurs vus en classe, de sélectionner correctement les hyperparamètre le cas échéant, et d'identifier le meilleur classifieur sur un jeu de données de validation. Je vous invite à suivre les étapes ci-bas pour ce faire

### Explorer les données 

Prenez quelques minutes pour analyser le jeu de données. Voici une liste de questions que vous pouvez explorer à l'aide de statistiques descriptives.

* Quelle est la taille du jeu de données?
* Quelle est la prévalence du diabète dans le jeu de données?
* De quel type sont les différentes variables explicatives?
* Semble-t-il y avoir une corrélation entre les différentes variables explicatives?
* Quelles variables semblent les plus utiles pour prédire le diabète?
* Y a-t-il des données manquantes?


```{r data_exploration}
library(mlbench)
library(data.table)
data("PimaIndiansDiabetes2")

#On renomme le jeu de données
dat = data.table(PimaIndiansDiabetes2)

summary(dat)
dim(dat)
table(dat$diabetes)
apply(apply(dat, 2, is.na), 2, sum)/nrow(dat)
```

### Tirer un jeu de données de validation 

On conservera 50 observations au hasard pour tester les différents classifieurs obtenus. J'ai déjà tiré ces 50 observations pour qu'on teste tous sur le même jeu de données.   

```{r }
#Le code que j'ai utilisé pour choisir les 50 observations; ne pas rouler
#set.seed(15)
#idx = sample(768,50,replace = F)
#a = cat(idx, sep = ",")
#a = paste(idx,collapse = ",") 
#b = paste( "c(", a, ")", collapse = "")

#Les index des 50 observations à conserver dans le jeu de validation
#Vous n'avez qu'a
index_test = c( 463,150,741,498,281,755,622,194,523,632,80,490,385,534,651,634,337,725,106,582,602,593,267,44,421,757,758,110,687,353,368,190,362,87,377,483,89,378,221,555,667,333,669,188,249,311,239,543,608,8 )

dat_test = dat[index_test,] #À garder pour la fin
dat_train = dat[-index_test,] #À utiliser pour obtenir les classifieurs
```

### Obtenir différents classifieurs 

L'idée ici est d'ajuster différents modèles vus en classe sur le jeu de données. Vous devrez parfois choisir un ou des hyperparamètres pour le modèle; utilisez dans ce cas la validation croisée. Voici une liste des classifieurs vus en classe:

* Méthode du plus proche voisin (k à choisir)
* Classifieur naif de Bayes (loi normale ou estimation de densité par noyau)
* Régression logistique (choisir variables à inclure dans le modèle et point de coupure)
* Analyse discriminante linéaire
* Analyse discriminante quadratique 
* Arbre de régression (nombre de divisions et critère de division à choisir)
* Forêt aléatoire (plusieurs choix possibles, incluant le nombre d'arbres, et le nombre de variables à tester à chaque division)

Je vous invite à commencer par les arbres et les forêts aléatoires puisque vous n'avez pas encore utilisé ces fonctions avec R. Testez ensuite les modèles de votre choix. 

Notez que dans certains cas vous aurez des problèmes avec les données manquantes. Vous pouvez simplement décider de ne pas utiliser ces classifieurs, ou alors penser à une façon de gérer les données manquantes. 

```{r cv_function}
create_folder <- function(k, data){
  
  n <- nrow(data)
  # Shuffle the data
  data <-  data[sample(1:n, replace = F),]
  # Create the groups
  data$group <- cut(seq(1, nrow(data)), breaks = 10, labels = FALSE)

  return(data)
}

```



```{r classfiers}

set.seed(666) #Number of the beast

# Create cross-validation data
k = 10
data_cv <- create_folder(k, dat_train)
vec_remove <- c("group")

# Classifieur 1 - Arbre classique 
library(rpart)

# Arbre 1 - Gini
err_tree_1 = rep(NA, k)
err_tree_1 <- sapply(1:k, function(x){
  model = rpart(diabetes ~ ., data = data_cv[group != x, -(vec_remove), with = F])
  out = predict(model, data_cv[group == x, -(vec_remove), with = F], type = "class")
  sum(out != data_cv[group == x, diabetes]) / nrow(data_cv[group == x,])
})

# Erreur test
model = rpart(diabetes ~ ., data = dat_train)
out = predict(model, dat_test, type = "class")
err_tree_1_test <- sum(out != dat_test[, diabetes]) / nrow(dat_test)



# Arbre 2 - Cross-entropy
err_tree_2 = rep(NA, k)
err_tree_2 <- sapply(1:k, function(x){
  model = rpart(diabetes ~ ., data = data_cv[group != x, -(vec_remove), with = F], parms = list(split = "information"))
  out = predict(model, data_cv[group == x, -(vec_remove), with = F], type = "class")
  sum(out != data_cv[group == x, diabetes]) / nrow(data_cv[group == x,])
})

# Erreur test
model = rpart(diabetes ~ ., data = dat_train, parms = list(split = "information"))
out = predict(model, dat_test, type = "class")
err_tree_2_test <- sum(out != dat_test[, diabetes]) / nrow(dat_test)


plot <- list()
for (x in 1:k) {
  model = rpart(diabetes ~ ., data = data_cv[group != x, -(vec_remove), with = F], parms = list(split = "information"))
  pdf(NULL)
  dev.control(displaylist="enable")
  plotcp(model)
  plot[[x]] <- recordPlot()
  invisible(dev.off())
}

# Arbre 3 - Cross-entropy pruned
err_tree_3 = rep(NA, k)
err_tree_3 <- sapply(1:k, function(x){
  model = rpart(diabetes ~ ., data = data_cv[group != x, -(vec_remove), with = F], parms = list(split = "information"))
  model = prune.rpart(model, cp = 0.0)
  out = predict(model, data_cv[group == x, -(vec_remove), with = F], type = "class")
  sum(out != data_cv[group == x, diabetes]) / nrow(data_cv[group == x,])
})

# Erreur test
model = rpart(diabetes ~ ., data = dat_train, parms = list(split = "information"))
model = prune.rpart(model, cp = 0.1)
out = predict(model, dat_test, type = "class")
err_tree_3_test <- sum(out != dat_test[, diabetes]) / nrow(dat_test)


# Classifieur 2 - Arbres bagging
library(randomForest)

# Impute the data very roughly
data_cv_imputed <- na.roughfix(data_cv)
data_test_imputed <- na.roughfix(dat_test)
data_train_imputed <- na.roughfix(dat_train)

err_bag = rep(NA, k)
err_bag <- sapply(1:k, function(x){
  model = randomForest(diabetes ~ ., data = data_cv_imputed[group != x, -(vec_remove), with = F], mtry = ncol(data_cv) - 1)
  out = predict(model, data_cv_imputed[group == x, -(vec_remove), with = F])
  sum(out != data_cv_imputed[group == x, diabetes]) / nrow(data_cv_imputed[group == x,])
})

# Test error
model = randomForest(diabetes ~ ., data = data_train_imputed, mtry = (ncol(data_train_imputed) - 1))
out = predict(model, data_test_imputed)
err_bag_test <- sum(out != data_test_imputed[, diabetes]) / nrow(data_test_imputed)

# Classifieur 3 - Random forrest
var_imp <- list()
err_rf = rep(NA, k)
err_rf <- sapply(1:k, function(x){
  model = randomForest(diabetes ~ ., data_cv_imputed[group != x, -(vec_remove), with = F], mtry = 4, ntree = 1000)
  var_imp[[x]] <<- importance(model)
  out = predict(model, data_cv_imputed[group == x, -(vec_remove), with = F])
  sum(out != data_cv_imputed[group == x, diabetes]) / nrow(data_cv_imputed[group == x,])
})


# Test error
model = randomForest(diabetes ~ ., data = data_train_imputed, mtry = 4, ntree = 500)
out = predict(model, data_test_imputed)
err_rf_test <- sum(out != data_test_imputed[, diabetes]) / nrow(data_test_imputed)

# Classifieur 4 - LDA
library(MASS)
err_lda = rep(NA, k)
err_lda <- sapply(1:k, function(x){
  model = lda(diabetes ~ ., data = data_cv_imputed[group != x, -(vec_remove), with = F])
  out = predict(model, data_cv_imputed[group == x, -(vec_remove), with = F])$class
  sum(out != data_cv_imputed[group == x, diabetes]) / nrow(data_cv_imputed[group == x,])
})

# Test error
model = lda(diabetes ~ ., data = data_train_imputed)
out = predict(model, data_test_imputed)$class
err_lda_test <- sum(out != data_test_imputed[, diabetes]) / nrow(data_test_imputed)

# Classifieur 5 - QDA
err_qda = rep(NA, k)
err_qda <- sapply(1:k, function(x){
  model = qda(diabetes ~ ., data = data_cv_imputed[group != x, -(vec_remove), with = F])
  out = predict(model, data_cv_imputed[group == x, -(vec_remove), with = F])$class
  sum(out != data_cv_imputed[group == x, diabetes]) / nrow(data_cv_imputed[group == x,])
})

# Test error
model = qda(diabetes ~ ., data = data_train_imputed)
out = predict(model, data_test_imputed)$class
err_qda_test <- sum(out != data_test_imputed[, diabetes]) / nrow(data_test_imputed)

# Classieur 6 - Logistic regression
err_log = rep(NA, k)
summary <- list()
thres <- 0.4 # choisi par CV
err_log <- sapply(1:k, function(x){
  model = glm(diabetes ~ ., data = data_cv_imputed[group != x, -(c(vec_remove, "triceps", "insulin", "pressure", "age")), with = F], family = "binomial")
  summary[[x]] <<- summary(model)
  out = predict(model, data_cv_imputed[group == x, -(c(vec_remove, "triceps", "insulin", "pressure", "age")), with = F])
  out2 = ifelse(out >= thres, "pos", "neg")
  sum(out2 != data_cv_imputed[group == x, diabetes]) / nrow(data_cv_imputed[group == x,])
})

# Test error
model = glm(diabetes ~ ., data = data_train_imputed[, -(c(vec_remove, "triceps", "insulin", "pressure", "age")), with = F], family = "binomial")
out = predict(model, data_test_imputed[, -(c(vec_remove, "triceps", "insulin", "pressure", "age")), with = F])
out2 = ifelse(out >= thres, "pos", "neg")
err_log_test <- sum(out2 != data_test_imputed[, diabetes]) / nrow(data_test_imputed)


# Classifieur 7- Naive Bayes
library(e1071)

err_bayes = rep(NA, k)
err_bayes <- sapply(1:k, function(x){
  model = naiveBayes(diabetes ~ ., data = data_cv_imputed[group != x, -(vec_remove), with = F])
  out = predict(model, data_cv_imputed[group == x, -(vec_remove), with = F])
  sum(out != data_cv_imputed[group == x, diabetes]) / nrow(data_cv_imputed[group == x,])
})

# Test error
model = naiveBayes(diabetes ~ ., data = data_train_imputed)
out = predict(model, data_test_imputed)
err_bayes_test <- sum(out != data_test_imputed[, diabetes]) / nrow(data_test_imputed)

```



### Testez vos modèles finaux sur le jeu de données de validation

Il s'agit ici de tester tous les modèles ajustés précédemment sur le jeu de données de validation avec 50 observations. 

* Quel est le meilleur modèle?
* Que pouvez-vous dire de la précision de votre estimation de l'erreur de classification obtenue avec le jeu de données de validation?

```{r resultats, echo=FALSE}
vec_modele <- c("Arbre Gini", "Arbre cross-entropy", "Arbre cross-entropy pruned", "Bagging", "Random forrest", "LDA", "QDA", "Reg. logistique", "Naive Bayes")
vec_validation <- c(mean(err_tree_1), mean(err_tree_2), mean(err_tree_3), mean(err_bag), mean(err_rf), mean(err_lda), mean(err_qda), mean(err_log), mean(err_bayes))
vec_test <- c(err_tree_1_test, err_tree_2_test, err_tree_3_test, err_bag_test, err_rf_test, err_lda_test, err_qda_test, err_log_test, err_bayes_test)

data_resultats <- data.table(
  modele = vec_modele,
  erreur_validation = vec_validation,
  erreur_test = vec_test
)[order(erreur_test)]

knitr::kable(data_resultats, digits = 3, row.names = FALSE, col.names = c("Modèle", "Erreur moyenne de validation", "Erreur test"), align = rep("c", 3), caption = "Performances des différents classifieurs sur les jeux de données CV et test")
```