Ce TP est basé sur un document initialement rédigé par Boris Hejblum et repris par Robin Genuer.

Toxicité du paracétamol chez le rat

Charger mixOmics et les données correspondant à notre étude

library(mixOmics)
## Loading required package: MASS
## Loading required package: lattice
## Loading required package: ggplot2
## 
## Loaded mixOmics 6.3.2
## 
## Thank you for using mixOmics!
## 
## How to apply our methods: http://www.mixOmics.org for some examples.
## Questions or comments: email us at mixomics[at]math.univ-toulouse.fr  
## Any bugs? https://bitbucket.org/klecao/package-mixomics/issues
## Cite us:  citation('mixOmics')
data(liver.toxicity)

(Bushel, Wolfinger, and Gibson 2007)

Le sch?ma de l’étude est le suivant:

Travail demandé

On ne va considérer que le problème de régression. Il faut donc expliquer, pour les 64 rats, les 10 variables biochimiques en fonction des 3116 gènes, ou plutôt des gènes intéressants.

Pour celà on propose le protocole suivant :

  • Tuner les différents param?tres de façon adequat afin de considérer un modèle précis.

  • Commenter les performances de ce modèle et donner les représentations classiques de ce dernier :

    • Représentation des individus, discriminer visuellement les individus grâce aux conditions expérimentales (doses de paracétamol et dates de prélèvement).

    • Représentation des variables.

    • Utiliser une cim si le temps vous le permet.

  • A chaque fois commenter les résultats obtenus.

Suivre ce protocole!

Extraction des données

On standardise les données extraites

X <- scale(liver.toxicity$gene) # On standardise les donn?es
y <- scale(liver.toxicity$clinic) # On standardise les donn?es

On regarde un peu leurs dimensions

p <- ncol(X)
q <- ncol(y)
n <- nrow(y)
print(dim(X))
## [1]   64 3116
print(dim(y))
## [1] 64 10
print(getwd())
## [1] "C:/Users/hl1/Dropbox/cours_m2/spls_regression"

Estimation des paramètres

On décide de faire une sPLS, car

On décide de raliser la paramétrisation du modèle grâce à un leave-one-out, car

Pour l’instant on ne fait de la sélection que selon X.

Premiere composante

On cherche le \(keep_X\) qui minimise l’erreur de prédiction moyenne sur la première composante.

keepX_s <- c(1,5,10,15,20,30,50)
error <- matrix(NA, length(keepX_s),q)
for(id in 1:length(keepX_s)){
  model <- spls(X,y,ncomp = 1,keepX = keepX_s[id],mode = "regression")
  res_validation_croisee <- perf(model, validation="loo")
  error[id,] <- res_validation_croisee$MSEP[,1]
}

matplot(keepX_s,error,type="l",lty=1,col=1:10,lwd=2)
points(keepX_s,rowMeans(error),lty=2,lwd=3,col=1,type="l")
abline(v=keepX_s[which(rowMeans(error)==min(rowMeans(error)))],lwd=3,lty=3)
legend("topright",
       legend = c(colnames(y),"Mean of MSEP errors","Min of Mean of MSEP errors"),
       col = c(1:10,1,1),lty =c(rep(1,10),2,3),lwd=c(rep(2,10),3,3))

Dans ce cas, il vient

\[ keep_X^{(1)} = 30 \]

pour une erreur de \(0.6794462\).

Deuxieme composante

keepx_1 <- 30
keepX_s <- c(1,5,10,15,20,30,50)
error_2 <- matrix(NA, length(keepX_s),q)
for(id in 1:length(keepX_s)){
  model <- spls(X,y,ncomp = 2,keepX = c(keepx_1,keepX_s[id]),mode = "regression")
  res_validation_croisee <- perf(model, validation="loo")
  error_2[id,] <- res_validation_croisee$MSEP[,2]
}

