NIR/Apprent-R-tecator.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 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",
    "## Tecator\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Résumé\n",
    "Modélisation (ou calibration) en grande dimension. Les variables explicatives sont les discrétisations de spectres en proche infra rouge (NIR) et la variable à expliquer la part de matière grasse dans des échantillons de viande. Différente méthodes de régression [lasso, krls](http://wikistat.fr/pdf/st-m-app-select.pdf), [pls](http://wikistat.fr/pdf/st-m-app-sparse-pls.pdf), [svm](http://wikistat.fr/pdf/st-m-app-svm.pdf), [mars](http://wikistat.fr/pdf/st-m-app-cart.pdf)...) sont utilisées et leur qualités prédictives comparées. La prise en compte de la forme des spectres par l'introduction des dérivées après lissage spline est utile sur ce jeu de données. La librairie `caret` est directement utilisée pour industrialiser la stratégie de choix de modèle et de méthode. Ce scénario vient compléter celui des données NIR de pâte à biscuit."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1 Introduction\n",
    "\n",
    "### 1.1 Contexte\n",
    "*These data are recorded on a Tecator Infratec Food and Feed Analyzer working in the wavelength range 850 - 1050 nm by the Near Infrared Transmission (NIT) principle. Each sample contains finely chopped pure meat with different moisture, fat and protein contents.}\\footnote{The data are available in the public domain with no responsability from the original data source. The data can be redistributed as long as this permission note is attached.*\n",
    "\n",
    "Chaque échantillon de viande est caractérisé par des mesures d'absorption pour 100 longueurs d'ondes  échantillonnées dans le proche infra rouge ainsi que par des mesures chimiques de la proportion (pourcentage) en eau, graisse et protéine. Les mesures d'absorption ont déjà été transformées ($\\log_{10}$). Nous nous attacherons à modéliser le pourcentage de matière grasse.\n",
    "\n",
    "Ces données ont été largement étudiées dans la littérature car elles servent de jeu de donnée test (*benchmark*) pour comparer des prétraitements des données (centrage, réduction sur données brutes ou après ACP, pls, décomposition sur base de splines et dérivations, réduction de dimension par ACP ou PLS) des méthodes de modélisation (régression linéaire, svm, réseaux de neurones, cart...), des techniques de sélection de variables (stepwise, par critère bayésien, informaiton mutuelle, par algorithme stochastique...) et toutes les combinaisons de ces différentes approches...\n",
    "\n",
    "Contrairement à la plupart des données NIR du domaine public, ces données présentent une *composante non linéaire* certes légère mais présente. De plus, la *prise en compte des dérivées* ou tout du moins des différences d'ordre 2 des spectres est aussi pertinente, d'où l'intérêt qu'elles ont sucité et suscitent toujours dans la littérature du métier, mais aussi celle statistique ou de *machine learning*, pour justifier l'utilisation de nouvelles méthodes: non-linéaires et/ou fonctionnelles.\n",
    "\n",
    "### 1.2 Données\n",
    "Les données  sont  accessibles dans R au sein du package `caret`. Elles sont contenues dans deux matrices: `absorp` de $n=215$ lignes et $p=100$ colonnes, les absorptions par longueur d'onde, et  `endpoints` de $n=215$ lignes et $p=3$ colonnes, les pourcentages d'eau, graisses, et protéines. Les 129 premières observations sont généralement utilisées comme échantillon d'apprentissage.\n",
    "\n",
    "\n",
    "### 1.3 Objectif\n",
    "L'objectif principal est la recherche d'un meilleur modèle de prévision du pourcentage de matière grasse. Il s'agit d'évaluer, comparer, différentes stratégies et méthodes pour aboutir à celle la plus efficace. Le travail est découpé en trois parties. La première explore les spectres, leur lissage et de leur dérivées; la deuxième modélise les spectres initiaux tandis que la suivante prend en compte les dérivées des spectre après lissage spline avant finalement d'automatiser les traitements pour optimiser les choix.\n",
    "\n",
    "## 2 Exploration\n",
    "### 2.1 Statistiques élémentaires\n",
    "Trois variables très coorélées entre elles peuvent être modélisées et prédites à partir des mesutes d'absorbtion ou spectres. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(caret)\n",
    "\n",
    "# Variables à expliquer\n",
    "data(tecator)\n",
    "summary(endpoints)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Spectre des absorptions ou variables explicatives\n",
    "dim(absorp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# corrélation des variables à expliquer\n",
    "options(repr.plot.width=4, repr.plot.height=4)\n",
    "splom(endpoints)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## tracer un échantillon aléatoire de 10 spectres\n",
    "# sélection des observations\n",
    "set.seed(1)\n",
    "inSubset = sample(1:dim(endpoints)[1], 10)\n",
    "absorpSubset = absorp[inSubset,]\n",
    "endpointSubset = endpoints[inSubset, 2]\n",
    "# tri des spectres selon la première valeure\n",
    "newOrder =  order(absorpSubset[,1])\n",
    "absorpSubset =  absorpSubset[newOrder,]\n",
    "endpointSubset =  endpointSubset[newOrder]\n",
    "# définition des couleurs\n",
    "plotColors =  rainbow(10)\n",
    "# tracé des échelles\n",
    "plot(absorpSubset[1,], type = \"n\", ylim = range(absorpSubset), xlim = c(0, 105), xlab = \"Wavelength Index\", ylab = \"Absorption\")\n",
    "# tracé des spectres et du taux de graisse\n",
    "for(i in 1:10){\n",
    "points(absorpSubset[i,],type = \"l\", \n",
    "   col = plotColors[i], lwd = 2)\n",
    "text(105,absorpSubset[i,100], endpointSubset[i], \n",
    "   col = plotColors[i])\n",
    "}\n",
    "title(\"Profiles for 10 Random Samples\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La moyenne de chaque spectre ne semble pas être une information pertinente pour l'estimation du taux de graisse. En revanche, ils présentent des formes remarquable dans certaines zones, d'où l'intérêt éventuel d'étudier la dérivée de ces courbes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# la variable à expliquer choisie, taux de matières grasses, est la deuxième\n",
    "gras=endpoints[,2]\n",
    "hist(gras)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# les diagrammes boîtes des variables explicatives\n",
    "options(repr.plot.width=8, repr.plot.height=4)\n",
    "boxplot(data.frame(absorp))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Analyse en composantes principales"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "acp=prcomp(absorp)\n",
    "options(repr.plot.width=3, repr.plot.height=3)\n",
    "plot(acp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=4, repr.plot.height=4)\n",
    "biplot(acp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(acp$x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot.ts(acp$rotation[,1:10])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Très grande importance de la première composante principale. Elle a pour forme celle moyenne des courbes dont les décalages génèrent une variance importante. \n",
    "\n",
    "### 2.3 Lissage spline et dérivées des spectres\n",
    "\n",
    "Choix du paramètre de lissage ou degré de liberté qui est la trace de la *hat matrix* ou matrice de lissage. Le premier spectre, son lissage et ses dérivées sont tracés pour différentes valeurs du paramètre {\\tt df}. Si celui-ci n'est pas fourni (dernier cas), il est optimisé par validation croisée afin de reconstruire au mieux la courbe. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(splines)\n",
    "options(repr.plot.width=6, repr.plot.height=6)\n",
    "par(mfcol=c(3,3))\n",
    "degre=c(3,10,20)\n",
    "for (i in 1:3){\n",
    "model=smooth.spline(x=1:100,y=absorp[1,],df=degre[i])\n",
    "lisse=predict(model)$y\n",
    "deriv1=predict(model,deriv=1)$y\n",
    "deriv2=predict(model,deriv=2)$y\n",
    "ts.plot(cbind(lisse,absorp[1,]))\n",
    "ts.plot(deriv1)\n",
    "ts.plot(deriv2)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model=smooth.spline(x=1:100,y=absorp[1,])\n",
    "lisse=predict(model)$y\n",
    "deriv1=predict(model,deriv=1)$y\n",
    "deriv2=predict(model,deriv=2)$y\n",
    "options(repr.plot.width=6, repr.plot.height=3)\n",
    "par(mfcol=c(1,3))\n",
    "ts.plot(cbind(lisse,absorp[1,]))\n",
    "ts.plot(deriv1)\n",
    "ts.plot(deriv2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Il semblerait raisonnable de lisser un peu plus que ce que suggère la validation croisée afin d'obtenir des dérivées plus lisses. L'idéal serait évidemment d'optimiser ce paramètre en cherchant à prévoir au mieux le pourcentage de matières grasses.\n",
    "\n",
    "Génération des bases de B-splines cubiques et tracer de quelques unes d'entre elles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "base=bs(1:100,df=model$fit$nk)\n",
    "options(repr.plot.width=3, repr.plot.height=3)\n",
    "ts.plot(base[,c(1,20,30,64)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour chaque spectre, les commandes ci-dessous calculent un lissage des spectres par B-splines cubiques : noeuds équirépartis, degré fixé à 15. Le programme fournit les dérivées première et seconde aux points de discrétisation, les coefficients dans la base de base de B-Splines. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "degre=15\n",
    "# initialisation des matrices\n",
    "coeff=matrix(ncol=model$fit$nk,nrow=215)\n",
    "lisse=matrix(ncol=100,nrow=215)\n",
    "deriv1=matrix(ncol=100,nrow=215)\n",
    "deriv2=matrix(ncol=100,nrow=215)\n",
    "\n",
    "# itération sur tous les spectres\n",
    "for (obs in 1:215){\n",
    "   model=smooth.spline(x=1:100,y=absorp[obs,],df=degre)\n",
    "   coeff[obs,]=model$fit$coef\n",
    "   lisse[obs,]=predict(model)$y\n",
    "   deriv1[obs,]=predict(model,deriv=1)$y\n",
    "   deriv2[obs,]=predict(model,deriv=2)$y\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## 3 Modélisation des spectres\n",
    "L'étude se limite à une utilisation de la librairie `caret` (Kunh, 2008) pour quelques unes des méthodes de modélisation proposées. \n",
    "\n",
    "### 3.1 Préparation des données\n",
    "Extraction des échantillons d'apprentissage et de test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dimnames(absorp)[[2]]=1:100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Indices usuelles des échantillons \n",
    "# d'apprentissage et de test\n",
    "inTrain = 1:129\n",
    "trainDescr=absorp[inTrain,]\n",
    "testDescr=absorp[-inTrain,]\n",
    "trainY=gras[inTrain]\n",
    "testY=gras[-inTrain]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Il est recommandé de centrer et réduire les variables dans plusieurs méthodes. C'est fait systématiquement et simplement en utilisant évidemment les mêmes transformations sur l'échantillon test que celles mises en place sur l'échantillon d'apprentissage. D'autres part, l'erreur de prévision est estimée par validation croisée."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xTrans=preProcess(trainDescr)\n",
    "trainDescr=predict(xTrans,trainDescr)\n",
    "testDescr=predict(xTrans,testDescr)\n",
    "cvControl=trainControl(method=\"cv\",number=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Estimation et optimisation des modèles\n",
    "La librairie intègre beaucoup plus de méthodes mais celles séelctionnées ci-dessous semblent les plus pertinentes. Exécuter successivement les blocs de commandes pour tracer séparamment chacun des graphes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(doParallel)\n",
    "cl <- makeCluster(4)\n",
    "registerDoParallel(cl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#1 Régression linéaire leapBackward\n",
    "set.seed(2)\n",
    "llmFit = train(trainDescr, trainY,method = \"leapBackward\", tuneLength = 20,trControl = cvControl)\n",
    "plot(llmFit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#2 régression pls  \n",
    "set.seed(2)\n",
    "plsFit = train(trainDescr, trainY,method = \"pls\", tuneLength = 30,trControl = cvControl)\n",
    "plot(plsFit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#3 Ridge regression\n",
    "set.seed(2)\n",
    "ridgeFit = train(trainDescr, trainY,method = \"ridge\", tuneLength = 8,trControl = cvControl)\n",
    "plot(ridgeFit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#4 Regression elastic net\n",
    "set.seed(2)\n",
    "enetFit = train(trainDescr, trainY,method = \"enet\", tuneLength = 8,trControl = cvControl)\n",
    "plot(enetFit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#5 Support vector machine noyau linéaire\n",
    "set.seed(2)\n",
    "svmFit = train(trainDescr, trainY, method = \"svmLinear\",trControl = cvControl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#6 Support vector machine noyau gaussien\n",
    "set.seed(2)\n",
    "svmgFit = train(trainDescr, trainY,method = \"svmRadial\", tuneLength = 8,trControl = cvControl)\n",
    "plot(svmgFit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#8 Kernel-based Regularized Least Squares (KRLS) \n",
    "set.seed(2)\n",
    "krlsFit = train(trainDescr, trainY,method = \"krlsRadial\", tuneLength = 10,trControl = cvControl)\n",
    "plot(krlsFit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "models=list(leaplm=llmFit,pls=plsFit,elasticnet=enetFit,svm=svmFit,svmg=svmgFit,krls=krlsFit)  \n",
    "testPred=predict(models, newdata = testDescr)\n",
    "lapply(testPred,function(x)sqrt(mean((x-testY)^2)))\n",
    "resPlot=extractPrediction(models, testX=testDescr, testY=testY)\n",
    "plotObsVsPred(resPlot)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remarquer quels sont les modèle prenant en compte la non linéarité des données, ceux acceptant des valeurs tests relativements atypiques."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4 Utilisation des dérivées\n",
    "### Extraction des échantillons\n",
    "Extraction des échantillons d'apprentissage et de test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mêmes indices de l'échantillon d'apprentissage\n",
    "dimnames(deriv1)[[2]]=1:100\n",
    "inTrain = 1:129\n",
    "trainDescr1=deriv1[inTrain, ]\n",
    "testDescr1 =deriv1[-inTrain,]\n",
    "# pré-process\n",
    "xTrans1=preProcess(trainDescr1)\n",
    "trainDescr1=predict(xTrans1,trainDescr1)\n",
    "testDescr1 =predict(xTrans1,testDescr1 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Le calcul est restreint au 3 meilleures méthodes trouvées précédemment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Régression PLS\n",
    "set.seed(2)\n",
    "plsFit1 = train(trainDescr1, trainY,method = \"pls\", tuneLength = 30,trControl = cvControl)\n",
    "plot(plsFit1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Kernel-based Regularized Least Squares (KRLS) \n",
    "set.seed(2)\n",
    "krlsFit1 = train(trainDescr1, trainY,method = \"krlsRadial\", tuneLength = 10,trControl = cvControl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(krlsFit1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Caclul des erreurs et tracé des graphes\n",
    "models=list(pls1=plsFit1,krls1=krlsFit1)  \n",
    "testPred1=predict(models, newdata = testDescr1)\n",
    "lapply(testPred1,function(x)sqrt(mean((x-testY)^2)))\n",
    "resPlot=extractPrediction(models, testX=testDescr1, testY=testY)\n",
    "plotObsVsPred(resPlot)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2 Automatisation de la  validation croisée *Monte Carlo*\n",
    "L'échantillon est de faible taille, les estimations des erreurs très dépendantes de l'échantillon test sont sujettes à caution et on peut s'interroger sur la réalité des différences entre les différentes méthodes. En regardant les graphiques, on observe qu'il suffit d'une observation pour influencer les résultats. Il est donc important d'itérer le processus sur plusieurs échantillons tests afin de décrire les distributions des erreurs.\n",
    "\n",
    "Evidemment le temps de calcul s'en ressent.\n",
    "\n",
    "Exécuter la fonction en annexe avant d'exécuter: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "models=c(\"pls\",\"krlsRadial\")\n",
    "noptim=c(10,10,10,10)\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=30 ; Init=11  \n",
    "# Appel de la fonction définie en annexe\n",
    "X=deriv1; Y= gras\n",
    "pred.tecat=pred.autom(X,Y,methodes=models,N=Niter,xinit=Init,size=noptim)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Puis calculer :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calcul des risques\n",
    "prev=pred.tecat$pred\n",
    "obs=pred.tecat$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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lapply(moy,mean)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Annexe: Programme de validation croisée *Monte Carlo*"
   ]
  },
  {
   "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",
    "}"
   ]
  },
  {
   "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}