projet_stat.Rmd

---
title: "projet_stat"
output:
  word_document: default
date: "2023-01-18"
---

```{r, echo=FALSE}
#setwd("~/Documents/3A/Stats/projet") #mettre votre directory
```

Dans les données « nutrimouse », nous avons les mesures de l'expression de 120 gènes potentiellement impliqués dans les problèmes nutritionnels des souris, les concentrations de 21 acides gras hépatiques et les génotypes de 40 souris soumises à 5 différents régimes.

"PPAR" est une abréviation pour "Peroxisome Proliferator-Activated Receptor", qui est un récepteur dans les cellules qui régule divers processus métaboliques, tels que la décomposition des graisses et la production d'énergie. 


## Librairies
```{r setup, include=FALSE}
library(corrplot)
library(FactoMineR)
library(factoextra)
library(glmnet)
library(MultiVarSel)
```

## Dataframe

```{r}
load('nutrimouse.rda')
```

## Pre-processing

```{r}
data_gene = nutrimouse[['gene']]
data_lipid = nutrimouse[['lipid']]
diet = nutrimouse[['diet']]
genotype = nutrimouse[['genotype']]
data_both = cbind(data_gene, data_lipid) # les gènes et les lipides 
data_both_gen = cbind(diet, data_both) # les dietes + les genes + les lipides 
data = cbind(genotype, data_both_gen) # toutes les donnés 
#rm(data_gene, data_lipid, diet, genotype, data_both, data_both_gen, nutrimouse)
```

```{r}

#summary(data)
#summary(data$diet)
#summary(data$genotype)

```

```{r}
X1 <- data[, 1]
X2 <- data[, 2]
table(X1, X2)
```

## Matrices de corrélation sur l'ensemble des individus

```{r}
data_gen <- data[, 1:122]
corrplot(cor(data_gen[, 3:122]))
```

```{r}
data_lipid <- cbind(data[,1:2], data[, 123:143])
corrplot(cor(data_lipid[, 3:23]))
```

## Matrices de corrélation sur chaque groupe de souris (wt ou ppar)

```{r}
data_lipid_wt <- cbind(data[1:20, 1:2], data[1:20, 123:143])
corrplot(cor(data_lipid_wt[,3:23]))
```

```{r}
data_lipid_ppar <- cbind(data[21:40,1:2], data[21:40, 123:143])
corrplot(cor(data_lipid_ppar[3:23]))
```

## PCA sur les gènes 


```{r}
res.pca <- PCA(nutrimouse$gene, graph = FALSE) #TRUE affiche 2 graphes directement
get_eigenvalue(res.pca) #donne un tableau avec les valeurs propres
fviz_pca_var(res.pca, axes=1:2)
nutrimouse$genotype <- as.factor(nutrimouse$genotype )
#data_gene$genotype <- as.factor(data_gene$genotype)

colors <- c("red", "blue")
fviz_pca_ind(res.pca, axes = 1:2, geom = "point",
             habillage = nutrimouse$genotype,
             col.hab = colors,
             ggtheme = theme_classic(),
             legend = "bottom")

fviz_eig(res.pca, addlabels = TRUE, ylim = c(0,50)) #graphe avec inertie
```

## PCA sur les lipides 