matplot(keepX_s,error_2,type="l",lty=1,col=1:10,lwd=2)
points(keepX_s,rowMeans(error_2),lty=2,lwd=3,col=1,type="l")
abline(v=keepX_s[which(rowMeans(error_2)==min(rowMeans(error_2)))],lwd=3,lty=3)
legend("topright",
       legend = c(colnames(y),"Mean of MSEP errors","Min of Mean of MSEP errors"),
       col = c(1:10,1,1),lty =c(rep(1,10),2,3),lwd=c(rep(2,10),3,3))

Dans ce cas, il vient

\[ keep_X^{(2)} = 10 \] pour une erreur de \(0.5778407\).

On voit bien que la seconde composante permet de prédire mieux la variable TBA.umol.L.. Il faut alors regarder l’évolution de l’erreur de prédiction entre un modèle à 1 composante et un modèle à 2 composantes :

\[ error_{1\ comp} = 0.6794462, \]

\[ error_{2\ comp} = 0.5778407. \]

La seconde composante permet de prédire mieux les dix variables. On la conserve donc. Dans la suite il faudrait regarder les composantes suivantes et évaluer l’évolution de cette erreur de prédiction. Une erreur qui augmente prouve que le modèle se concentre sur des détails spécifiques à l’individu, soit une carctéristique personnelle soit même du bruit, on parle alors de surraprentissage ou overfitting en anglais.

Donc, deux cas peuvent se présenter à la construction d’une nouvelle composante, l’erreur peut:

  • augmenter : c’est qu’il ne faut pas considérer la composante courante: On fait du surraprentissage. En plus celà implique nécessairement qu’il faut arrêter là la construction du modèle PLS.

  • diminuer : alors la composante courante est à inclure dans le modèle. Il faut donc refaire le même travail pour la composante suivante.

On trouve les erreurs suivantes :

\[ error_{3\ comp} = 0.5484291, \ error_{4\ comp} = 0.5311632,\ error_{5\ comp} = 0.5489561. \]

On voit bien que l’erreur réaugment à partir de la cinquièmee composante. Un choix raisonnable serait donc de ne conserver que \(4\) composantes dans le modèle.

Création du modèle à 2 composantes

On crée le modèle et on observe les gènes impliqués à droite et les vraiables prédites à droite.

model <- spls(X,y,ncomp=2,keepX = c(30,10))

plotLoadings(model, comp = 1, method = 'mean', contrib = 'max',legend=F)

plotLoadings(model, comp = 2, method = 'mean', contrib = 'max',legend=F)

On voit que les Y les mieux prédits par la première composante, AST.IU.L. et ALT.IU.L., sont les deux variables les plus corrélées avec la première composante. Ce sont aussi les variables les moins corrélées avec la seconde composante.

Interprétation des axes

Par rapport à la dose administrée

Y_dose <- as.factor(liver.toxicity$treatment$Dose.Group)
plot(model$variates$X,col=Y_dose,pch=16,
     main="Two first components according to dose")
legend("topleft",legend=levels(Y_dose),fill = 1:4)

On voit très nettement ques les deux premières composantes sont très sensibles au niveau de dose.

Par rapport au temps de prélèvement

Y_time <- as.factor(liver.toxicity$treatment$Time.Group)
plot(model$variates$X,col=as.numeric(Y_time)+4,pch=16,cex=1.5,
     main="Two first components according to time of necropsies")
legend("topleft",legend=levels(Y_time),fill = 1:4+4)

La première composante est sensible aux individus prélevés aux temps 18 et 24 alors que la seconde composante aux individus prélevés au temps 48. On peut dès lors supposer que l’état des souris aux temps 18 et 24 est partiellement le même et que celui-ci est bien différent des états aux temps 3 et 48.

References

Bushel, Pierre R, Russell D Wolfinger, and Greg Gibson. 2007. “Simultaneous Clustering of Gene Expression Data with Clinical Chemistry and Pathological Evaluations Reveals Phenotypic Prototypes.” BMC Systems Biology 1 (1). BioMed Central: 15.