---
title: "Projet étude de régression et sélection de gènes pour le cancer de la vessie"
output: html_document
---

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

```


```{r}
library('pls')
library('curatedBladderData')
library('Biobase')
library('glmnet')
data(package="curatedBladderData")

```

```{r}
data(GSE89_eset)
expressionData <- exprs(GSE89_eset)
```
Préparation de données :


```{r}
X <- t(expressionData)
otherData <- pData(GSE89_eset@phenoData)
Y <- (otherData$summarystage == "invasive")*1
```


```{r}
Data = data.frame(X = I(t(expressionData)), Y = Y)
train <- rbinom(length(Data$Y),1,0.8)
Data.train <- c()
Data.test <- c()
Data.train$X <- Data$X[train==1,]
Data.train$Y <- Data$Y[train==1]
Data.test$X <- Data$X[train==0,]
Data.test$Y <- Data$Y[train==0]
Data.train <- as.data.frame(Data.train)
Data.test <- as.data.frame(Data.test)


```

```{r}
gene_names <- colnames(Data$X)
```

2. Fonction Bootstrap qui génère un échantillon bootstrap à partir du jeu initial.

```{r}
bootstrap <- function(Data) {
  indice <- sample(1:nrow(Data) ,nrow(Data) ,replace = TRUE)
  return(Data[indice, ])
}

```

On note n le nombre de valeur que l'on va retenir.
```{r}
fit_function <- function(Data){
  pcrbladder <- pcr(Data$Y ~ Data$X , validation="CV")
  ncomponent <- pcrbladder$ncomp
  reduction.matrix <- pcrbladder$loadings[,1:ncomponent]
  Data$reducedX <- Data$X%*%reduction.matrix
  pcr.model <- glm(Data$Y~Data$reducedX ,data = Data
                   ,family =
                  binomial(link="logit"))
  matricecoeff <- pcr.model$coefficients
  return(matricecoeff)
}
```

###Méthode PCR : ###
```{r}
Selection_Pcr <- function(Data, n){
  pcrbladder <- pcr(Data$Y ~ Data$X , validation="CV")
  ncomponent <- pcrbladder$ncomp
  reduction.matrix <- pcrbladder$loadings[,1:ncomponent]
  Data$reducedX <- Data$X%*%reduction.matrix
  pcr.model <- glm(Data$Y~Data$reducedX ,data = Data
                   ,family =
                     binomial(link="logit"))
  matricecoeff <- pcr.model$coefficients
  matricecoeff[is.na(matricecoeff)] <- 0
  coeff_finaux <- reduction.matrix%*%matricecoeff[seq(2,ncomponent+1)]
  result <- order(abs(coeff_finaux) ,decreasing
                  = TRUE)
  result <- result[seq(1, n)]
  gene <- colnames(pcrbladder$model$`Data$X`)
  Data$reducedX <- NULL
  gene_selected <- gene[result]
  return(gene_selected)
}

```



### Méthode PLS : ###

```{r}
Selection_Pls <- function(Data, n){
  plsbladder <- plsr(Data$Y ~ Data$X ,
                     validation="CV")
  ncomponent <- plsbladder$ncomp
  reduction.matrix <- plsbladder$loadings[,1:ncomponent]
  Data$reducedXpls <- Data$X%*%reduction.matrix
  plsr.model <- glm(Data$Y~Data$reducedXpls ,data = Data
                    ,family =
                     binomial(link="logit"))
  matricecoeff <- plsr.model$coefficients
  matricecoeff[is.na(matricecoeff)] <- 0
  coeff_finaux <- reduction.matrix%*%matricecoeff[seq(2, ncomponent+1)]
  result <- order(abs(coeff_finaux) ,decreasing
                  = TRUE)
  result <- result[seq(1, n)]
  gene <- colnames(plsbladder$model$`Data$X`)
  Data$reducedX <- NULL
  gene_selected <- gene[result]
  return(gene_selected)
}


