BARNABE_DAVY.Rmd

---
title: "Projet AnaStat"
output: html_document
---

## Chargement des données et construction des jeux d'apprentissage et de test

Les données représentent 184 femmes qui ont été touchées par le cancer du sein ainsi que leur taux d'expression pour 4654 gènes. Nous disposons également d'une variable "label" indiquant la rechute ou non pour chaque femme. L'objectif ici est de déterminer la meilleure méthode d'apprentissage statistique pour prédire la rechute du cancer du sein. Tout d'abord les données sont chargées et mises en forme. Nous avons testé les méthodes avec un jeu de données d'apprentissage et un jeu test. Au vu du faible nombre d'observations, nous avons réalisé une boucle pour chaque méthode avec 5 seed différents puis moyenné les erreurs de classification afin de comparer plus justement les méthodes entre elles. 
Nous avons donc 184 observations pour 4654 covariables et une variable à expliquer. Nous avons donc beaucoup plus de paramètres que d'observations ce qui complique la construction d'un bon modèle. 

```{r}
# setwd("~/Documents/M2-AMI2B/Apprentissage_Statistique/projet/")
train <- t(read.table("xtrain.txt", sep="\t", as.is = NA))
data = numeric()
for (i in 2:nrow(train)){
  data <- rbind(data, as.numeric(train[i,]))
}
data <- cbind(data, read.table("ytrain.txt"))
colnames(data) <- c(train[1,], "label")
data <- within(data, {
  label = factor(label, labels = c("-1","1"))
})
cat("Dimension des données :",dim(data), "\n")
str(data)
table(data$label)
```

## Régression logistique multiple

La méthode de régression logistique multiple vise à sélectionner les paramètres en testant l'influence de chaque covariable sur le modèle. Si l'absence de la covariable réduit la qualité du modèle alors sa présence est importante. Ainsi, le meilleur modèle est sélectionné. 
```{r}
data.glm <- glm(label~., data=data, family = binomial)
summary(data.glm)
```

On observe que l'algorithme ne converge pas sur les données surement du au fait du grand nombre de covariables (4654) par rapport au nombre d'observations (184). En effet, la méthode glm ne permet d'estimer que n-1 paramètres où n est le nombre d'observations. Ce modèle ne correspond donc pas aux données. 

## Arbre CART

La méthode CART vise à construire un arbre de décision à partir des covariables. A chaque noeud construit est associé une condition dépendant des covariables. On construit dans un premier temps un arbre maximal à partir des données d'apprentissages avec la fonction rpart puis on réalise une étape d'élagage pour réduire le nombre de feuilles. L'élagage optimal de l'arbre est réalisé en choisissant la valeur de cp pour laquelle l'erreur de validation croisée du rpart est minimale.
```{r}
error.cart <- numeric() ; error.cartprune <- numeric()
library(rpart)
for(i in c(2018,1234,2017,1996,1107)) {
  set.seed(i)
  
  test <- sample(1:nrow(data),round(nrow(data)/3))
  train <- (1:nrow(data))[-test]
  train <- data[train,]
  test <- data[test,]
  t.max = rpart(label~., data=train, control=rpart.control(cp = 0))
  
  # Elagage de l'arbre en choisissant la valeur de cp pour laquelle
  # erreur de validation croisée du rpart est minimale
  t.prune = prune(t.max, cp = t.max$cptable[which.min(t.max$cptable[,4]), 1])
  
  # cat("Arbre maximal\n")
  # print(t.max)
  # plot(t.max)
  # text(t.max, cex=1)
  # 
  # cat("\nArbre après élagage\n")
  # print(t.prune)
  # plot(t.prune)
  # text(t.prune, cex=1)
  
  pred_t = predict(t.max, newdata = test, type="class")
  pred_tprune = predict(t.prune, newdata = test, type="class")
  
  error.cart <- c(error.cart ,sum(test$label != pred_t)/length(test$label))
  error.cartprune <- c(error.cartprune, sum(test$label != pred_tprune)/length(test$label) )
}
cat("Erreur de classification moyenne pour l'arbre maximal :", mean(error.cart))
cat("\nErreur de classification moyenne pour l'arbre élagué :", mean(error.cartprune))
```

On observe une erreur d'apprentissage d'en moyenne 47% pour l'arbre maximal et 36,7% pour l'arbre élagué avec les seed que nous avons utilisé. 
L'erreur importante et les inconvénients de l'arbre CART nous ont poussé à étudier d'autres méthodes. En effet, la construction de l'arbre dépend du découpage des données et n'est pas constante. Cette méthode n'est donc pas robuste. Une méthode se basant également sur les arbres permet de contourner ce problème :  il s'agit des forêts aléatoires.

## Forêt aléatoire

