NIR/Apprent-R-cookies.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=\"float:right; max-width: 250px; display: inline\"  alt=\"Wikistat\"/></a>\n",
    "\n",
    "</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [Scénarios d'Apprentissage Statistique](https://github.com/wikistat/Apprentissage)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Discrétisation de spectres en proche infra-rouge (NIR) 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> \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Résumé\n",
    "Modélisation (ou calibration) en grande dimension $p>n$ où les variables explicatives sont les discrétisations de spectres en proche infra rouge (NIR) et la variable à expliquer la part de sucre dans une pâte à gâteau. Différentes méthodes de régression ([ridge, lasso, lars, elasticnet](http://wikistat.fr/pdf/st-m-app-select.pdf), [svm](http://wikistat.fr/pdf/st-m-app-svm.pdf), [pls](http://wikistat.fr/pdf/st-m-app-sparse-pls.pdf), sur composantes principales}...) adaptées à cette situation sont utilisés et leur qualités prédictives comparées. La librairie `caret` est ensuite utilisée pour industrialiser la stratégie de choix de modèle et de méthode. \n",
    "\n",
    "Un deuxième scénario sur des données NIR ([Tecator](http://wikistat.fr/pdf/st-scenar-app-tecator.pdf)) aborde des méthodes non linéaires et l'usage des splines pour lisser et calculer les dérivées."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1 Introduction\n",
    "### 1.1 Contexte\n",
    "Le travail présenté s'intéresse à un problème de contrôle de qualité sur une chaîne de fabrication de biscuits (cookies). Il est nécessaire de contrôler le mélange des ingrédients avant cuisson afin de s'assurer que les proportions en lipides, sucre, farine, eau, sont bien respectées c'est-à-dire proches des valeurs nominales de la recette qui a fait la réputation de l'entreprise avant l'industrialisation de la production. Il s'agit de savoir s'il est possible de dépister au plus tôt une dérive afin d'intervenir sur les équipements concernés. Les mesures et analyses, faites dans un laboratoire classique de chimie, sont relativement longues et coûteuses ; elles ne peuvent être entreprises pour un suivi régulier ou même en continue de la production. C'est pour cela qu'il à été décidé l'utilisation d'un spectromètre en proche infrarouge (NIR). \n",
    "\n",
    "Ce type de problème est très classique en agro-alimentaire et plus généralement dans l'industrie et correspond à une étude dite de calibration en chimiométrie. Une très abondante littérature lui est consacrée. \n",
    "### 1.2 Données\n",
    "Les données originales sont dues à Osbone et al. (1984) et ont été souvent utilisées pour la comparaison de méthodes (Stone et al. 1990, Brown et al. 2001, Krämer et al. 2008). Elles sont  accessibles dans R au sein du package `ppls`. \n",
    "Les mesures ont été faites sur deux échantillons, l'un de taille 40 prévu pour l'apprentissage, l'autre de taille 32 pour les tests. Pour chacun de ces 72 biscuits, les compositions en lipides, sucre, farine, eau, sont mesurées par une approche classique tandis que le spectre est observé  sur toutes les longueurs d'ondes entre 1100 et 2498 nanomètres, régulièrement espacés de 2 nanomètres. Nous avons donc 700 valeurs observées, ou variables potentiellement explicatives, par échantillon de pâte à biscuit. \n",
    "\n",
    "Typiquement, cette étude se déroule dans un contexte de grande dimension avec $p>n$.\n",
    "\n",
    "### 1.3 Objectif\n",
    "L'objectif principal est de répondre à la question suivante : est-il possible de prévoir ces quantités à partir des spectres ? En cas de réponse positive, le gain de temps, donc d'argent, serait immédiat. L'étude est restreinte à  la seule modélisation du taux de sucre pour la recherche d'un meilleur modèle de prévision. Il s'agit donc d'évaluer, comparer différentes stratégies et méthodes pour aboutir à celle la plus efficace. Le travail est découpé en deux parties, dans la première sont testées différentes méthodes en construisant les modèles à la main ou pas à pas tandis que la 2ème partie fait appel à une librairie `caret` qui propose une comparaison automatique d'optimisation des paramètres pour un vaste ensemble de méthodes.\n",
    "\n",
    "**Répondre** aux questions **Q**.\n",
    "\n",
    "#### Références\n",
    "\n",
    "Nicole Krämer, Anne Laure Boulesteix et Gerhard Tutz (2008). Penalized Partial Least Squares with applications to B-spline transformations and functional  data, *Chemometrics  and  Intelligent  Laboratory  Systems* 94, 60–69.\n",
    "\n",
    "Max Kuhn (2008). Building Predictive Models in R Using the caret Package, *Journal of Statistical Software*, 28, 5.\n",
    "\n",
    "B. G. Osborne, T. Fearn, A. R. Miller et S. Douglas (1984). Application of Near Infrared Reflectance spectroscopy to the compositional analysis of biscuits and biscuit doughs, *J. Sci. Food Agric.* 35, 99–105.\n",
    "\n",
    "## 2 Approche exploratoire\n",
    "### 2.1 Prise en charge des données"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(ppls)\n",
    "data(cookie)\n",
    "# extraire le taux de sucre et les spectres\n",
    "cook = data.frame(cookie[,702],cookie[,1:700]) \n",
    "names(cook)= c(\"sucre\",paste(\"X\",1:700,sep=\"\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Statistiques élémentaires"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# la variable à expliquer\n",
    "options(repr.plot.width=4, repr.plot.height=4)\n",
    "hist(cook[,\"sucre\"], main=\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** La figure ci-dessus représente la distribution des valeurs de concentration en sucre. Commentez en quelques mots.\n",
    "\n",
    "**Q** Les figures ci-dessous représentent les différents spectres (l'absorbance) en fonction de la fréquence (pas du temps) pour chaque biscuit, puis leur distribution. La couleur est fonction de la concentration de sucre. Que dire de ces courbes, de leur apparence commune, de leur régularité et de *l'influence sur les corrélations* des variables après discrétisation ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# les spectres colorés par l'intensité en sucre\n",
    "coul=rainbow(20)[as.integer(as.factor(as.integer(cook[,1]-10)))]\n",
    "ts.plot(t(cook[,-1]),col=coul)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Les diagrammes boîtes puis variances des variables explicatives\n",
    "options(repr.plot.width=4, repr.plot.height=4)\n",
    "boxplot(cook[,-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hist(var(cook[,-1]), main=\"\") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 Analyse en Composantes Principales\n",
    "\n",
    "**Q** Les graphes ci-dessous (éboulis, biplot, vecteurs propres) sont produits par une analyse en composantes principales des spectres discrétisés. Combien d'axes retenir ? Interpréter brièvement le premier axe, le 2ème."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "acp=prcomp(cook[,-1])\n",
    "options(repr.plot.width=4, repr.plot.height=3)\n",
    "plot(acp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=5, repr.plot.height=5)\n",
    "biplot(acp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=4, repr.plot.height=4)\n",
    "plot(acp$x,col=coul,pch=19)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot.ts(acp$rotation[,1:10], main=\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Quel traitement préalable sous contrôle de l'acp serait envisageable sur ces courbes pour une étude exploratoire? Néanmoins, serait-il pertinent au vu de l'objectif de modélisation de ces données?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Une fonction sera utilisée pour une représentation homogène des résidus des modèles :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plot.res=function(x,y,titre=\"titre\")\n",
    "{\n",
    "plot(x,y,col=\"blue\",ylab=\"Résidus\",xlab=\"Valeurs predites\",main=titre,pch=19)\n",
    "abline(h=0,col=\"green\")\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3 Modélisation artisanale\n",
    "\n",
    "La modélisation est construite méthode par méthode pour mieux en assimiler le déroulement. Chaque méthode de modélisation possède des spécificités, notamment dans la manière d'optimiser plus ou moins facilement les valeurs des paramètres. Il est indispensable de bien se familiariser avec ces différents éléments.\n",
    "### 3.1 Les échantillons\n",
    "Les 40 premières lignes sont usuellement considérées pour l'apprentissage tandis que les 32 dernières sont l'échantillon test. C'est fait pour éventuellement comparer avec des résultats de la littérature. \n",
    "\n",
    "**Q** Compte tenu de l'objectif de modélisation et de la structure des données, peut-on envisager l’estimation d’un modèle linéaire gaussien? Pourquoi?\n",
    "\n",
    "**Q** Une approche classique de type \"sélection de variables\" en régression permettrait-elle de résoudre ce problème? Pourquoi?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# extraire apprentissage et test\n",
    "cook.app=cook[1:40,]\n",
    "cook.test=cook[41:72,]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Régression ridge\n",
    "**Q** Quel est l'objectif de cette régression?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(MASS)\n",
    "# paramètres en fonction de la pénalisation\n",
    "plot(lm.ridge(sucre ~ ., data=cook.app,lambda = seq(0,1,0.01)))\n",
    "cook.ridge=lm.ridge(sucre ~ ., data=cook.app,lambda = seq(0,0.2,0.001))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Comment optimiser la valeur de `lambda`?\n",
    "\n",
    "**Q** Quelle stratégie est mise en place ci-dessous? Que fait la fonction `choix.kappa`?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sélection des indices de validation croisée\n",
    "library(pls)\n",
    "set.seed(87)\n",
    "cvseg=cvsegments(nrow(cook.app),k=4,type=\"random\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "choix.kappa=function(kappamax,cvseg,nbe=100){\n",
    "press=rep(0,nbe)\n",
    "for (i in 1:length(cvseg)){\n",
    "valid=cook.app[unlist(cvseg[i]),]\n",
    "modele=lm.ridge(sucre~.,data=cook.app[unlist(cvseg[-i]),],lambda=seq(0,kappamax,length=nbe))\n",
    "coeff=coef(modele)\n",
    "prediction=matrix(coeff[,1],nrow(coeff),nrow(valid))+coeff[,-1]%*%t(data.matrix(valid[,-1])) \n",
    "press=press+rowSums((matrix(valid[,1],nrow(coeff),nrow(valid),byrow=T)-prediction)^2)\n",
    "      }\n",
    "kappaet=seq(0,kappamax,length=nbe)[which.min(press)]\n",
    "return(list(kappaet=kappaet,press=press))\n",
    "   }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# exécution\n",
    "res=choix.kappa(0.5,cvseg,nbe=1000)\n",
    "options(repr.plot.width=3, repr.plot.height=3)\n",
    "plot(seq(0,0.5,length=1000),res$press,pch=19)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Comment interpréter le résulat?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# valeur optimale\n",
    "kappaet=res$kappaet\n",
    "cook.ridgeo=lm.ridge(sucre~.,data=cook.app,lambda=kappaet)\n",
    "coeff=coef(cook.ridgeo)\n",
    "# Calcul des valeurs ajustées et des résidus\n",
    "fit.rid=rep(coeff[1],nrow(cook.app))+as.vector(coeff[-1]%*%t(data.matrix(cook.app[,-1])))\n",
    "plot(fit.rid,cook.app[,\"sucre\"],pch=19)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** La régression mettant en oeuvre le paramètre optimal fournit le résultat ci-dessus. Commenter la qualité de l'ajustement. Que vaut approximativement le $R^2$?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res.rid=fit.rid-cook.app[,\"sucre\"]\n",
    "options(repr.plot.width=4, repr.plot.height=4)\n",
    "plot.res(fit.rid,res.rid,titre=\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Commenter le graphique précédent.\n",
    "\n",
    "**Q** Quels sont les calculs ci-dessous? Noter le résultat pour comparer avec les autres méthodes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ychap=rep(coeff[1],nrow(cook.test))+as.vector(coeff[-1]%*%t(data.matrix(cook.test[,-1])))\n",
    "options(repr.plot.width=4, repr.plot.height=4)\n",
    "mean((cook.test[,1]-ychap)^2) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 Régression PLS\n",
    "\n",
    "**Q** Quel paramètre faut-il optimiser pour ce modèle?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cook.pls=plsr(sucre~.,ncomp=28,data=cook.app,validation=\"CV\",segments=cvseg)\n",
    "options(repr.plot.width=4, repr.plot.height=4)\n",
    "plot(cook.pls,pch=19)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Commenter la précédure ci-dessous et le résultat obtenu."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msepcv.pls=MSEP(cook.pls,estimate=c(\"train\",\"CV\"))\n",
    "plot(msepcv.pls,type=\"l\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ncompo=which.min(msepcv.pls$val[\"CV\",,])-1\n",
    "# modèle optimal \n",
    "cook.plso=plsr(sucre~.,ncomp=ncompo,data=cook.app)\n",
    "# résidus \n",
    "fit.pls=predict(cook.plso,ncomp=1:ncompo)[,,ncompo]\n",
    "plot(fit.pls,cook.app[,\"sucre\"], pch=19)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res.pls=fit.pls-cook.app[,\"sucre\"]\n",
    "plot.res(fit.pls,res.pls)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Commenter les deux graphiques ci-dessus.\n",
    "\n",
    "**Q** Commenter les résultats ci-dessous."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ychap=predict(cook.plso,newdata=cook.test)[,1,ncompo]\n",
    "plot(ychap,cook.test[,1])\n",
    "mean((cook.test[,\"sucre\"]-ychap)^2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Attention*, une observation (la 23ème) semble se démarquer de l'échantillon. Ce n'est pas flagrant mais elle est signalée comme atypique (*outlier*) dans la littérature. La repérer sur le graphique; elle sera prise en compte dans la section suivante.\n",
    "\n",
    "Trois graphiques analogues à ceux de l'ACP viennent compléter les résultats mais ceux-ci sont moins utiles ou pertinents qu'en version 2 de la PLS quand la variable à expliquer $Y$  est multidimensionnelle. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scoreplot(cook.pls)\n",
    "corrplot(cook.pls)\n",
    "loadingplot(cook.pls)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.4 Régression sur composantes principales\n",
    "**Q** Quel paramètre doit-il être optimisé?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cook.pcr=pcr(sucre~.,ncomp=28,data=cook.app,validation=\"CV\",segments=cvseg)\n",
    "msepcv.pcr=MSEP(cook.pcr,estimate= c(\"train\",\"CV\"))\n",
    "plot(msepcv.pcr,type=\"l\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ncompo=which.min(msepcv.pcr$val[\"CV\",,])-1\n",
    "cook.pcro=pcr(sucre~.,ncomp=ncompo,data=cook.app)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Commenter le résultat ci-dessous."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit.pcr=predict(cook.pcro,ncomp=1:6)[,,6]\n",
    "res.pcr=fit.pcr-cook.app[,\"sucre\"]\n",
    "chap=predict(cook.pcro,newdata=cook.test)[,1,ncompo]\n",
    "mean((cook.test[,\"sucre\"]-ychap)^2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.5 Random forest\n",
    "**Q** Quel paramètre serait à optimiser?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(randomForest)\n",
    "# estimation\n",
    "cook.rf=randomForest(sucre~.,data=cook.app,xtest=cook.test[,-1],ytest=\n",
    "   cook.test[,\"sucre\"],do.trace=50,mtry=50,corr.bias=TRUE)\n",
    "pred.rfr=cook.rf$test$predicted\n",
    "mean((pred.rfr-cook.test[,\"sucre\"])^2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred.rfr=cook.rf$test$predicted\n",
    "mean((pred.rfr-cook.test[,\"sucre\"])^2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Que dire de la qualité de ce modèle non linéaire"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.6 Régression Lasso par algorithme lars\n",
    "Le déroulement suit les mêmes étapes.\n",
    "\n",
    "**Q** Quel paramètre est à optimiser?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(lars)\n",
    "frac.delta=seq(0,1,length=1000)\n",
    "mse.cv=cv.lars(data.matrix(cook.app[,-1]),cook.app[,1],K=4,se=F,index=frac.delta,use.Gram=F)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "frac.delta.o=frac.delta[which.min(mse.cv$cv)]\n",
    "cook.lasso=lars(data.matrix(cook.app[,-1]),cook.app[,1],use.Gram=F)\n",
    "fit.lasso=predict(cook.lasso,data.matrix(cook.app[,-1]),s=frac.delta.o,mode=\"fraction\")$fit\n",
    "plot(fit.lasso,cook.app[,\"sucre\"], pch=19)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred.lasso=predict(cook.lasso,data.matrix(cook.test[,-1]),s=frac.delta.o,mode=\"fraction\")$fit\n",
    "mean((pred.lasso-cook.test[,\"sucre\"])^2)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Commenter le résultat précédent.\n",
    "\n",
    "**Q** Que représente le graphique ci-dessous?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(cook.lasso,breaks=F)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4 Modélisation industrielle\n",
    "Kunh (2008) a développé tout un environnement (librairie `caret`) qui facilite modélisation et optimisation des paramètres pour quasiment toutes les méthodes de classification et/ou régression disponibles dans des librairies R et, il y en a beaucoup. Chacune possède une logique ou une stratégie particulière de gestion des paramètres; cet outil permet d'en uniformiser la syntaxe. Cela concerne la transmission des données, du ou des paramètres de complexité à optimiser lors de l'appel des fonctions et le calcul des prévisions.\n",
    "\n",
    "La librairie possède plusieurs fonctionnalités qui sont mises en oeuvre et illustrées successivement dans le cas présent de la régression:\n",
    "- extraction des échantillons d'apprentissage et test,\n",
    "- transformation et sélection préliminaire des variables,\n",
    "- estimation des modèles,\n",
    "- choix du mode d'estimation de l'erreur de prévision (validation croisée, bootstrap, .632 bootstrap, \\emph{out of bag}...)\n",
    "- optimisation des paramètres (complexité, pénalisation),\n",
    "- calcul des prévisions de l'échantillon test ou d'un autre jeu de données,\n",
    "- graphes de diagnostic (résidus). \n",
    "\n",
    "Les mêmes fonctionnalités sont proposées en situation de classification supervisée avec des spécificités pour le type de l'erreur (mauvais classement, courbe ROC).\n",
    "\n",
    "\n",
    "De façon générale, le défi pour la suite, est d'arriver à faire mieux en terme de qualité de prévision que la régression PLS.\n",
    "\n",
    "### 4.1 Préparation des données\n",
    "Osbone et al. (1984) et leur successeurs s'accordent sur le fait que l'observation 23 est atypique. Ils font remarquer de plus que les extrémités des spectres sont bruitées ; en les retirant, les erreurs de prévision sont réduites. Enfin, Goutils et Fearn (1996) suggèrent également que la discrétisation est trop fine. Que devient la prévision avec quatre fois moins de variables ?\n",
    "\n",
    "Il s'agit donc de reprendre les calculs des erreurs de prévisions avec l'extraction ci-dessous qui supprime la 23ème observation, les 50 premières et dernières valeurs du spectre tout en ne considérant qu'une valeur sur 4:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "cook2=cook[-23,c(1,seq(50,650,4))]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Que font les commandes suivantes?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "inTrain = 1:39\n",
    "trainDescr=cook2[inTrain,-1]\n",
    "testDescr=cook2[-inTrain,-1]\n",
    "trainY=cook2[inTrain,1]\n",
    "testY=cook2[-inTrain,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Intérêt des commandes suivantes?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(doParallel)\n",
    "cl <- makeCluster(4)\n",
    "registerDoParallel(cl) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Que calculent les transformations suivantes? Purquoi est-ce important pour une méthode comme la régression PLS?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(caret)\n",
    "xTrans=preProcess(trainDescr)\n",
    "trainDescr=predict(xTrans,trainDescr)\n",
    "testDescr=predict(xTrans,testDescr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Que définissent les commandes suivantes?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "btControl=trainControl(method=\"boot\",number=100)\n",
    "cvControl=trainControl(method=\"cv\",number=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.3 Estimation et optimisation des modèles\n",
    "La librairie intègre beaucoup plus de méthodes mais celles sélectionnées ci-dessous semblent les plus pertinentes. \n",
    "\n",
    "**Q** Commenter les paramètres optimisés et les résultats de chacune des méthodes ci-dessous.\n",
    "\n",
    "**Q** A quoi sert la commande `set.seed`?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# régression ridge\n",
    "set.seed(2)\n",
    "ridgeFit = train(trainDescr, trainY,method = \"ridge\", tuneLength = 8,trControl = cvControl)\n",
    "ridgeFit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=4, repr.plot.height=3)\n",
    "plot(ridgeFit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# régression pls  \n",
    "set.seed(2)\n",
    "plsFit = train(trainDescr, trainY,method = \"pls\", tuneLength = 12,trControl = cvControl)\n",
    "plsFit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(plsFit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# régression lasso\n",
    "set.seed(2)\n",
    "lassoFit = train(trainDescr, trainY,method = \"lasso\", tuneLength = 8,trControl = cvControl)\n",
    "lassoFit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(lassoFit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# régression lars\n",
    "set.seed(2)\n",
    "larsFit = train(trainDescr, trainY,method = \"lars\", tuneLength = 8,trControl = cvControl)\n",
    "larsFit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(larsFit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Quelles est plus précisément la méthode ci-dessous? Que sont ses paramètres?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# elasticnet\n",
    "set.seed(2)\n",
    "enetFit = train(trainDescr, trainY,method = \"enet\", tuneLength = 8,trControl = cvControl)\n",
    "enetFit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(enetFit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# support vector machine \n",
    "set.seed(2)\n",
    "svmFit = train(trainDescr, trainY,method = \"svmLinear\",trControl = cvControl)\n",
    "svmFit\n",
    "p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.4 Prévisions, graphes et erreurs\n",
    "La librairie offre la possibilité de gérer directement une liste des modèles et donc une liste des résultats.\n",
    "\n",
    "**Q** Commenter les résultats. Quelles méthodes semblent sortir du lot?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "models=list(ridge=ridgeFit,pls=plsFit,lasso=lassoFit,lars=larsFit,elasticnet=enetFit,svm=svmFit)  \n",
    "testPred=predict(models, newdata = testDescr)\n",
    "lapply(testPred, function(x) mean((x-testY)^2))\n",
    "resPlot=extractPrediction(models, testX=testDescr, testY=testY)\n",
    "plotObsVsPred(resPlot)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5 Automatisation\n",
    "\n",
    "**Q** Quelle procédure est mise ne place ci-dessous? Pourquoi?\n",
    "\n",
    "**Q** Que contient la variable `noptim`?\n",
    "\n",
    "**Q** Quelle valeur convient-il de donner à `Niter` une fois que le programme est testé?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Attention: exécuter la fonction en annexe\n",
    "models=c(\"ridge\",\"pls\",\"lasso\",\"lars\",\"enet\",\"svmLinear\")\n",
    "noptim=c(10,10,10,10,5,5)\n",
    "# Initialiser le générateur et fixer le nombre d’itérations\n",
    "# Changer ces valeurs. Attention au temps de calcul! Être patient!\n",
    "Niter=3 ; Init=11  \n",
    "# Appel de la fonction définie en annexe\n",
    "X=cook2[,-1]; Y= cook2[,1]\n",
    "pred.cook=pred.autom(X,Y,methodes=models,N=Niter,xinit=Init,size=noptim)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calcul des risques\n",
    "prev=pred.cook$pred\n",
    "obs=pred.cook$obs\n",
    "# carrés des différences\n",
    "dif=lapply(prev,function(x)(x-obs)^2)\n",
    "# risque pour chaque échantillon\n",
    "moy=lapply(dif,function(x)apply(x,2,mean))\n",
    "options(repr.plot.width=6, repr.plot.height=4)\n",
    "boxplot(data.frame(moy))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Quel est le graphique ci-dessus. Comment l'interpréter?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lapply(moy,mean)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Conclusion finale. Choix de la méthode.\n",
    "## Annexe\n",
    "**Q** Quelle est la fonction ci-dessous?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred.autom=function(X,Y,p=1/2,methodes=c(\"knn\",\"rf\"),size=c(10,2),xinit=11,N=10,typerr=\"cv\",number=4,type=\"raw\") {\n",
    "# Fonction de prévision de N échantillons tests\n",
    "# par une liste de méthodes de régression\n",
    "# ou classification (uniquement 2 classes)\n",
    "# Optimisation des paramètres par validation\n",
    "# croisée (défaut) ou bootstrap ou... (cf. caret)\n",
    "# X : matrice ou frame des variables explicatives\n",
    "# Y : variable cible quantitative ou qualitative\n",
    "# p : proportion entre apprentissage et test\n",
    "# methodes : liste des méthodes de rdiscrimination\n",
    "# size : e grille des paramètres à optimiser\n",
    "# xinit : générateur de nombres aléatoires\n",
    "# N : nombre de réplications apprentissage/test\n",
    "# typerr : \"cv\" ou \"boo\" ou \"oob\"\n",
    "# number : nombre de répétitions CV ou bootstrap\n",
    "# pred : liste des matrices de prévision\n",
    "# type d’erreur\n",
    "Control=trainControl(method=typerr,number=number)\n",
    "# initialisation du générateur\n",
    "set.seed(xinit)\n",
    "# liste de matrices stockant les prévisions\n",
    "# une par méthode\n",
    "inTrain=createDataPartition(Y,p=p,list=FALSE)\n",
    "ntest=length(Y[-inTrain])\n",
    "pred=vector(\"list\",length(methodes))\n",
    "names(pred)=methodes\n",
    "pred=lapply(pred,function(x)x=matrix(0,nrow=ntest,ncol=N))\n",
    "obs=matrix(0,ntest,N)\n",
    "set.seed(xinit)\n",
    "for(i in 1:N) {\n",
    " # N itérations\n",
    " # indices de l’échantillon d’apprentissage\n",
    " inTrain=createDataPartition(Y,p=p,list=FALSE)\n",
    " # Extraction des échantillons\n",
    " trainDescr=X[inTrain,]\n",
    " testDescr=X[-inTrain,]\n",
    " trainY=Y[inTrain]\n",
    " testY=Y[-inTrain]\n",
    " # stockage des observés de testY\n",
    " obs[,i]=testY\n",
    " # centrage et réduction des variables\n",
    " xTrans=preProcess(trainDescr)\n",
    " trainDescr=predict(xTrans,trainDescr)\n",
    " testDescr=predict(xTrans,testDescr)\n",
    " # estimation et optimisation des modèles\n",
    " # pour chaque méthode de la liste\n",
    " for(j in 1:length(methodes)) {\n",
    "  # modélisation\n",
    "  modFit = train(trainDescr, trainY,method = methodes[j], tuneLength = size[j],trControl = Control)\n",
    "  # prévisions\n",
    "  if (type==\"prob\")  pred[[j]][,i]=predict(modFit,newdata = testDescr,type=type)[,1]\n",
    "  else pred[[j]][,i]=predict(modFit,\n",
    "  newdata = testDescr)\n",
    "  }}\n",
    "list(pred=pred,obs=obs)\n",
    "# résultats\n",
    "}"
   ]
  }
 ],
 "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}