ModLinHD2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>\n",
    "<a href=\"http://www.insa-toulouse.fr/\" ><img src=\"http://www.math.univ-toulouse.fr/~besse/Wikistat/Images/logo-insa.jpg\" style=\"float:left; max-width: 120px; display: inline\" alt=\"INSA\"/></a> \n",
    "\n",
    "<a href=\"http://wikistat.fr/\" ><img src=\"http://www.math.univ-toulouse.fr/~besse/Wikistat/Images/wikistat.jpg\" style=\"max-width: 250px; display: inline\"  alt=\"Wikistat\"/></a>\n",
    "\n",
    "<a href=\"http://www.math.univ-toulouse.fr/\" ><img src=\"http://www.math.univ-toulouse.fr/~besse/Wikistat/Images/logo_imt.jpg\" style=\"float:right; max-width: 250px; display: inline\" alt=\"IMT\"/> </a>\n",
    "</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modèles linéaires et sélection de variables  en grande dimension avec <a href=\"https://cran.r-project.org/\"><img src=\"https://cran.r-project.org/Rlogo.svg\" style=\"max-width: 40px; display: inline\" alt=\"R\"/></a> : Etude de données d'expression de gènes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Objectifs **\n",
    "\n",
    "Sélection de variables en grande dimension dans le cadre du modèle linéaire. Le nombre de variables est ici très supérieur au nombre d'observations\n",
    "\n",
    "Comparaison sur le même jeu de données des méthodes suivantes :\n",
    "\n",
    "-  Algorithmes de sélection de modèles (critères AIC, BIC)\n",
    "- Régularisation ( Lasso, Elastic Net)\n",
    "- Régression sur composantes (PLS)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Les données\n",
    "\n",
    "Le jeu de  données `liver.toxicity`  contient, pour 64 individus :\n",
    "- 10 variables cliniques parmi lesquelles nous retiendrons une variable à expliquer $Y$ :  niveau de toxicité du foie chez le rat (il s'agit d'un niveau de cholestérol)\n",
    "-  3116 variables explicatives qui sont  des mesures d'expression de gènes (logarithme du rapport des niveaux d'expression entre deux conditions)\n",
    "\n",
    "\n",
    "** Source : Bushel, P., Wolfinger, R. D. and Gibson, G. (2007). Simultaneous clustering of gene expression\n",
    "data with clinical chemistry and pathological evaluations reveals phenotypic prototypes. BMC\n",
    "Systems Biology.**\n",
    "\n",
    "On est donc clairement ici dans un contexte de grande dimension : le nombre de variables est très supérieur au nombre d'individus. \n",
    "\n",
    "L'objectif est double : \n",
    "- Identifier les gènes qui ont un lien avec le niveau de cholestérol\n",
    "- Construire un modèle permettant de faire de la prédiction. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lecture des données"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(mixOmics)\n",
    "data(liver.toxicity)\n",
    "\n",
    "\n",
    "#load(\"liverToxicity.rdata\")\n",
    "\n",
    "# décrit les données\n",
    "gene=liver.toxicity$gene\n",
    "dim(gene)\n",
    "clinic=liver.toxicity$clinic\n",
    "dim(clinic)\n",
    "names(clinic)\n",
    "\n",
    "\n",
    "nomrats=dimnames(gene)[[1]]\n",
    "nomrats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# On retient les mesures d'expression des gènes et la variable clinique associée au taux de cholestérol : \n",
    "\n",
    "cholesterol=clinic[,10]\n",
    "liver=data.frame(cholesterol,gene)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Description des données"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "summary(liver[,1:5])\n",
    "\n",
    "hist(cholesterol)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hist(cor(gene,cholesterol))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Distribution de quelques covariables : \n",
    "\n",
    "par(mfrow=c(2,2))\n",
    "hist(liver[,2],main=\"Histogramme\",xlab=\"\")\n",
    "hist(liver[,3],main=\"Histogramme\",xlab=\"\")\n",
    "hist(liver[,4],main=\"Histogramme\",xlab=\"\")\n",
    "hist(liver[,5],main=\"Histogramme\",xlab=\"\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Matrice de correlations complete \n",
    "cormat<-cor(liver)\n",
    "dim(cormat)\n",
    "max(abs(cormat[1,-1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Variables les plus corrélées à la réponse\n",
    "which.max(abs(cormat[1,-1]))\n",
    "cormat[1,1158]\n",
    "\n",
    "colnames(liver)[1158]\n",
    "\n",
    "which.max(abs(cormat[1,c(-1,-1158)]))\n",
    "cormat[1,1837]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Graphes vs. la réponse correspondants\n",
    "par(mfrow=c(1,2))\n",
    "plot(liver[,1]~liver[,1158])\n",
    "plot(liver[,1]~liver[,1837])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sélection de  variables\n",
    "## Sélection par méthode  forward"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(leaps)\n",
    "#help(regsubsets)\n",
    "\n",
    "modselect=regsubsets(cholesterol~.,data=liver,method=\"forward\",nvmax=6,nbest=1,intercept=TRUE,really.big=TRUE)\n",
    "attributes(modselect)\n",
    "# afficher les variables selectionnees\n",
    "\n",
    "coef(modselect,1:6)\n",
    "\n",
    "### Selection forward: autre méthode pour afficher les variables selectionnees : pas très lisible ..\n",
    "varselect1=summary(modselect)$which\n",
    "varselect1\n",
    "\n",
    "mod1<-lm(cholesterol~.,data=liver[,varselect1[6,]])\n",
    "summary(mod1)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "selectforward=which(varselect1[6,]==TRUE)\n",
    "attributes(selectforward)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sélection par méthode \"sequential replacement\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(leaps)\n",
    "\n",
    "modselect=regsubsets(cholesterol~.,data=liver,method=\"seqrep\",nvmax=6,nbest=1,intercept=TRUE,really.big=TRUE)\n",
    "\n",
    "# nbest : nombre de modèle de chaque taille affiché\n",
    "\n",
    "# afficher les variables selectionnees\n",
    "\n",
    "coef(modselect,1:6)\n",
    "\n",
    "### Selection forward: autre méthode pour afficher les variables selectionnees\n",
    "varselect2=summary(modselect)$which\n",
    "\n",
    "mod2=lm(cholesterol~.,data=liver[,varselect2[6,]])\n",
    "summary(mod2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "selectforward=which(varselect1[6,]==TRUE)\n",
    "attributes(selectforward)\n",
    "\n",
    "selectseqrep=which(varselect2[6,]==TRUE)\n",
    "attributes(selectseqrep)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "library(gplots)\n",
    "\n",
    "data=list(X=selectforward,Y=selectseqrep)\n",
    "\n",
    "v=venn(data)\n",
    "attr(v, \"intersections\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Q ** Commenter ces résultats. Que concluez vous ? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a=selectforward[-1]\n",
    "a\n",
    "gene1=gene[,a-1]\n",
    "gene1\n",
    "\n",
    "b=selectseqrep[-1]\n",
    "b\n",
    "gene2=gene[,b-1]\n",
    "gene2\n",
    "\n",
    "cor(gene1,gene2)\n",
    "\n",
    "\n",
    "cor(cholesterol,gene1)\n",
    "cor(cholesterol,gene2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sélection de modèle par pénalisation Lasso avec la  Librairie glmnet\n",
    "\n",
    "L'utilisation de la librairie **`glmnet`** fournit des résultats plus rapides, ce qui peut s'avérer\n",
    "important pour des données de grande dimension. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Construction du modèle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(glmnet)\n",
    "\n",
    "x<-as.matrix(liver[,-1])\n",
    "y<-liver$cholesterol\n",
    "lasso<-glmnet(x,y)\n",
    "\n",
    "par(mfrow=c(1,2))\n",
    "plot(lasso,xvar=\"lambda\")\n",
    "plot(lasso,xvar=\"norm\")\n",
    "\n",
    "a=cv.glmnet(x,y)\n",
    "\n",
    "lambda.opt=a$lambda.min\n",
    "lambda.opt\n",
    "\n",
    "app=glmnet(x,y,lambda=lambda.opt)\n",
    "app\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Liste des variables sélectionnées : "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "### Lasso: variables selectionnees pour la valeur de lambda sélectionnée par validation croisée\n",
    "\n",
    "listevar<-predict(lasso,s=lambda.opt,type=\"nonzero\")# liste des variables avec un coefficient non-nul, pour un lambda donne\n",
    "colnames(liver)[listevar[1:(dim(listevar)[1]),]+1]# liste des noms de ces variables selectionnees\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modèle sélectionné"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geneselect=liver[,listevar[1:(dim(listevar)[1]),]+1]\n",
    "cholesterol=liver$cholesterol\n",
    "liverselect=data.frame(cholesterol,geneselect)\n",
    "modselect=lm(cholesterol~.,data= liverselect)\n",
    "summary(modselect)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colnames(geneselect)\n",
    "cor(cholesterol,geneselect)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sélection de modèle par pénalisation avec Elastic Net"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# on peut faire varier avec le paramètre alpha de glmnet\n",
    "out.elnet <- glmnet(x,y,alpha=0.5)\n",
    "\n",
    "a.elnet=cv.glmnet(x,y,alpha=0.5)\n",
    "lambda.opt=a.elnet$lambda.min\n",
    "app=glmnet(x,y,lambda=lambda.opt)\n",
    "\n",
    "\n",
    "app\n",
    "\n",
    "listevar<-predict(app,s=lambda.opt,type=\"nonzero\")# liste des variables avec un coefficient non-nul, pour un lambda donne\n",
    "colnames(liver)[listevar[1:(dim(listevar)[1]),]+1]# liste des noms de ces variables selectionnees\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a=colnames(geneselect)\n",
    "b=colnames(liver)[listevar[1:(dim(listevar)[1]),]+1]\n",
    "a\n",
    "b\n",
    "\n",
    "data=list(X=a,Y=b)\n",
    "\n",
    "v=venn(data)\n",
    "attr(v, \"intersections\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Q ** Commenter les résultats obtenus"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Régression PLS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "library(pls) \n",
    "#nombre optimal de composantes par validation croisée\n",
    "simpls= mvr(cholesterol~.,data=liver, ncomp=20, validation=\"CV\", method=\"simpls\")\n",
    "summary(simpls)\n",
    "\n",
    "#noter le nombre optimal de composantes \n",
    "plot(simpls,\"val\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Calcul des prévisions\n",
    "\n",
    "predapp.pls=predict(simpls,ncomp=10)\n",
    "resapp.pls=predapp.pls-cholesterol\n",
    "\n",
    "#Erreur d'apprentissage\n",
    "mean(resapp.pls**2)\n",
    "\n",
    "#graphique\n",
    "plot(predapp.pls,cholesterol)\n",
    "abline(0,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Q **  La méthode ne permet pas directement de selectionner les variables : on peut  étudier les corrélations des composantes PLS avec les variables initiales. Il est possible aussi d'utiliser la \"sparse PLS\". "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.3.3"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}