```


Un nombre B fois :
<ol>
<li> on génére un jeu de données bootstrap <li>
<li> on applique la fonction de la question précédente <li>
On retient pour chaque gène un score d'apparition.

On retourne une liste qui contient tous les gènes retenue + les scores des n gènes séléctionner.
écrire une fonction teste qui vérifie si un gène appartient à mon data frame ou pas.

En appliquant GLM après le bootstrap, on observe l'apparition de valeur NA dans notre matrice de coeffiscient, Le modèle est inutilisable. Certains coefficients n’ont pas été estimés et correspondent à la valeur NA (not available). Nous ne pouvons ni les interpréter, ni les déployer sur des
individus supplémentaires néanmoins on comprends que cela veut dire que les directions de valeurs associées NA sont  colinèaire à une autre direction, elle n'apporte pas de nouvelle information. On choisit de leurs associer un coeffiscient nul.
Puisqu'elle n'apporte aucunes informations et de procèder. Ainsi, on tue l'information apporté par l'axe puisque dejà présente dans un autre
Néanmoins, on ne comprend pas, pourquoi cette direction a été retenu par l'acp...


Gènes retenue par la méthode pls et leurs scores respectifs :

```{r}
B <- 20
test_pls <- NULL
for(i in seq(1:B)){
  Data_boots <- bootstrap(Data)
  test_pls <- c(test_pls,
          Selection_Pls(Data_boots,50))
}
result_pls <- table(test_pls)/B
gene_pls <- result_pls[result_pls > 0.85]
sort(gene_pls, decreasing = TRUE)
```

Gènes retenue par la méthode pcr et leurs scores respectifs :

```{r}
B <- 20
test_pcr <- NULL
for(i in seq(1:B)){
  Data_boots <- bootstrap(Data)
  test_pcr <- c(test_pcr,
            Selection_Pcr(Data_boots,50))
}
result_pcr <- table(test_pcr)/B
gene_pcr <- result_pcr[result_pcr > 0.95]
sort(gene_pcr, decreasing = TRUE)
```



Certains gènes vous semblent-ils candidats pour
faire partie d’une signature à valider biologiquement? Proposer une manière
de sélectionner les candidats potentiels.

Les gènes qui ont été sélectionnées par les deux méthodes, pourraient éventuellement être intéressant à étudier.


5. Que peut-on reprocher à cette méthode en terme d’interprétation?
En termes de prise en compte de la colinéarité?

On peut repprocher à cette méthode son manque d'interprétabilité, comment fixer le seuil. si X_1 et X_2 sont colinéaire ben il y en a qu'une seule qui va apparaitre dans la selection de composante. Comment privilégier l'un ou l'autre.
Certains gènes sont toujours sélectionner, malgré le bootstrap
QUESTION :
Que peut-on de plus reprocher à l’utilisation de la PLS dans le cas d’une variable à expliquer binaire ?

ça reste très heurestique comme approche, on détruit les directions colinéraires dans glm. On fait pas mal d'hypothèses et appriori


\textbf{Question 6}

L’introduction de pénalisations en norme L1 induit une sélection de variables optimale en régression Cette démarche conduit à la définition de versions parcimonieuses de l’Analyse en Composantes Principales et de la régression PLS pour différents objectifs : exploration, comparaison ou intégration de deux jeux de données en version régression ou canonique, analyse discriminante PLS.



### 5. Sans réduction de dimension préalable:###

_La librairie "glmnet"_ est connue pour son efficacité (rapidité
des calculs, qualité des résultats).


7. Une troisième manière de faire est de s'affranchir de la réduction de dimension en utilisant une régression pénalisée en utilisant la fonction _glmnet_ du package _glmnet_
Pour avoir des résultats parcimonieux.

Lors de l'appel de la fonction, on précise la valeur 1 pour alpha  pour avoir une régression de Lasso ce qui permettra de fixer selectivement les coeffiscient à la valeur 0 et réaliser une selection de variable et nous laissons R choisir la valeur adéquat pour lambda

```{r}
fit <- glmnet(Data.train$X, Data.train$Y,family = "binomial", standardize = FALSE, alpha = 1)
```

```{r}
fit
```

```{r}
plot(fit, xvar = "lambda")
```

Identification de la valeur « optimale » de λ. Nous avons les logarithmes des λi en abscisse
de notre graphique. Ils varient de log(0.008982)=-4.7124 à log(89.825)=4.4978. Une règle
empirique consiste à choisir la valeur de λ à partir de laquelle les coefficients commencent
« à se stabiliser ».
La solution viable consiste à passer par la validation croisée,
qui ne fait intervenir que l’échantillon d’apprentissage. Nous faisons appel à la fonction
cv.glmnet()

Fonction sélection de gène par méthode glmnet :

```{r}
Selection_glmnet <- function(Data){
  model.glmnet <- cv.glmnet(Data$X, Data$Y,family = "binomial", type.measure="class",
 nfolds=10,alpha=1,keep=TRUE)
  lambda_min <- model.glmnet$lambda.min
  Coeff_sparse <- coef(model.glmnet, s = lambda_min)
  result_l <- order(abs(Coeff_sparse), decreasing = TRUE)
  l <- length(Coeff_sparse[Coeff_sparse >0])
  resultt_l <- result_l[seq(2,l)] # On récupère les coeffiscient non nulles de la matrice
  gene <- colnames(Data$X)
  gene_selection <- gene[resultt_l]
  return(gene_selection)
}
```


Selection de gènes par méthode bootstrap :

```{r results='hide', message=FALSE, warning=FALSE}
B <- 20
test_glmnet <- NULL
result_glmnet <- NULL
gene_glmnet <- NULL

for(i in seq(1:B)){
  Data_boots <- bootstrap(Data)
  test_glmnet <- c(test_glmnet,
          Selection_glmnet(Data_boots))
}

result_glmnet <- table(test_glmnet)/B
quantile_glmnet <- quantile(result_glmnet)[4]
gene_glmnet <- result_glmnet[result_glmnet > quantile_glmnet]
gene_glmnet


```