```{r}
res.pca <- PCA(nutrimouse$lipid, graph = FALSE) #TRUE affiche 2 graphes directement
get_eigenvalue(res.pca) #donne un tableau avec les valeurs propres
fviz_pca_var(res.pca, axes=1:2)
#nutrimouse$genotype <- as.factor(nutrimouse$genotype )
#data_gene$genotype <- as.factor(data_gene$genotype)

colors <- c("red", "blue")
fviz_pca_ind(res.pca, axes = 1:2, geom = "point",
             habillage = nutrimouse$diet,
             col.hab = colors,
             ggtheme = theme_classic(),
             legend = "bottom")

fviz_eig(res.pca, addlabels = TRUE, ylim = c(0,50)) #graphe avec inertie
```
```{r}
#install.packages('plotly', dependencies = TRUE)
#install.packages('gridExtra', dependencies = TRUE)
#library(plotly)
library(gridExtra)
library(ggplot2)

par(mfrow = c(2, 2))

fig1 <- boxplot(nutrimouse$lipid[,3] ~ nutrimouse$diet, xlab = "régime", ylab = "lipide n° 3")
fig2 <- boxplot(nutrimouse$lipid[,13] ~ nutrimouse$diet, xlab = "régime", ylab = "lipide n° 13")
fig3 <- boxplot(nutrimouse$lipid[,3] ~ nutrimouse$genotype, xlab = "génotype", ylab = "lipide n° 3")
fig4 <- boxplot(nutrimouse$lipid[,13] ~ nutrimouse$genotype, xlab = "génotype", ylab = "lipide n° 13")
title("Influence du régime ou du génotype sur l'expression lipidique", outer = TRUE, line = -0.9)
#ggsave("boxplot.pdf")



```
```{r}
heatmap(abs(cor(nutrimouse$lipid)), symm = TRUE)
ggsave("heatmap_lipides.pdf")
```


```{r}
boxplot(data_gen[,4]~data_gen[,2], col = "green")
```

```{r}
heatmap(abs(cor(data_gene[, -c(1,2)])), symm = TRUE)
```

```{r}
heatmap(abs(cor(data[, -c(1,2)])), symm = TRUE)
```

```{r}
ref_data <- subset(data, diet == "ref")
```

## Selection de variables 

### Les acides gras selon le génotype 

Soit Y la matrice des acides gras
En ligne les individus en colonnes les differents acides gras.
Soit X la matrice des variables explicatives : les gènes. 

On suppose une indépendance des lignes sur les Ei.

```{r}
#Y <- as.matrix(data_lipid[, -c(1, 2)])
Y <- as.matrix(nutrimouse$lipid)
X1 <- as.matrix(nutrimouse$genotype)
X2 <- as.matrix(nutrimouse$diet)
#X1 <- as.matrix(data[1])
#X2 <- as.matrix(data[2])
dim(Y)
dim(X1)
dim(X2)
```

```{r}
X <- model.matrix(lm(Y~X1+0)) # on crée la matrice de design
p <- ncol(X) #2
n <- nrow(X) #40
q <- dim(Y)[2] #21
```

```{r}
Y <- scale(Y)
Y
```

Sélection des variables

```{r}
residuals <- lm(Y~X-1)$residuals
pvalue = whitening_test(residuals)
pvalue
```

P-value très faible, on rejette H0: nos colonnes sont dépendantes.
H0 = colonnes indépendantes,
H1 = colonnes dépendantes.

Nous ne pouvons pas vectoriser notre Y tel quel, il faut rendre indépendantes nos colonnes.

Nous allons blanchir nos données. Il faut choisir la méthode de blanchiement la plus adaptée.

```{r}
result=whitening_choice(residuals,c("AR1","nonparam","ARMA"),pAR=1,qMA=1)
result
```

Il faut connaître les différents principes associés à ces méthodes. Nous choisissons la méthode non-paramétrique.

```{r}
square_root_inv_hat_Sigma=whitening(residuals,"nonparam",pAR=1,qMA=0)
square_root_inv_hat_Sigma
```

--> Matrice par laquelle nous multiplions notre formule Y = XB + E afin de rendre indépendantes nos colonnes.

```{r}
# Blanchiement / vectorisation / LASSO / Stability selection 
Frequencies=variable_selection(Y,X,square_root_inv_hat_Sigma,nb_repli=500,parallel=FALSE)
head(Frequencies)
```

Dans la technique de "stability selection", nous allons sous-echantillonner plusieurs fois la moitié de nos observations. Ensuite, il calcule l'estimateur B à partir de ce sous-echantillon (en LASSO). À chaque itération on va regarder quelle variable explicative a été considérée comme informative. À la fin, on obtient le dataframe "Frequencies" avec la "fréquence de sélection" associé à chaque lipide et chaque génotype.

```{r}
Frequencies[Frequencies$frequency>=1,]
```


