ModLinHD.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 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>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Objectifs **\n",
    "Comparaison sur le même jeu de données des qualités de prévision de plusieurs modèles obtenus par:\n",
    "- Modèle linéaire\n",
    "-  Algorithmes de sélection de modèles (critères AIC, BIC)\n",
    "- Régularisation (ridge, Lasso, Elastic Net)\n",
    "- Régression sur composantes (PCR, PLS)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Les données\n",
    "\n",
    "Nous allons utiliser les données `Prostate` du package **`lasso 2 `** de R.\n",
    "\n",
    " These data come from a study that examined the correlation between the\n",
    "  level of prostate specific antigen and a number of clinical measures\n",
    "  in men who were about to receive a radical prostatectomy.  It is data\n",
    "  frame with 97 rows and 9 columns.\n",
    "\n",
    "** Source : Stamey, T.A., Kabalin, J.N., McNeal, J.E., Johnstone, I.M., Freiha,\n",
    "  F., Redwine, E.A. and Yang, N. (1989). *Prostate specific antigen in the diagnosis and treatment of\n",
    "  adenocarcinoma of the prostate: II. radical prostatectomy treated\n",
    "  patients* ,  Journal of Urology, {141}(5), 1076--1083. **\n",
    "\n",
    "Liste des 9 variables la dernière est à expliquer : \n",
    "- **`lcavol`** :  log(cancer volum)\n",
    "- **`lweight`** :  log(prostate weight)\n",
    "- **`age`** :   age\n",
    "- **`lbph `** :  log(benign prostatic hyperplasia amount)\n",
    "- **`svi`** :  seminal vesicle invasion\n",
    "- **`lcp`** :  log(capsular penetration)\n",
    "- **`lweight`** :  log(prostate weight)\n",
    "- **`gleason`** : Gleason score\n",
    "- **`pgg45 `** :   percentage Gleason scores 4 or 5\n",
    "\n",
    "- **`lpsa `** :  log(prostate specific antigen)\n",
    "\n",
    "\n",
    "Les variables **`svi`**  et **`gleason`** sont qualitatives.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Lecture des données"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(lasso2)\n",
    "data(Prostate)\n",
    "# Une seule valeur 8 du score de Gleason ; on regroupe  les modalités 8 et 9 ensembles\n",
    "Prostate$gleason[Prostate$gleason==8]=9\n",
    "# On transforme svi et gleason en variables qualitatives\n",
    "Prostate$svi=as.factor(Prostate$svi)\n",
    "Prostate$gleason=as.factor(Prostate$gleason)\n",
    "summary(Prostate)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Description des données"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dim(Prostate)\n",
    "names(Prostate)\n",
    "summary(Prostate)\n",
    "cor(Prostate[,-c(5,7)])\n",
    "hist(Prostate$lpsa,10)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Création d'un échantillon d'apprentissage et d'un échantillon test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Les valeurs de lpsa sont rangées par ordre croissant\n",
    "#On conserve 1/4 des données pour l'échantillon test \n",
    "ind.test=4*c(1:22)\n",
    "Prostate.app=Prostate[-ind.test,]\n",
    "Prostate.test=Prostate[c(ind.test),]\n",
    "dim(Prostate.test)\n",
    "dim(Prostate.app)\n",
    "ntest=length(Prostate.test$lpsa)\n",
    "napp=length(Prostate.app$lpsa)\n",
    "summary(Prostate.app)\n",
    "summary(Prostate.test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modèle linéaire complet"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Estimation du modèle et graphes des résidus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Une fonction utile de graphe des résidus.\n",
    "\n",
    "plot.res=function(x,y,titre=\"\")\n",
    "{\n",
    "plot(x,y,col=\"blue\",ylab=\"Résidus\",\n",
    "   xlab=\"Valeurs predites\",main=titre)\n",
    "abline(h=0,col=\"green\")\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "modlin=lm(lpsa~., data=Prostate.app)\n",
    "summary(modlin) # noter les p-valeurs\n",
    "\n",
    "#Residus \n",
    "res=residuals(modlin)\n",
    "\n",
    "#Regroupement des graphiques sur la meme page\n",
    "par(mfrow=c(1,2))\n",
    "hist(residuals(modlin))\n",
    "qqnorm(res)\n",
    "qqline(res, col = 2)\n",
    " \n",
    "par(mfrow=c(1,1)) # retour au graphique standard\n",
    "plot.res(predict(modlin),res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Erreur d'apprentissage et erreur de généralisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Erreur d'apprentissage\n",
    "\n",
    "mean(res**2)\n",
    "\n",
    "#Erreur sur l'échantillon test\n",
    "\n",
    "pred.test=predict(modlin, newdata=Prostate.test)\n",
    "res.test=pred.test-Prostate.test$lpsa\n",
    "mean(res.test**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Nouvelle paramétrisation \n",
    "\n",
    "Afin de faciliter l'interprétation des résultats concernant les variables qualitatives, on introduit une nouvelle paramétrisation à l'aide de contrastes. \n",
    "\n",
    "Par défaut, la référence est prise pour la valeur 0 de  **`svi`** et 6 de  **`gleason`**, qui sont les plus petites valeurs. Les paramètres indiqués pour les variables svi1, gleason 7, 8 et 9 indiquent l'écart estimé par rapport à cette référence. Il est plus intéressant en pratique de se référer à la moyenne des observations sur toutes les modalités des variables qualitatives, et d'interpréter les coefficients comme des écarts à cette moyenne. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "contrasts(Prostate.app$svi)=\n",
    "contr.sum(levels(Prostate.app$svi)) \n",
    " \n",
    "contrasts(Prostate.app$gleason)=\n",
    "contr.sum(levels(Prostate.app$gleason)) \n",
    "\n",
    "modlin2=lm(lpsa~., Prostate.app)\n",
    "summary(modlin2)\n",
    " \n",
    "# Attention au nom des variables : gleason1 =6,  gleason2 =7, gleason3 =8-9 (pas affiché),\n",
    "#la somme des coefficients associés à ces variables est nulle\n",
    "#svi1=0, svi2=1 (pas affiché), la somme des deux coefficients est nulle. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sélection de modèle par sélection de variables\n",
    "### Sélection par AIC et backward"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(MASS)\n",
    "\n",
    "modselect_b=stepAIC(modlin2,~.,trace=TRUE,\n",
    "direction=c(\"backward\"))\n",
    "\n",
    "summary(modselect_b)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sélection par AIC et forward"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod0=lm(lpsa~1,data=Prostate.app)\n",
    "\n",
    "modselect_f=stepAIC(mod0,lpsa~lcavol+lweight\n",
    "+age+lbph+svi+lcp+gleason+pgg45,\n",
    "data=Prostate.app,trace=TRUE,direction=c(\"forward\"))\n",
    "\n",
    "summary(modselect_f)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sélection par AIC et stepwise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "modselect=stepAIC(modlin2,~.,trace=TRUE,\n",
    "direction=c(\"both\"))\n",
    "#both est l'option par défaut\n",
    "\n",
    "summary(modselect)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sélection par BIC et stepwise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#  k=log(napp) pour BIC au lieu de AIC.\n",
    "\n",
    "modselect_BIC=stepAIC(modlin2,~.,trace=TRUE,\n",
    "direction=c(\"both\"),k=log(napp))\n",
    "\n",
    "summary(modselect_BIC)\n",
    "#Le modèle sélectionné est plus parcimonieux\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calcul de l'erreur d'apprentissage "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Modèle stepwise AIC\n",
    "\n",
    "mean((predict(modselect)-Prostate.app[,\"lpsa\"])**2)\n",
    "\n",
    "#Modèle stepwise BIC\n",
    "mean((predict(modselect_BIC)-Prostate.app[,\"lpsa\"])**2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Q.** Commenter ces résultats : laquelle des 2 erreurs est la plus grande.  Pouvait on le prévoir ?  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calcul de l'erreur de test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Modèle stepwise AIC\n",
    "\n",
    "mean((predict(modselect,newdata=Prostate.test)-Prostate.test[,\"lpsa\"])**2)\n",
    "\n",
    "#Modèle stepwise BIC\n",
    "\n",
    "mean((predict(modselect_BIC,newdata=Prostate.test)-Prostate.test[,\"lpsa\"])**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Q.** Commenter ces résultats : laquelle des 2 erreurs est la plus grande.  Pouvait on le prévoir ? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sélection de modèle par pénalisation Ridge"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comportement des coefficients\n",
    "\n",
    "Calcul des coefficients pour différentes valeurs du paramètre lambda"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(MASS)\n",
    "mod.ridge=lm.ridge(lpsa~.,data=Prostate.app,\n",
    "lambda=seq(0,20,0.1))\n",
    "par(mfrow=c(1,1))\n",
    "plot(mod.ridge)\n",
    "\n",
    "matplot(t(mod.ridge$coef),lty=1:3,type='l',col=1:10)\n",
    "legend(\"top\",legend=rownames(mod.ridge$coef),\n",
    "col=1:10,lty=1:3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pénalisation optimale par validation croisée"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "select(mod.ridge) # noter la valeur puis estimer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mod.ridgeopt=lm.ridge(lpsa ~  .,data=Prostate.app,lambda=9.8)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Q. ** Automatiser le calcul du modèle avec la valeur optimale de lambda  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prévision  sur l'échantillon d'apprentissage\n",
    "\n",
    "La  fonction `predict.ridgelm` n'existe pas, il faut calculer les valeurs prédites à partir des coefficients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Coefficients du modèle sélectionné :\n",
    "\n",
    "coeff=coef(mod.ridgeopt)\n",
    "\n",
    "#On crée des vecteurs pour les variables  qualitatives\n",
    "\n",
    "svi0.app=1*c(Prostate.app$svi==0)\n",
    "svi1.app=1-svi0.app\n",
    "gl6.app=1*c(Prostate.app$gleason==6)\n",
    "gl7.app=1*c(Prostate.app$gleason==7)\n",
    "gl9.app=1*c(Prostate.app$gleason==9)\n",
    "\n",
    "#variables quantitatives\n",
    "\n",
    "lcavol.app=Prostate.app$lcavol\n",
    "lweight.app=Prostate.app$lweight\n",
    "age.app=Prostate.app$age\n",
    "lbph.app=Prostate.app$lbph\n",
    "lcp.app=Prostate.app$lcp\n",
    "pgg45.app=Prostate.app$pgg45\n",
    "\n",
    "\n",
    "#Calcul des valeurs prédites\n",
    "fit.rid=rep(coeff[1],napp)+coeff[2]*lcavol.app+coeff[3]*lweight.app+coeff[4]*age.app+coeff[5]*lbph.app+coeff[6]*svi0.app-coeff[6]*svi1.app+coeff[7]*lcp.app+coeff[8]*gl6.app+coeff[9]*gl7.app-(coeff[8]+coeff[9])*gl9.app+coeff[10]*pgg45.app\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Tracé des valeurs prédites en fonctions  des valeurs observées : **"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(Prostate.app$lpsa,fit.rid)\n",
    "abline(0,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calcul et tracé des résidus\n",
    "res.rid=fit.rid-Prostate.app[,\"lpsa\"]\n",
    "plot.res(fit.rid,res.rid,titre=\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Erreur d'apprentissage\n",
    "mean(res.rid**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prévision sur l'échantillon test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "svi0.test=1*c(Prostate.test$svi==0)\n",
    "svi1.test=1-svi0.test\n",
    "gl6.test=1*c(Prostate.test$gleason==6)\n",
    "gl7.test=1*c(Prostate.test$gleason==7)\n",
    "gl9.test=1*c(Prostate.test$gleason==9)\n",
    "\n",
    "ntest=length(Prostate.test$lpsa)\n",
    "prediction=rep(coeff[1],ntest)+coeff[2]*Prostate.test$lcavol+coeff[3]*Prostate.test$lweight+coeff[4]*Prostate.test$age+coeff[5]*Prostate.test$lbph+coeff[6]*svi0.test-coeff[6]*svi1.test+coeff[7]*Prostate.test$lcp+coeff[8]*gl6.test+coeff[9]*gl7.test-(coeff[8]+coeff[9])*gl9.test+coeff[10]*Prostate.test$pgg45\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Erreur sur l'échantillon test **"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean((Prostate.test[,\"lpsa\"]-prediction)^2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "par(mfrow=c(1,1))\n",
    "plot(mod.ridge)\n",
    "abline(v=9.8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Q. ** Commenter les résultats et comparer avec les modèles précédents"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sélection de modèle par pénalisation Lasso avec la  Librairie Lasso2 \n",
    "\n",
    "### Construction du modèle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "library(lasso2)\n",
    "l1c.P <- l1ce(lpsa ~ ., Prostate.app, bound=(1:100)/100, absolute.t=FALSE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La borne est ici relative, elle correspond à une certaine proportion de la norme $\\mathbb{L}_1$ du vecteur des \n",
    " coefficients des moindres carrés. Une borne égale à $1$ correspond donc  à l'absence de pénalité, on retrouve \n",
    "l'estimateur des  moindres carrés. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualisation des coefficients : chemins de régularisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coefficients=coef(l1c.P)\n",
    "plot(l1c.P,col=1:11,lty=1:3,type=\"l\")\n",
    "legend(\"topleft\",legend=colnames(coefficients),col=1:11,lty=1:3)\n",
    "\n",
    "#On supprime le terme constant\n",
    "penalite_relative=c(1:100)/100\n",
    "matplot(penalite_relative,coefficients[,-1],lty=1:3,type='l',col=1:10)\n",
    "legend(\"topleft\",legend=colnames(coefficients[,-1]),col=1:10,lty=1:3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sélection de la pénalité par validation croisée"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "vc=gcv(l1c.P)\n",
    "crit.vc=vc[,\"gcv\"]\n",
    "bound_opt=vc[which.min(crit.vc),\"rel.bound\"]\n",
    "\n",
    "l1c.opt <- l1ce(lpsa ~ ., Prostate.app, \n",
    "bound=bound_opt, absolute.t=FALSE)\n",
    "coef=coef(l1c.opt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "### Erreur d'apprentissage et erreur de généralisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#erreur apprentissage\n",
    "fit=fitted(l1c.opt)\n",
    "mean((fit-Prostate.app[,\"lpsa\"])^2)\n",
    "\n",
    "\n",
    "#Erreur sur l'échantillon test}\n",
    "prediction=predict(l1c.opt,newdata=Prostate.test)\n",
    "mean((prediction-Prostate.test[,\"lpsa\"])^2)"
   ]
  },
  {
   "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. Par contre, on ne peut pas traiter à priori des variables qualitatives.\n",
    " Nous allons donc devoir créer des vecteurs avec des variables indicatrices des diverses modalités \n",
    "pour les variables qualitatives. Nous ne prendrons pas en compte les contrastes. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mise en forme des variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#on construit une matrice xx.app d'apprentissage et xx.test de test\n",
    "# On recharge les données qui sont quantitatives\n",
    "data(Prostate)\n",
    "Prostate$gleason[Prostate$gleason==8]=9\n",
    "Prostate.app=Prostate[-ind.test,]\n",
    "Prostate.test=Prostate[c(ind.test),]\n",
    "\n",
    "x.app=Prostate.app[,-9]\n",
    "y.app=Prostate.app[,9]\n",
    "x.app=as.matrix(x.app)\n",
    "\n",
    "#on construit une matrice avec les vecteurs indicatrices \n",
    "gl7.app=1*c(Prostate.app$gleason==7)\n",
    "gl9.app=1*c(Prostate.app$gleason==9)\n",
    "\n",
    "xx.app=matrix(0,ncol=9,nrow=nrow(x.app))\n",
    "xx.app[,1:6]=as.matrix(x.app[,1:6])\n",
    "xx.app[,7:8]=cbind(gl7.app,gl9.app)\n",
    "xx.app[,9]=cbind(x.app[,8])\n",
    "\n",
    "#on nomme les colonnes avec le noms des variables\n",
    "colnames(xx.app)=c(\"lcavol\",\"lweight\", \"age\",\"lbph\",\"svi1\",\"lcp\",\"gl7\",\"gl9\",\"pgg45\")\n",
    "\n",
    "#on fait de meme pour l'echantillon test\n",
    "x.test=Prostate.test[,-9]\n",
    "y.test=Prostate.test[,9]\n",
    "x.test=as.matrix(x.test)\n",
    "gl7.test=1*c(Prostate.test$gleason==7)\n",
    "gl9.test=1*c(Prostate.test$gleason==9)\n",
    "\n",
    "#on construit une matrice avec les vecteurs indicatrices \n",
    "xx.test=matrix(0,ncol=9,nrow=nrow(x.test))\n",
    "xx.test[,1:6]=x.test[,1:6]\n",
    "xx.test[,7:8]=cbind(gl7.test,gl9.test)\n",
    "xx.test[,9]=x.test[,8]\n",
    "colnames(xx.test)=colnames(xx.app)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Construction du modèle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(glmnet)\n",
    "out.lasso <- glmnet(xx.app,y.app)\n",
    "l=length(out.lasso$lambda)\n",
    "b=coef(out.lasso)[-1,1:l]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualisation des coefficients  : chemins de régularisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "matplot(t(as.matrix(out.lasso$beta)),type='l',\n",
    "col=1:10,lty=1:3)\n",
    "legend(\"topleft\",legend=colnames(xx.app),\n",
    "col=1:10,lty=1:3)\n",
    "title(\"Lasso\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sélection de la pénalité par validation croisée"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "a=cv.glmnet(xx.app,y.app)\n",
    "\n",
    "lambda.opt=a$lambda.min\n",
    "app=glmnet(xx.app,y.app,lambda=lambda.opt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "### Erreur d'apprentissage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "appr=predict(app,newx=xx.app)\n",
    "mean((appr-Prostate.app[,9])^2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "### Erreur sur l'échantillon test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred=predict(app,newx=xx.test)\n",
    "mean((pred-Prostate.test[,9])^2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Q. ** Comparer les résultats  et le modèle obtenu avec ce qui a été obtenu par la librairie Lasso 2. "
   ]
  },
  {
   "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(xx.app,y.app,alpha=0.5)\n",
    "\n",
    "a.elnet=cv.glmnet(xx.app,y.app,alpha=0.5)\n",
    "lambda.opt=a.elnet$lambda.min\n",
    "app=glmnet(xx.app,y.app,lambda=lambda.opt)\n",
    "\n",
    "#erreur apprentissage\n",
    "app.elnet=predict(a.elnet,newx=xx.app)\n",
    "mean((app.elnet-Prostate.app[,9])^2)\n",
    "\n",
    "#erreur de prédiction\n",
    "predi.elnet=predict(a.elnet,newx=xx.test)\n",
    "mean((predi.elnet-Prostate.test[,9])^2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sélection de modèle projection sur composantes orthogonales\n",
    "\n",
    "### Régression PLS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data(Prostate)\n",
    "Prostate$gleason[Prostate$gleason==8]=9\n",
    "Prostate$svi=as.factor(Prostate$svi)\n",
    "Prostate$gleason=as.factor(Prostate$gleason)\n",
    "Prostate.app=Prostate[-ind.test,]\n",
    "Prostate.test=Prostate[c(ind.test),]\n",
    "\n",
    "library(pls) \n",
    "#nombre optimal de composantes par validation croisée\n",
    "simpls= mvr(lpsa~.,data=Prostate.app, ncomp=9, validation=\"CV\", method=\"simpls\")\n",
    "summary(simpls)\n",
    "#graphique\n",
    "plot(simpls)\n",
    "abline(0,1)\n",
    "\n",
    "#noter le nombre optimal de composantes \n",
    "plot(simpls,\"val\")\n",
    "\n",
    "#Avec leave-one -out\n",
    "simplsloo= mvr(lpsa~.,data=Prostate.app, ncomp=9, \n",
    "   validation=\"LOO\", method=\"simpls\")\n",
    "summary(simplsloo)\n",
    "\n",
    "#On sélectionne 6 composantes\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Calcul des prévisions\n",
    "predapp.pcr=predict(simpls,ncomp=5)\n",
    "resapp.pcr=predapp.pcr-Prostate.app$lpsa\n",
    "\n",
    "#Erreur d'apprentissage\n",
    "mean(resapp.pcr**2)\n",
    "\n",
    "#Valeurs prédites sur l'échantillon test\n",
    "pred.pls=predict(simpls,newdata=Prostate.test,ncomp=5)\n",
    "#Erreur de test\n",
    "mean((pred.pls-Prostate.test[,\"lpsa\"])**2) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Régression sur composantes principales"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "mod.pcr = pcr(lpsa~.,data=Prostate.app, ncomp=9, \n",
    "validation=\"CV\")\n",
    "summary(mod.pcr) \n",
    "# noter le nombre optimal\n",
    "plot(mod.pcr,\"val\")\n",
    "#On sélectionne 8 composantes\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Calcul des prévisions\n",
    "predapp.pls=predict(mod.pcr,ncomp=8)\n",
    "resapp.pls=predapp.pls-Prostate.app$lpsa\n",
    "#Erreur d'apprentissage\n",
    "mean(resapp.pls**2)\n",
    "\n",
    "\n",
    "#Valeurs prédites sur l'échantillon test\n",
    "pred.pcr=predict(mod.pcr,newdata=Prostate.test,ncomp=8)\n",
    "#Erreur de test\n",
    "mean((pred.pcr-Prostate.test[,\"lpsa\"])**2) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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
}