StatElem-R-MLG.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'Initiation à la Statistique](https://github.com/wikistat/Exploration)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modèle Linéaire Général  (MLG), exemples en <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",
    "\n",
    "**Résumé** Ce scénario présente trois exemples de modèle linéaire général. Le premier est une [analyse de variance / covariance](http://wikistat.fr/pdf/st-m-modlin-anacova.pdf) (modèle gaussien) pour une enquête marketing. Les deux suivants des exemples de [régression logistique](http://wikistat.fr/pdf/st-m-modlin-logit.pdf) (modèle binomial). le deuxième exemple analyse une enquête sur le travail (absence / présence) des femmes, le troisième propose de modéliser puis prévoir une variable clinique binaire: occurrence d'une pathologie cardiaque, par un ensemble de variables qualitatives et quantitatives. C'est l'illustration de différentes stratégies pour optimiser une [régression logistique](http://wikistat.fr/pdf/st-m-modlin-reglog.pdf) (modèle binomial) avec sélection de variables avec  minimisation du critère AIC, estimation de l'erreur de prévision, courbe ROC.\n",
    "\n",
    "\n",
    "**Répondre aux questions en s'aidant des résultats des exécutions**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1 Modèle gaussien: AnOVa *vs.* AnCoVa\n",
    "L'étude de ces données marketing est déjà réalisée dans un [scénario avec SAS](http://www.math.univ-toulouse.fr/~besse/Wikistat/pdf/st-scenar-reg-GLM-SAS.pdf). Elle n'est pas complètement détaillée ici, seule les modèles complets sont repris à titre de comparaison. \n",
    "\n",
    "### 1.1 Données d'une enquête marketing\n",
    "Les données [milkR.dat](https://github.com/wikistat/StatElem/blob/master/milkR.dat), extraites de Jobson (1991), sont issues d'une étude marketing visant à étudier l'impact de campagnes publicitaires sur les ventes de différents aliments. Un échantillon ou \"panel\" de familles a été constitué en tenant compte du lieu d'habitation ainsi que de la constitution de la famille. Chaque semaine, chacune de ces familles ont rempli un questionnaire décrivant les achats réalisés. Nous nous limitons ici à l'étude de l'impact sur la *consommation de lait* de quatre campagnes diffusées sur des chaînes locales de télévision. Quatre villes, une par campagne publicitaire, ont été choisies dans cinq différentes régions géographiques. Les consommations en lait par chacune des six familles par ville alors été mesurées (en dollars) après deux mois de campagne. \n",
    "\n",
    "Les données initiales se présentent sous la forme d'un tableau à 6 variables~:\n",
    "la région géographique, les 4 consommations pour chacune des  villes ou campagnes publicitaires diffusées, la taille de la famille. \n",
    "\n",
    "Elles sont lues et réorganisées pour satisfaire à la structure classique d'un plan factoriel équilibré (donc complet) croisant trois facteurs~: région, taille de la famille, type de campagne. L'objectif est d'étudier l'effet du type de publicité sur la consommation des ménages. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lecture des données\n",
    "dat1=read.table(\"milkR.dat\",header=T,sep=\",\")\n",
    "summary(dat1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " Mise en forme des données par construction de la matrice du plan d'expérience (*design matrix*). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mise en forme\n",
    "design = expand.grid(\n",
    "    Pub=factor(1:4, labels=c(\"Pub1\",\"Pub2\",\"Pub3\",\"Pub4\")),\n",
    "    Taille=factor(1:6),\n",
    "    Region=factor(1:5, labels=c(\"R1\",\"R2\",\"R3\",\"R4\",\"R5\")))\n",
    "# vectorisation des données de consommation\n",
    "# La taille est un facteur\n",
    "Conso=as.vector(t(dat1[,2:5]))\n",
    "dat2=data.frame(design,Conso)\n",
    "summary(dat2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Même chose mais en considérant que la variable \"Taille\" (de la famille) est une variable quantitative."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# La taille est quantitative\n",
    "dat3=data.frame(dat2[,-2],\"Taille\"=as.numeric(dat2[,2]))\n",
    "summary(dat3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualisation des données sous cette forme. \n",
    "\n",
    "**Q** Est-il possible de conclure à des effets des facteurs sur la consommation?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "car=as.integer(dat3[,\"Region\"])\n",
    "plot(dat3[,\"Taille\"],dat3[,\"Conso\"],col=dat3[,\"Pub\"],xlab=\"Taille\",ylab=\"Consommation\",pch=car)\n",
    "legend(\"topleft\",pch=1:5,legend=c(\"Reg1\",\"Reg2\",\"Reg3\",\"Reg4\",\"Reg5\"))\n",
    "legend(\"bottomright\",legend=c(\"Pub1\",\"Pub2\",\"Pub3\",\"Pub4\"),text.col=1:4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** et de choisir une campagne de publicié \"optimale\" par région?\n",
    "### 1.3 Modèle d'AnOVa\n",
    "La variable \"taille\" est considérée qualitative comme un facteur à 6 nivaux. L'approche rigoureuse consiste à prendre en compte directement toutes les variables du plan d'expérience. \n",
    "\n",
    "Une version simplifiée sans interaction:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod1=glm(Conso~Taille+Pub+Region,data=dat2,family=gaussian)\n",
    "print(mod1)\n",
    "anova(mod1,test=\"F\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Une version complète avec les interactions d'ordre 2; celle d'ordre 3 ne peut pas être estimée car les observations ne sont pas répliquées."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod2=glm(Conso~(Pub+Taille+Region)**2,data=dat2,family=gaussian)\n",
    "print(mod2)\n",
    "anova(mod2,test=\"F\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Quel modèle retenir ? \n",
    "\n",
    "**Q** Pouvez vous conclure sur le choix de la campagne de publicité optimale par région?\n",
    "\n",
    "### 1.4 Modèle d'AnCoVa complet\n",
    "Même chose en considérant la variable \"taille\" comme quantitative. \n",
    "\n",
    "**Q** Quel avantage à prendre cette option?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod3=glm(Conso~(Pub+Taille+Region)**2,data=dat3,family=gaussian)\n",
    "print(mod3)\n",
    "anova(mod3,test=\"F\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Comment choisir la campagne de publicité optimale? retrouver les ordonnées à l'origine et les pentes de chaque droite de régression par couple pub x région.\n",
    "\n",
    "**Q** Comparer avec les résultats (p-valeurs) de l'anova."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2 Modèle binomial: régression logistique de données d'enquête\n",
    "### 2.1 Les données \n",
    "Les données étudiées (Jobson 1992) sont issues d’une enquête réalisée auprès de 200 femmes mariées du Michigan. Les variables qualitatives sont les suivantes:\n",
    "- `THISYR`: la variable à expliquer, (`Woui`) si la femme travaille l'année en cours (`Wnon`) sinon;\n",
    "- `LASTYR`: (`Loui`) si la femme travaille l'année précédente (`Lnon`) sinon;\n",
    "- `CHILD1`: code la présence (`Boui`) ou l'absence (`Bnon`) d’un enfant de moins de 2 ans;\n",
    "- `CHILD2`: présence (`Eoui`) ou absence (`Enon`) d’un enfant entre 2 et 6 ans;\n",
    "- `ASCEND`: ascendance noire (`Anoi`) ou blanche  (`Abla`).\n",
    "\n",
    "Les  autres  variables,  âge (`AGE`),  nombre  d’années  d’études (`EDUC`), revenu du mari (`HUBINC`) sont quantitatives.\n",
    "\n",
    "**Attention** les données ne sont pas issues d'un plan expérimental orthogonal. Elles sont disponibles dans le fichier [jobpanel.dat](https://github.com/wikistat/StatElem/blob/master/jobpanel.dat).\n",
    "\n",
    "Lire le fichier puis recoder les facteurs. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lecture des données:\n",
    "type=c(\"character\",\"character\",\"numeric\",\"numeric\",\"numeric\",\"character\",\"character\",\"character\")\n",
    "panel=read.table(\"jobpanel.dat\",colClasses=type,header=TRUE)\n",
    "summary(panel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Codage explicite des facteurs\n",
    "panel[,\"THISYR\"]=factor(panel[,1],levels=c(\"0\",\"1\"),labels=c(\"Wnon\",\"Woui\"))\n",
    "panel[,\"LASTYR\"]=factor(panel[,2],levels=c(\"0\",\"1\"),labels=c(\"Lnon\",\"Loui\"))\n",
    "panel[,\"HUBINC\"]=panel[,3]\n",
    "panel[,\"AGE\"]=panel[,4]\n",
    "panel[,\"EDUC\"]=panel[,5]\n",
    "panel[,\"CHILD1\"]=factor(panel[,6],levels=c(\"0\",\"1\"),labels=c(\"Bnon\",\"Boui\"))\n",
    "panel[,\"CHILD2\"]=factor(panel[,7],levels=c(\"0\",\"1\"),labels=c(\"Enon\",\"Eoui\"))\n",
    "panel[,\"ASCEND\"]=factor(panel[,8],levels=c(\"0\",\"1\"),labels=c(\"Abla\",\"Anoi\"))\n",
    "panel=panel[,-c(1:8)]\n",
    "summary(panel)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Exploration\n",
    "#### Unidimensionnelle\n",
    "Vérifier les distributions des variables quantitatives, justifier la transformation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hist(panel[,\"HUBINC\"])\n",
    "panel[,\"LHUBINC\"]=log(1+panel[,\"HUBINC\"])\n",
    "hist(panel[,\"LHUBINC\"])\n",
    "hist(panel[,\"AGE\"])\n",
    "hist(panel[,\"EDUC\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Bidimensionnelle\n",
    "Observer les liaisons entre les variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(panel[,\"AGE\"],panel[,\"LHUBINC\"])\n",
    "mosaicplot(THISYR~ASCEND,data=panel)\n",
    "boxplot(AGE~THISYR,data=panel)\n",
    "boxplot(LHUBINC~THISYR,data=panel)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Multidimensionnelle\n",
    "Transformer  les  variables  quantitatives  en  qualitatives  avant  de  calculer l’analyse multiple des correspondances à l’aide de la librairie FactoMineR."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Transformations\n",
    "panel[,\"AGEq\"]=cut(panel$AGE,breaks=quantile(panel[,\"AGE\"],probs = seq(0, 1, 1/3)),\n",
    "               labels=c(\"Ag1\",\"Ag2\",\"Ag3\"),include.lowest = TRUE)\n",
    "panel[,\"HUBINCq\"]=cut(panel$HUBINC,breaks=quantile(panel[,\"HUBINC\"],probs=seq(0,1,1/3)),\n",
    "               labels=c(\"HI1\",\"HI2\",\"HI3\"),include.lowest=TRUE)\n",
    "panel[,\"EDUCq\"]=cut(panel$EDUC,breaks=c(0,11,13,20),labels=c(\"Edu0\",\"Edu1\",\"Edu2\"),include.lowest=TRUE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "summary(panel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# AFCM \n",
    "library(FactoMineR)\n",
    "afcm=MCA(panel[,c(1,2,6,7,8,10,11,12)],quali.sup=c(1,2),graph=F)\n",
    "plot(afcm,habillage=\"quali\",invisible=\"ind\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tenter une interprétation des axes.\n",
    "### 2.3 Modélisation par régression logistique\n",
    "La  prévision  de  la  variable `THISYR` n’est pas l’objectif, celle-ci est d'ailleurs très mauvaise; l’objectif est plutôt l'explication de cette variable par les autres afin de mettre en évidence l’influence des différents facteurs.\n",
    "\n",
    "Par ailleurs la variable `LASTYR` très liée à celle à expliquer n'est pas utilisée dans le smodèles pour se focaliser sur les factuers socio-économiques.\n",
    "\n",
    "Le modèle complet sans interaction est estimé en utilisant soit les variables quantitatives, soit celles quantitatives transformées en qualitatives par découpage en classes.\n",
    "\n",
    "#### Avec les variables quantitatives"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "panel.glm=glm(THISYR~LHUBINC+AGE+EDUC+CHILD1+CHILD2+ASCEND,family=binomial,data=panel)\n",
    "summary(panel.glm)\n",
    "anova(panel.glm,test=\"Chi\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Que dire de la qualité d'ajustement du modèle?\n",
    "\n",
    "**Q** Il ne s'agit pas d'un plan d'expérience orthogonal. Quelle conséquence pour les p-valeurs et l'interprétation des effets?\n",
    "\n",
    "**Q** Interpréter l'influence des facteurs.\n",
    "\n",
    "#### Avec tout qualitatif"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "panel.glm=glm(THISYR~HUBINCq+AGEq+EDUCq+CHILD1+CHILD2+ASCEND,family=binomial,data=panel)\n",
    "summary(panel.glm)\n",
    "anova(panel.glm,test=\"Chi\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** L'ajustement est-il meilleur?\n",
    "\n",
    "**Q** L'interprétation de l'influence des facteurs reste-t-elle la même?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "#### Avec interactions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Supprimer les warnings pour simplifier la sortie\n",
    "options(warn=-1)\n",
    "panel.glm=glm(THISYR~(LHUBINC+AGE+EDUC+CHILD1+CHILD2+ASCEND)**2,family=binomial,data=panel)\n",
    "panel.step=step(panel.glm,direction=\"backward\", trace=0)\n",
    "options(warn=0)\n",
    "summary(panel.step)\n",
    "anova(panel.step,test=\"Chi\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Quel problème (*warnings*) apparaît?\n",
    "\n",
    "**Q** Que devient l'ajustement du modèle?\n",
    "\n",
    "*Remarque*: Si l'ojectif est l'explication plutôt que la prévision, il serait utile de continuer à éliminer des termes du modèles (un par un) car beaucoup restent \"approximativement\" non significatifs et cela éviterait les *warnings*. Cette option utilisant le test de Ward ou celui du maximum de vraisemblance, plutôt que le critère AIC, est présente dans SAS, pas dans R. Il faudrait le faire \"à la main\".\n",
    "\n",
    "**Q** Une interaction reste présente entre la variable `CHILD1` et celle `ASCEND`. Comment l'interpréter?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3 Modèle binomial : régression logistique pour prévision d'une coronopathie\n",
    "### 3.1 Gestion des données\n",
    "#### Données cliniques\n",
    "Des données publiques disponibles sur le site [UCI repository](http://archive.ics.uci.edu/ml/) décrivent des facteurs de risque et résultats cliniques: 13 parmi 75 de l’étude originale de [Detrano et al. (1989)](http://www.ajconline.org/article/0002-9149%2889%2990524-9/abstract), liés à une maladie coronarienne (athérosclérose). Celle-ci est jugée présente lorsque tous les vaisseaux coronariens sont obstrués à plus de 50% par des athéromes. Les variables étudiées sont observées sur un échantillon de 270 patients admis dans une clinique de Cleveland (Ohio) à la suite de douleurs thoraciques pouvant être dues à une angine de poitrine. Elles sont décrites dans le tableau ci-dessous:\n",
    "\n",
    "Num| Code | Libellé | Valeurs\n",
    "-|--|--|--\n",
    "1|`Age`   |  \t\t\n",
    "2|`Sexe`|   |sxF, sxM\n",
    "3|`Douleur`|\tThoracique|\tdlA (angine typique), dlb(atypique) dlc(différent) dlD(asymptom.)\n",
    "4|`Tension`|\tSystolique|\tmmHg à l’admission et au repos\n",
    "5|`Cholest`|\tTaux|\tmg/dl (préférable<200, limite entre 200 et 240, risqué au-delà)\n",
    "6|`Sucre`|\tTaux à jeun|\tscN (<120mg/dl), scO (>120mg/dl)\n",
    "7|`Cardio`|\tECG au repos|\tcdA (Normal) cdB (ST/T anormal) cdC (hypertrophie ventr. gauche)\n",
    "8|`FreqM`|\tFréquence|\tcardiaque maximum lors du test d’effort\n",
    "9|`AngInd`|\tAngine induite|\tpar l’effort :  tmA (oui), tmB (non)\n",
    "10|`PicInd`|\tDépression ST|\tInduite par effort / repos\n",
    "11|`PentInd`|\tSegment ST| \tInduit à l’effort piA(ascendante), piB(plate), piC(descendante)  \n",
    "12|`Nvais`|\tNombre de|\tvaisseaux fl0, fl1, fl2, fl3 majeurs colorés par fluoroscopie\n",
    "13|`Thal`|\tScintigraphie|\tthN(normal) thF(défaut fixé) thR(défaut révers.) avec effort\n",
    "14|`Coro`|\tCoronaropathie|\tCoA(Absence), CoP(Présence)\n",
    "\n",
    "Certaines sont associées à des risques potentiels et d’autres sont résultats d’examens cliniques au repos ou à la suite d’un test d’effort. Les variables 1, 4, 5, 8, 10 sont quantitatives, les autres sont qualitatives dont certaines binaires : 2, 6, 9, 14. Le diagnostic (variable `Coro`) a été établi par une *angiographie* permettant de mesurer l’obstruction des artères coronariennes. \n",
    "\n",
    "Après la [phase exploratoire](https://github.com/wikistat/Exploration/blob/master/Diag-coro/Explo-R-DiagCoro.ipynb) l’objectif sur ces données est de construire un modèle de prévision de la variable `Coro` à partir de l’observation des autres, pas ou peu invasives, car l’angiographie est un **examen invasif** comportant des risques. C'est en fait plus qu'un examen puisqu'il permet par exemple de simultanément poser des *stents* pour faciliter le débit sanguin. \n",
    "\n",
    "\n",
    "*N.B.* Ce tutoriel est une version plus courte, limité à la régression logistique, de [celui](https://github.com/wikistat/Apprentissage/blob/master/Diag-coro/Apprent-R-DiagCoro.ipynb) qui propose une comparaison systématique de différentes méthodes et algorithmes d'apprentissage machine."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Préparation des données\n",
    "Cette étape résume rapidement la suite des commandes qui découlent de la phase [phase exploratoire](https://github.com/wikistat/Exploration/blob/master/Diag-coro/Explo-R-DiagCoro.ipynb) des données qu'il est vivement conseillé d'avoir pratiquée pour se familiariser avec celles-ci.\n",
    "\n",
    "Lecture des données et préparation pour construire le *dataFrame*. Les données sont lues et les codes des modalités (niveaux des facteurs) sont renommés de façon explicite."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path=\"http://www.math.univ-toulouse.fr/~besse/Wikistat/data/\"\n",
    "#path=\"\"\n",
    "heart=read.table(paste(path,\"heart.dat\",sep=\"\"))\n",
    "# recodage des classes et nom des variables\n",
    "heart=data.frame(\n",
    "    Age=heart[,1],\n",
    "    Sexe=factor(as.factor(heart[,2]),labels=c(\"sxF\",\"sxM\")),\n",
    "    Douleur=factor(as.factor(heart[,3]),labels=c(\"dlA\",\"dlB\",\"dlC\",\"dlD\")),\n",
    "    Tension=heart[,4],\n",
    "    Cholest=heart[,5],\n",
    "    Sucre=factor(as.factor(heart[,6]),labels=c(\"scN\",\"scO\")),\n",
    "    Cardio=factor(as.factor(heart[,7]),labels=c(\"cdA\",\"cdB\",\"cdC\")),\n",
    "    FreqM=heart[,8],\n",
    "    AngInd=factor(as.factor(heart[,9]),labels=c(\"tmA\",\"tmB\")),\n",
    "    PicInd=heart[,10],\n",
    "    PenteInd=factor(as.factor(heart[,11]),labels=c(\"piA\",\"piB\",\"piC\")),\n",
    "    Nvais=factor(as.factor(heart[,12]),labels=c(\"fl0\",\"fl1\",\"fl2\",\"fl3\")),\n",
    "    Thal=factor(as.factor(heart[,13]),labels=c(\"thN\",\"thF\",\"thR\")),\n",
    "    Coro=factor(as.factor(heart[,14]),labels=c(\"CoA\",\"CoP\")))\n",
    "summary(heart)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Découpage des variables quantitatives en trois classes pour permettre une analyse factorielle multiple des correspondances."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "AgeQ=cut(heart$Age,breaks=quantile(heart$Age,c(0,.33,.66,1)),labels=c(\"AgeA\",\"AgeB\",\"AgeC\"),include.lowest = TRUE)\n",
    "TensQ=cut(heart$Tension,breaks=quantile(heart$Tension,c(0,.33,.66,1)),labels=c(\"TensA\",\"TensB\",\"TensC\"),include.lowest = TRUE)\n",
    "CholQ=cut(heart$Cholest,breaks=quantile(heart$Cholest,c(0,.33,.66,1)),labels=c(\"CholA\",\"CholB\",\"CholC\"),include.lowest = TRUE)\n",
    "FreQ=cut(heart$FreqM,breaks=quantile(heart$FreqM,c(0,.33,.66,1)),labels=c(\"FreqA\",\"FreqB\",\"FreqC\"),include.lowest = TRUE)\n",
    "PicQ=cut(heart$PicInd,breaks=quantile(heart$PicInd,c(0,.33,.66,1)),labels=c(\"PicA\",\"PicB\",\"PIcC\"),include.lowest = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppression d'une modalité trop rare."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "heart[heart$Cardio==\"cdB\",\"Cardio\"]=\"cdA\"\n",
    "heart$Cardio=factor(heart$Cardio,exclude=NULL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Création finale du *dataFrame*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "heartT=data.frame(heart,AgeQ,TensQ,CholQ,FreQ,PicQ)\n",
    "summary(heartT)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Analyse exploratoire des variables rendues qualitatives\n",
    "L'étude exploratoire montre l'importance des variables qualitatives. en particulier, l'âge n'intervient pas comme risque de la même façon en fonction du genre des patients. Aussi, une nouvelle variable croisant sexe et âge est construite et une version simplifiée, avec 3 modalités au lieu de 6, est ajoutée à la base."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sélection des seules variables qualitatives\n",
    "heartQ=heartT[,-c(1,4,5,8,10)]\n",
    "summary(heartQ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Nouvelle variable pour intégrer l'interaction Sexe x Age\n",
    "SexAge=heartQ$Sexe:heartQ$AgeQ\n",
    "heartQ=data.frame(heartQ,\"SexAge\"=SexAge)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simplifier la variable en regroupant les jeunes, hommes et femmes, dans la même classe\n",
    "SxAg=as.factor(rep(c(\"sHFageA\",\"sFageBC\",\"sHageBC\"),90))\n",
    "SxAg[heartQ$AgeQ==\"AgeA\"]=\"sHFageA\"\n",
    "SxAg[heartQ$SexAge==\"sxF:AgeB\" | heartQ$SexAge==\"sxF:AgeC\"]=\"sFageBC\"\n",
    "SxAg[heartQ$SexAge==\"sxM:AgeB\" | heartQ$SexAge==\"sxM:AgeC\"]=\"sHageBC\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Complétion du data frame\n",
    "heartQ2=data.frame(heartQ[,-15],SxAg)\n",
    "summary(heartQ2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ajout de variables ou caractéristiques (*features*)\n",
    "Cette étude teste la stratégie suivante: ajouter des nouvelles variables ou *features* sous la forme des composantes principales issues de l'AFCM du tableau disjoinctif complet. Si toutes les variables étaient quantitatives, cela reviendrait à ajouter les composantes principales d'une analyse en composantes principales mais en prenant toutes celles qualitatives plus les quantitatives découpées en classes, c'est une AFCM qui convient.\n",
    "\n",
    "Cette approche a un autre avantage, elle permet de considérer des méthodes de classification supervisée qui n'autorisent que des variables quantitatives comme variables prédictives tout en conservant l'information apportée par les variables qualitatives. \n",
    "\n",
    "**Q** Quelle autre approche peut être mise en oeuvre pour contourner ce problème?\n",
    "\n",
    "La pratique qui consiste à ajouter de nouvelles caractéristiques des données est habituelles en apprentissage machine, notamment dans les concours de type [kaggle](https://www.kaggle.com/) pour tenter d'améliorer la qualité de prévision."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Construction des composantes principales\n",
    "Les nouvelles variables (composantes principales) sont construites à l'aide d'une AFCM disponible dans la librairie *FactoMineR*. Elles sont encore les composantes principales de l'ACP, avec la métrique dite du Chi2, du tableau disjonctif complet.\n",
    "\n",
    "Une première exécution fournit la représentation graphique de toutes les modalités dans le premier plan factoriel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(FactoMineR)\n",
    "afc=MCA(heartQ2,quali.sup=c(1,9,10),graph=F)\n",
    "options(repr.plot.width=5, repr.plot.height=5)\n",
    "par(mar=c(2.1, 2.1, .1, .1))\n",
    "plot.MCA(afc,invisible=c(\"ind\"),cex=0.8,habillage=\"quali\",palette=palette(rainbow(14)),title=\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "L'interprétation de ce graphique était un des objectifs de la [phase exploratoire](https://github.com/wikistat/Exploration/blob/master/Diag-coro/Explo-R-DiagCoro.ipynb). Remarquer simplement la position des modalités liées à la présence ou non de la pathologie: `CoP, CoA`, et celles représentant les facteurs de risque, dont `sHageBC` (hommes plus âgés), alors que le cholestérol (`CholA, B, C`) augmente avec l'âge (tous genres confondus) sans être nécessairement lié avec l'aggravation du risque."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Regrouper toutes les variables dans un même dataFrame\n",
    "heartT=data.frame(heart,heartQ2[,c(10:15)])\n",
    "summary(heartT)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2  Construction des échantillons\n",
    "Un objectif est ici de déterminer si l'ajout de nouvelles \"*features*\": les composantes principales de l'afcm, ont un intérêt ou non pour améliorer la prévision.\n",
    "\n",
    "#### Séparation des échantillons apprentissage et test\n",
    "**Q** Qu'est ce qui nécessite cette démarche?\n",
    "\n",
    "Choisir ci-dessous une initialisation différente du générateur de nombre aléatoire (XX)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "XX=11"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# séparation apprentissage / test\n",
    "set.seed(XX)\n",
    "npop=nrow(heartT)\n",
    "\n",
    "ntest=100\n",
    "testi=sample(1:npop,ntest)\n",
    "appri=setdiff(1:npop,testi)\n",
    "\n",
    "datapp=heartT[appri,]\n",
    "datest=heartT[testi,]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ajout des **features**\n",
    "Les étapes suivantes calculent les composantes principales de l'AFCM ou coordonnées des individus dans la base des \n",
    "vecteurs propres. \n",
    "\n",
    "**Q** Les individus, patients, de l'étude appartenant à l'échantillon test sont considérés comme *individus supplémentaires*. Pourquoi?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(FactoMineR)\n",
    "summary(heartQ)  \n",
    "afc=MCA(heartQ2[,-9],ind.sup=testi,ncp=9,graph=F)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "summary(afc$ind$coord)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "summary(afc$ind.sup$coord)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Construction des bases de données de l'échantillon d'apprentissage et de test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Ajout des nouvelles variables\n",
    "datapp=data.frame(datapp,afc$ind$coord)\n",
    "datest=data.frame(datest,afc$ind.sup$coord)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "summary(datapp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 Modélisation par régression logistique linéaire\n",
    "Le problème à résoudre est la prévision de la présence d'une coronopathie. C'est donc un problème de discrimination binaire que beaucoup de méthodes peuvent aborder. L'objectif est explicatif pour mieux comprendre la maladie et identifier les facteurs de risque mais, dans l'intérêt du patient, il peut n'être que prédictif par l'utilisation d'un modèle ou algorithme \"boîte noire\" si celui-ci fournit de meilleures prévisions.\n",
    "\n",
    "Les premiers modèles sont linéaires, sans interactions.\n",
    "\n",
    "La sélection de variables et donc la complexité du modèle sont optimisées avant de calculer l'erreur de prévision sur l'échantillon test et ce pour différents modèles. Retenir le \"meilleur\".\n",
    "\n",
    "**Q** Quels sont les principaux critères et algorithmes de sélection de variables utilisables pour cet objectif? \n",
    "\n",
    "**Q** Quels sont ceux mis en oeuvre ci-dessous?\n",
    "\n",
    "#### Sur les seules variables initiales quantitatives et qualitatives"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res0.log=glm(Coro~1,data=datapp,family=\"binomial\")\n",
    "# Changer le paramètre \"trace\" pour avoir la liste des modèles testés\n",
    "res0.log.step=step(res0.log,scope=list(lower=~1,upper=~(Age+Sexe+Douleur+Tension+Cholest+\n",
    "                  AngInd+PicInd+PenteInd+Nvais+Thal+Sucre+Cardio+FreqM+SxAg),\n",
    "                  direction=\"both\"),trace=0)\n",
    "# Résultat des tests\n",
    "anova(res0.log.step, test=\"Chisq\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Commenter la liste des variables sélectionées, les p-valeurs des tests."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(res0.log.step)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Interpréter les facteurs de risque. Est-ce compatible avec la représentation de l'AFCM?\n",
    "\n",
    "La qualité de prévision du modèle est estimée par la prévision de l'échantillon test et mise en évidence par le calcul de la *matrice de confusion*. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prévision de l'échantillon test\n",
    "table(predict(res0.log.step,datest,type=\"response\")>0.5,datest[,14]==\"CoP\") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Que dire de la qualité de ce modèle? De la forme de la matrice? Du choix du seuil égal à 0.5?\n",
    "\n",
    "#### Sur l'ensemble des variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res.log=glm(Coro~1,data=datapp,family=\"binomial\")\n",
    "res.log.step=step(res.log,scope=list(lower=~1,upper=~(Age+Sexe+Douleur+Tension+Cholest+\n",
    "                AngInd+PicInd+PenteInd+Nvais+Thal+Sucre+Cardio+FreqM+SxAg+Dim.1+Dim.2+Dim.3+\n",
    "                                    Dim.4+Dim.5+Dim.6+Dim.8+Dim.9)),direction=\"both\",trace=0)\n",
    "anova(res.log.step, test=\"Chisq\")\n",
    "table(predict(res.log.step,datest,type=\"response\")>0.5,datest[,14]==\"CoP\") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** L'ajout des composantes principales semble-t-il pertinent?\n",
    "#### Sur les seules composantes de l'AFCM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res.log=glm(Coro~1,data=datapp,family=\"binomial\")\n",
    "res.log.step=step(res.log,scope=list(lower=~1,upper=~(Dim.1+Dim.2+Dim.3+Dim.4+\n",
    "                            Dim.5+Dim.6+Dim.8+Dim.9)),direction=\"both\",trace=0)\n",
    "anova(res.log.step, test=\"Chisq\")\n",
    "table(predict(res.log.step,datest,type=\"response\")>0.5,datest[,14]==\"CoP\") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Que pensez de la capacité des composantes pricipales de l'AFCM à résumer l'information pertinente apportée par les variables? Justifier.\n",
    "\n",
    "**Q** Que perd on dans cette démarche?\n",
    "\n",
    "Identifier le meilleur modèle en terme de prévision."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.4  Modèlisation par régression logistique quadratique\n",
    "Les modèles précédents en prennent pas en compte de possibles interactions entre les variables à l'exception de celle `age x sexe` introduite à la main. Les commandes suivantes recherchent directement un modèle avec interactions au plus d'ordre 2.\n",
    "\n",
    "Caractériser la stratégie mise en oeuvre."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res.log=glm(Coro~1,data=datapp,family=\"binomial\")\n",
    "# Supprimer les warnings pour simplifier la sortie\n",
    "options(warn=-1)\n",
    "res.log.step=step(res.log,scope=list(lower=~1,upper=~(Age+Sexe+Douleur+Tension+Cholest+\n",
    "                  AngInd+PicInd+PenteInd+Nvais+Thal+Sucre+Cardio+FreqM+SxAg)**2,\n",
    "                  direction=\"both\"),trace=0)\n",
    "# Résultat des tests\n",
    "anova(res.log.step, test=\"Chisq\")\n",
    "# remettre les warnings\n",
    "options(warn=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "table(predict(res.log.step,datest,type=\"response\")>0.5,datest[,14]==\"CoP\") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Commentaire sur les interactions sélectionnées?\n",
    "\n",
    "**Q** Le modèle quadratique est-il  meilleur?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.5 Arbre de décision \n",
    "A titre illustratif, un autre modèle est testé et qui prend en compte, par construction, les iteractions. Il s'agit d'un [arbre binaire de discrimination](http://wikistat.fr/pdf/st-m-app-cart.pdf).\n",
    "\n",
    "Ce modèle est estimé sur les données initiales"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(rpart)\n",
    "library(partykit)\n",
    "res.tree=rpart(Coro~.,data=datapp[,1:20],cp=0)\n",
    "par(mar=c(2.1, 2.1, .1, .1))\n",
    "plot(as.party(res.tree),gp = gpar(fontsize = 11))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Que dire de la facilité d'interprétation, des facteurs de risque identifiés par ce modèle? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "table(predict(res.tree,datest,type=\"class\"),datest[,14])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Que dire de la qualité de la prévision? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.6 Comparaison des courbes ROC\n",
    "Retenir le meilleur modèle de régression logistique et l'introduire dans le code ci-dessous. Attention ce résultat dépend de la séparation initiale de l'échantillon entre apprentissage et test; il ne correspond pas nécessairement au choix réalisé ci-dessous."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "library(ROCR)\n",
    "roclogit=predict(res0.log.step,newdata=datest,type=\"response\")\n",
    "predlogit=prediction(roclogit,datest[,14])\n",
    "perflogit=performance(predlogit,\"tpr\",\"fpr\")\n",
    "\n",
    "roctree=predict(res.tree,newdata=datest,type=\"prob\")[,2]\n",
    "predtree=prediction(roctree,datest[,14])\n",
    "perftree=performance(predtree,\"tpr\",\"fpr\")\n",
    "\n",
    "\n",
    "options(repr.plot.width=4, repr.plot.height=3)\n",
    "par(mar=c(2.1, 2.1, .1, .1))\n",
    "palette(\"default\")\n",
    "plot(perflogit,col=1,lwd=2,print.cutoffs.at=c(0.1, 0.5))\n",
    "plot(perftree,col=2,lwd=2,add=T,print.cutoffs.at=c(0.1, 0.5, 0.98))\n",
    "\n",
    "legend(\"bottomright\",legend=c(\"logit\",\"rpart\"),\n",
    "       y.intersp=2,col=c(1,2),lty=1,lwd=2,text.width=0.15,ncol=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Que dire de façon générale sur le choix de la meilleure stratégie pour déterminer le paquet de variables à utiliser? Distinguer selon l'objectif de prévision ou d'interprétation. \n",
    "\n",
    "**Q** Que signifient les valeurs présentes sur les courbes ROC? \n",
    "\n",
    "**Q** Quelle conclusion tirer de ce graphique en fonction du taux de faux positif acceptable?\n",
    "\n",
    "**Remarque** importante: l'échantillon test est de **petite taille** et donc l'estimation de l'erreur sujette à caution (faible précision) car affectée d'une variance importante. Une estimaiton plus précise de la distribution de cette erreur (validation croisée *Monte Carlo*) ainsi qu'une comparaison systématique avec les autres méthodes et algorithmes d'apprentissage (*k*-nn, neurones, *boosting, random forest*, SVM) sont proposées dans un autre [tutoriel](https://github.com/wikistat/Apprentissage/blob/master/Diag-coro/Apprent-R-DiagCoro.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(randomForest)\n",
    "res.tree=randomForest(Coro~.,data=datapp[,1:20])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "table(predict(res.tree,datest,type=\"class\"),datest[,14])"
   ]
  }
 ],
 "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
}