Le fait d'être un ppar va énormément jouer l'expression de l'acige gras C203n9 et C203n6.
```{r}
colnames(Frequencies)<-c('Names_of_Y','Names_of_X','frequency')
Frequencies$Names_of_X<-sub('X1','',Frequencies$Names_of_X)
Frequencies$Names_of_X
Frequencies$Names_of_Y
```

```{r}
p<-ggplot(data=Frequencies[Frequencies$frequency>=1,],aes(x=Names_of_Y,y=Names_of_X,color=frequency)) + geom_tile( aes(fill = frequency))+ylab('Levels of X')+xlab('m/z')+theme_bw() 
p

# +scale_color_gradient2(midpoint=0.99,mid ='orange')

```

```{r}
p<-ggplot(data=Frequencies[Frequencies$frequency>=1,],aes(x=Names_of_Y,y=Names_of_X,color=Names_of_X))+geom_point(size=1)+theme_bw()+ylab('Levels of X')+xlab('m/z') 
p
```

```{r}
#new_lipid <- data_lipid[, -c(1, 2)][,Frequencies[Frequencies$frequency==1,]$Names_of_Y]
#new_lipid <- nutrimouse$lipid[,Frequencies[Frequencies$frequency==1,]$Names_of_Y]
#new_X <- X[,Frequencies[Frequencies$frequency==1,]$Names_of_X]
#new_genotypes <- nutrimouse$genotype[,Frequencies[Frequencies$frequency==1,]$Names_of_Y]


```

```{r}
resultats <- lm(as.matrix(new_lipid) ~ X)
summary(resultats)
```


Nous pouvons refaire la meme chose avec les régimes cette fois ci, et non plus les génotypes. 

### Les acides gras selon le régime

Soit Y la matrice des acides gras
En ligne les individus en colonnes les differents acides gras.
Soit X2 la matrice de design des variables explicatives : les régimes. 

On suppose une indépendance des lignes sur les Ei.


```{r}
# NB : Y et X2 ont déjà été definis plus tot, pas besoin de les redefinir ne de rescale Y 
X <- model.matrix(lm(Y~X2+0)) # on crée la matrice de design à partir de X2 
p <- ncol(X) #2
n <- nrow(X) #40
q <- dim(Y)[2] #21

```


```{r}
residuals <- lm(Y~X-1)$residuals
pvalue = whitening_test(residuals)
pvalue
```
La pvalue est grande : il n'y a pas besoin de blanchir les données 


```{r}

# Donc pas besoin d'executer les lignes suivantes : 
#result=whitening_choice(residuals,c("AR1","nonparam","ARMA"),pAR=1,qMA=1)
#result

# Mais alors que vaut square_root_inv_hat_Sigma ??? --> matrice identité 
# On stocke la valeur de Sigma ^-1/2
#square_root_inv_hat_Sigma=whitening(residuals,"nonparam",pAR=1,qMA=0)
#square_root_inv_hat_Sigma

# Blanchiement / vectorisation / LASSO / Stability selection 
Frequencies_diet=variable_selection(Y,X, diag(21), nb_repli=500,parallel=FALSE)
head(Frequencies_diet)
```


```{r}
Frequencies_diet[Frequencies_diet$frequency>=1,]
```

```{r}
Frequencies_diet[Frequencies_diet$frequency>=0.998,]
```

Nous pouvons refaire la meme chose avec les gènes cette fois ci, et non plus les diet. 

### Les acides gras selon les gènes 

Soit Y la matrice des acides gras
En ligne les individus en colonnes les differents acides gras.
Soit X3 la matrice de design des variables explicatives : les gènes. 

On suppose une indépendance des lignes sur les Ei.


```{r}
 
X3 <- as.matrix(nutrimouse$gene)
dim(X3)

#X3 <- model.matrix(lm(Y~X3+0)) # on crée la matrice de design à partir de X2 
p <- ncol(X) #2
n <- nrow(X) #40
q <- dim(Y)[2] #21

```


```{r}
residuals <- lm(Y~X3)$residuals
residuals
pvalue = whitening_test(residuals)
pvalue
```