La méthode de forêt aléatoire est une méthode de classification adaptée aux gros jeux de données. Elle permet de construire plusieurs arbres CART maximals et les moyenner afin de diminuer la variance élévée d'un arbre seul. On observe ici un taux d'erreur de 39,34% avec les paramètres par défaut. 

```{r}
error.rf <- numeric()
library(glmnet)
library(randomForest)
for(i in c(2018,1234,2017,1996,1107)) {
  set.seed(i)
  
  test <- sample(1:nrow(data),round(nrow(data)/3))
  train <- (1:nrow(data))[-test]
  train <- data[train,]
  test <- data[test,]
  
  data.rf = randomForest(label~., data=train)
  
  pred.rf <- predict(data.rf, newdata=test)
  error.rf <- c(error.rf, mean(test$label != pred.rf))
}
mean(error.rf)
```


On cherche à optimiser les paramètres du randomForest. Pour cela on va chercher le nombre optimal de nombre d'arbres dans la forêt visuellement, à l'aide d'un graphe, en observant l'erreur estimée (OOB). On observe que l'erreur se stabilise à partir 1000 arbres. On choisit donc cette taille pour la forêt afin de minimiser l'erreur. Néanmoins, il ne faut pas trop augmenter ntree pour ne pas réaugmenter l'erreur de classification. On utilise ensuite la fonction tuneRF qui permet de trouver la meilleure valeur pour mtry, le nombre de variables sélectionnées à chaque noeud. On lance tuneRF avec un nombre d'arbre de 2000. Le mtry avec la plus petite erreur est 136. En utilisant ces nouveaux paramètres dans le randomForest, l'erreur de classification passe de 39,34% à 36,7% pour les mêmes jeux de données. 
```{r}
library(randomForest)
test <- sample(1:nrow(data),round(nrow(data)/3))
train <- (1:nrow(data))[-test]
train <- data[train,]
test <- data[test,]
test.rf = randomForest(label~., data=train, mtry=136, ntree=3000)
plot(1:test.rf$ntree, test.rf$err.rate[,1], type='l')
    
Y <- train$label
tunerf <- tuneRF(train[,-ncol(train)],Y, ntreeTry = 2000)
tunerf
test.pred.rf <- predict(test.rf, newdata=test)
mean(test$label != test.pred.rf)

error.rf <- numeric()
for(i in c(2018,1234,2017,1996,1107)) {
  set.seed(i)
  
  test <- sample(1:nrow(data),round(nrow(data)/3))
  train <- (1:nrow(data))[-test]
  train <- data[train,]
  test <- data[test,]
  data.rf
  
  pred.rf <- predict(data.rf, newdata=test)
  error.rf <- c(error.rf, mean(test$label != pred.rf))
  
  # Courbe ROC
  pred.rf = predict(data.rf, newdata = test, type = "prob")
  
  s = seq(0,1,0.01)
  
  absc = numeric(length(s))
  ordo = numeric(length(s))
  
  for (i in 1:length(s)) 
  {
    ordo[i] = sum(pred.rf[,"1"] >= s[i] & test$label == "1") / sum(test$label == "1")  
    absc[i] = sum(pred.rf[,"1"] >= s[i] & test$label == "-1") / sum(test$label == "-1")
  }
  
  lines(absc,ordo,type = "l",col = i)
}
mean(error.rf)
```


## Méthode de bagging

Il s'agit d'une méthode similaire au random forest mais on sélectionne toutes les covariables pour la construction de chaque noeud. Cette méthode se base sur le boostrap (ré-échantillonnage) pour générer ensuite un arbre maximal. Cette étape est répétée jusqu'à atteindre la condition d'arrêt où les arbres seront moyennés. Cette méthode est améliore souvent la qualité de la prédiction. 
```{r}
error.bag <- numeric()
for(i in c(2018,1234,2017,1996,1107)) {
  set.seed(i)
  
  test <- sample(1:nrow(data),round(nrow(data)/3))
  train <- (1:nrow(data))[-test]
  train <- data[train,]
  test <- data[test,]
  
  fit.bag <- randomForest(label~., data = train, mtry=ncol(data)-1, ntree=2000)
  fit.bag
  
  pred.bag <- predict(fit.bag, newdata=test)
  error.bag <- c(error.bag, mean(test$label != pred.bag))
}
mean(error.bag)
```
Ici, on peut voir que l'erreur de classification sur le bagging est de 38,4% donc supérieure au randomForest. Nous allons donc utiliser une nouvelle méthode permettant de réduire le modèle afin d'utiliser le glm. 

## Méthode LASSO

LASSO est une méthode de sélection de variables permettant de réduire les modèles avec de nombreuses covariables. Elle va tester la corrélation des covaribales et en choisir une comme représentative pour un groupe très proche. Cependant, la variable choisit est aléatoire donc limite l'interprétation. En effet, dans notre cas, un gène sélectionné comme représentatif pourrait ne pas réellement agir sur la rechute du cancer du sein mais simplement représenté un groupe de covariables corrélées contenant la vraie variable explicative. 
Cependant, cette méthode nous a permis de sélectionner une vingtaine de covariables et donc de lancer un glm. L'optimisation de LASSO, donc du paramètre lambda, a été faite visuellement en utilisant le graphique représentant le modèle. On observe une courbe diminuant jusqu'à arriver à un plateau. Le nombre optimal de variables dans le modèle se situe au niveau du coude entre la pente et le plateau. On fixe donc la valeur de lambda correspondant au coude. 
On peut ensuite récupérer les variables sélectionnées et les utiliser dans le glm. On obtient alors une erreur de classification de 40%. 
```{r}
error.lasso <- numeric() ; nbvar <- c()
error.rf.la <- numeric()
error.svm.la <- numeric()
library(glmnet)
library(randomForest)
library(e1071)
for(i in c(2018,1234,2017,1996,1107)) {
  set.seed(i)
  
  test <- sample(1:nrow(data),round(nrow(data)/3))
  train <- (1:nrow(data))[-test]
  train <- data[train,]
  test <- data[test,]
  
  X <- model.matrix(label~., data)[,-(ncol(data))]
  Y <- data$label
  
  fit.lasso = glmnet(X,Y,family = "binomial")
  cv.out = cv.glmnet(X,Y, family = "binomial", lambda = seq(exp(-6),1,length = 500))
  cv.out.zoom = cv.glmnet(X,Y, family = "binomial", lambda = seq(exp(-6),exp(-2.5),length = 500))
  
  # par(mfrow=c(1,2))
  # plot(cv.out)
  # plot(cv.out.zoom)
  
  min = cv.out.zoom$lambda.min
  
  Dpred.lasso = predict(fit.lasso, type="coefficients", s=min)
  beta = predict(fit.lasso, type="nonzero", s=min)   # type 'nonzero' retourne les indices des coefficients non nuls
  
  nbvar <- c(nbvar, length(beta[,])) # La longueur de beta correspond donc a la dimension du modèle séléctionné
  
  ### GLM
  select = which(coef(cv.out.zoom)!=0)
  data.glm <- glm(label~., data=train[,c(select,ncol(train))], family = binomial)
  summary(data.glm)
  
  proba.glm = predict(data.glm, newdata=test[,c(select,ncol(test))], type="response")
  pred.glm <- rep("-1", length(test$label))
  pred.glm[proba.glm > 1/2] = "1"
  
  error.lasso <- c(error.lasso, mean(test$label != pred.glm))

}
```
On observe une erreur de 40% sur le glm. 


## Méthode SVM

La méthode SVM (Support Vector Machine) a pour objectif de définir un hyperplan dans l'espace des covariables afin séparer les valeurs de la variable àexpliquer en deux sous-ensemble.
La fonction tune permet de tester les différents couples de valeur des paramètres "cost" et "gamma" et de trouver le meilleur couple de valeur pour la prédiction.
```{r}
library(e1071)
error.svm <- numeric()
for(i in c(2018,1234,2017,1996,1107)) {
  set.seed(i)
  
  test <- sample(1:nrow(data),round(nrow(data)/3))
  train <- (1:nrow(data))[-test]
  train <- data[train,]
  test <- data[test,]
  tune.out = tune(svm, label~., data = train, kernel = "radial",
                  ranges = list(
                    cost = c(1e-15,1e-7,1,10),
                    gamma = c(1e-15,1e-7,0.1,1),
                    # le nombre d'observation est faible
                    # nous avons donc choisi de réduire le nombre de partion
                    # pour la cross validation du tune
                    tunecontrol = tune.control(cross = 3)
                  ))
  summary(tune.out)
  
  fit.svm = svm(label~., train, kernel="radial", probability=TRUE, cost=tune.out$best.parameters$cost, gamma=tune.out$best.parameters$gamma)
  
  cat("Erreur de classifaction :",mean(test$label != predict(fit.svm,newdata = test)))
  
  error.svm = c(error.svm, mean(test$label != predict(fit.svm,newdata = test)))
}
mean(error.svm)
```

# Prediction du jeu de donnée test
## Chargement des données
```{r}
donnee_test <- t(read.table("xtest.txt", sep="\t", as.is = NA))
data_test = numeric()
for (i in 2:nrow(donnee_test)){
  data_test <- rbind(data_test, as.numeric(donnee_test[i,]))
}
colnames(data_test) <- donnee_test[1,]
cat("Dimension des données :",dim(data_test), "\n")
```

## Construction du modele avec la methode random forest
Pour la construction du modèle, les valeurs optimisées des parametres mtry et ntree sont utilisées
```{r}
fit <-  randomForest(label~., data=data, mtry=136, ntree=2000)

pred <- predict(fit, newdata = data_test, type = "prob")

vecteur <- pred [,2]

write(x = vecteur, file = "ytest.txt", sep = "\n")
```