TutosRudim/Cal5-R-Classif.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",
    "</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [Scénarios d'Exploration Statistique](https://github.com/wikistat/Exploration)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Classification non supervisée 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",
    "# Représentations des classes dans le [MDS](http://wikistat.fr/pdf/st-m-explo-mds.pdf), l'[ACP](http://wikistat.fr/pdf/st-m-explo-acp.pdf) ou l'[AFCM](http://wikistat.fr/pdf/st-m-explo-afcm.pdf)\n",
    "#### Résumé \n",
    "[Classification non supervisée](http://wikistat.fr/pdf/st-m-explo-class.pdf) (CAH, *kmeans*, *pam*) et représentation des classes dans un plan factoriel avec des données sous différentes formes: tableau de distance et [MDS](http://wikistat.fr/pdf/st-m-explo-mds.pdf), variables quantitatives et [ACP](http://wikistat.fr/pdf/st-m-explo-acp.pdf), variables qualitatives et [AFCM](http://wikistat.fr/pdf/st-m-explo-afcm.pdf). Application à des exemples élémentaires."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1 Introduction\n",
    "Le nombre de méthodes de classification non supervisée (*clustering*) et le nombre d'options dans chacune d'elles conduisent à une combinatoire de possibilités assez impressionnante. Le but est de mettre en évidence, de façon heuristique et dans des cas simples, quels sont les paramètres qui influencent, de visu, les classes obtenues. Un rôle important est donc attribué à la représentation des classes dans un plan factoriel. Le choix de la méthode: [analyse en composantes principales](http://wikistat.fr/pdf/st-m-explo-acp.pdf), [analyse factorielle multiple des correspondances](http://wikistat.fr/pdf/st-m-explo-afcm.pdf)  ou encore [positionnement multidimensionnel](http://wikistat.fr/pdf/st-m-explo-mds.pdf), dépend du type (quantitatives, qualitatives, distances) des données analysées.\n",
    "\n",
    "La plupart des analyses ci-dessous peuvent être reproduites avec SAS mais de façon nettement moins souple et pas toutes (PAM n'est pas implémenté). \n",
    "\n",
    "## 2 Positionnement multidimensionnelle ou *MDS*\n",
    "### 2.1 Objectif\n",
    "le *MultiDimensional Scaling, MDS* ou encore *ACP d'un tableau de distances* s'applique à des données archivées sous la forme d'une matrice ($n\\times n$) de distances ou dissimilarités. Cet algorithme est appliqué sur un exemple simple comparant des distances entre villes. \n",
    "\n",
    "Les données se présentent donc sous la forme du triangle inférieur d'une matrice symétrique, par construction, et contenant les distances kilométriques routières, donc non-euclidiennes, de 47 villes française ou proches prises 2 à 2 (Source : carte IGN). \n",
    "\n",
    "### 2.2 Lecture des données sous la forme d'une matrice triangulaire\n",
    "La lecture d'une matrice triangulaire inférieure en tant que matrice de distance pose quelques difficultés dans R dans la gestion du type des objets. Ceci nécessite une séquence de commandes pour la bonne mise en forme."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mdsville=read.table(\"Data/mdsville.dat\",fill=TRUE)\n",
    "# extraction des noms des villes\n",
    "villes=as.character(mdsville[2:48,1])\n",
    "# extraction des valeurs des distances\n",
    "m=mdsville[2:48,2:48]\n",
    "# transformation en une matrice alphanumérique\n",
    "m=as.matrix(m)\n",
    "# retour au numérique avec données manquantes\n",
    "m=as.numeric(m)\n",
    "# reformattage en une matrice\n",
    "m=matrix(m,47,47)\n",
    "# adjonction des noms des villes en ligne et colonne\n",
    "dimnames(m)[[1]]=villes\n",
    "dimnames(m)[[2]]=villes\n",
    "# transformation en un objet de type distance\n",
    "d=as.dist(m,diag=TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 MDS des villes\n",
    "En revanche, le MDS est immédiat. On remarque que le choix de la dimension, ici $k=2$ est à traiter comme en ACP. Une visite du manuel permettrait de voir comment calculer toutes les valeurs singulières afin d'aider ce choix par leur représentation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mds = cmdscale(d, k=2)\n",
    "plot(mds, type=\"n\", xlab=\"cp1\", ylab=\"cp2\")\n",
    "text(mds,villes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Commenter la représentation obtenue par rapport à la carte géographique.\n",
    "\n",
    "## 3 Classification à partir d'un tableau de distances\n",
    "\n",
    "Les mêmes données sont reprises en vue de regrouper les villes en classes homogènes au regard de leurs distances respectives.\n",
    "\n",
    "### 3.1 CAH et représentation par MDS\n",
    "**Q** Quel est le graphe ci-dessous? Commenter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "chv = hclust(d, method=\"ward.D\")\n",
    "options(repr.plot.width=7, repr.plot.height=5)\n",
    "plot(chv,main=NULL,sub=\"\",xlab=\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Comment est définie la corrélation cophénétique? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cophenetic correlation \n",
    "cor(cophenetic(chv),d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Quel est le graphe ci-dessous? Quelle décision en tirer?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=4, repr.plot.height=4)\n",
    "plot(chv$height[46:30],xlab=\"nb de classes\",ylab=\"Hauteur\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Quel est le graphe ci-dessous? Commentaire."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "color=cutree(chv,k=5)\n",
    "library(cluster)\n",
    "options(repr.plot.width=6, repr.plot.height=5)\n",
    "plot(silhouette(color,d))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Il est possible de calculer puis représenter la silhouette moyenne pour plusieurs valeurs du nombre de classes et tenter de minimiser cette valeur; cette stratégie conduisant à un choix ``trivial'' $(k=3)$ est laissée de côté.\n",
    "\n",
    "Les données provenant d'une matrice de distances, le MDS s'impose pour une représentation fatorielle des villes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=5, repr.plot.height=5)\n",
    "plot(mds, type=\"n\", xlab=\"cp1\", ylab=\"cp2\")\n",
    "# représentation avec des couleurs\n",
    "text(mds,villes,col=color)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Influence déterminante du choix de la distance entre groupe.\n",
    "\n",
    "**Q** Quelles sont les options possibles de l'option `method` de la fonction `hclust`?\n",
    "\n",
    "**Q** Que signifie le choix `single` ci-dessous?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "chv <- hclust(d, method=\"single\")\n",
    "cor(cophenetic(chv),d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Que devient la corrélation cophénétique?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=7, repr.plot.height=5)\n",
    "plot(chv,main=NULL,sub=\"\",xlab=\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=5, repr.plot.height=4)\n",
    "plot(chv$height[46:30],xlab=\"nb de classes\",ylab=\"Hauteur\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Le choix du nombre de classes est-il facile?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "color=cutree(chv,k=5)\n",
    "options(repr.plot.width=5, repr.plot.height=5)\n",
    "plot(mds, type=\"n\", xlab=\"cp1\", ylab=\"cp2\")\n",
    "text(mds,villes,col=color) # représentation avec des couleurs "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tester les autres options de distances entre groupes. \n",
    "\n",
    "**Q** Conclusion: le critère de plus grande corrélation cophénétique est-il à considérer pour choisir l'option de distance entre les groupes?\n",
    "\n",
    "### 3.2 PAM et représentation par MDS\n",
    "L'algorithme de réallocation k-means n'est pas adapté à une matrice de distances ou de dissimilarités; en revanche, PAM (partition around medoïds) est opérationnel si le nombre d'observations n'est pas trop important; sinon utiliser l'adaptation {\\tt clara}. Ces fonctions sont disponibles dans la librairie `cluster`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(cluster)\n",
    "# faire varier le nombre de classes\n",
    "plot(silhouette(pam(d,6)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pamv=pam(d,6)\n",
    "color=pamv$clustering\n",
    "plot(mds, type=\"n\", xlab=\"cp1\", ylab=\"cp2\")\n",
    "# représentation avec des couleurs\n",
    "text(mds,villes,col=color)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# le même avec des ellipses\n",
    "clusplot(d,pamv$clustering,diss=TRUE,labels=2,color=TRUE,col.txt=pamv$clustering,main=\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Remarque*: le critère de `clusGap` pour la sélection du nombre de classe ne s'applique pas à une matrice de distances mais il pourrait s'appliquer aux coordonnées principales issues du MDS."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4 Classification de données quantitatives\n",
    "### 4.1 Données OCDE\n",
    "Les données du fichier `ocdeR.dat` sont celles étudiées comme exemple d'[ACP cubique](http://wikistat.fr/pdf/st-m-explo-acp.pdf) avec SAS. Comme elles sont connues par une représentation euclidienne, la distance considérée par défaut est euclidienne mais d'autres choix sont possibles (manhattan ou valeurs absolues...). \n",
    "\n",
    "\n",
    "### 4.2 CAH et représentation par ACP\n",
    "Lecture des données puis construction de la matrice de distance.\n",
    "\n",
    "**Attention** Il est important de centrer réduire des variables quantitatives avant de faire une classification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ocde=read.table(\"Data/ocdeR.dat\")\n",
    "ds=dist(scale(ocde))\n",
    "# classification hiérarchique\n",
    "hc.ds = hclust(ds,method=\"ward.D\")\n",
    "options(repr.plot.width=8, repr.plot.height=5)\n",
    "plot(hc.ds,main=\"\")  # dendogramme"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# choix du nombre de classes\n",
    "options(repr.plot.width=4, repr.plot.height=4)\n",
    "plot(hc.ds$height[67 :58],type=\"b\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "color= cutree(hc.ds,k=4) # couleurs des classes\n",
    "# silhouette\n",
    "library(cluster)\n",
    "plot(silhouette(color,ds))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# corrélation cophénétique\n",
    "cor(cophenetic(hc.ds),ds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# acp pour représentations\n",
    "library(FactoMineR)\n",
    "acp=PCA(ocde,ncp=13,graph=F) \n",
    "# graphe de l'acp\n",
    "options(repr.plot.width=5, repr.plot.height=5)\n",
    "plot(acp,choix=\"ind\", habillage=\"ind\",col.hab=rep(1:17, c(rep(4,17))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# graphe de l'acp avec les couleurs des classes de la cah\n",
    "plot(acp,choix=\"ind\", habillage=\"ind\",col.hab=color)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Refaire tourner l'algorithme en retirant la réduction des variables (fonction scale), même chose en remplaçant `method=\"ward.D\"` par `method=\"single\"`.\n",
    "\n",
    "**Q** Noter l'importance de la normalisation des variables.\n",
    "\n",
    "**Q** Que choisir comme méthode?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.3 Algorithme de ré-allocation et représentation par ACP\n",
    "Les données étant quantitatives, k-means est opérationnel. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# version avec réduction\n",
    "#  choix de 4 classes suggérées par la cah\n",
    "kocde=kmeans(scale(ocde),4) \n",
    "color=kocde$cluster\n",
    "plot(acp,choix=\"ind\", habillage=\"ind\",col.hab=color)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# version sans réduction\n",
    "kocde=kmeans(ocde,4) \n",
    "color=kocde$cluster\n",
    "plot(acp,choix=\"ind\", habillage=\"ind\",col.hab=color)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# utilisation de clusGap\n",
    "library(cluster)\n",
    "clusGap(ocde, FUN = kmeans, K.max = 8, B = 60)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plusieurs critères sont concurrents pour sélectionner un $k$ optimal une fois les calculs de `clusGap` réalisés par *bootsrap*. le choix reste complexe.\n",
    "\n",
    "### 4.4 Classification sur composantes principales\n",
    "Dans des situaitons de très grande dimension, il peut être intéressant de réduire celle-ci par ACP avant d'opérer une classification. \n",
    "\n",
    "La librairie `FactoMineR` propose cette possibilité avec la production de graphiques élaborées mais sans doute pas tous utiles et qui bloquent jupyter. Un autre paramètre est susceptible d'avoir un impact, le nombre `ncp` de composantes retenues. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "acp=PCA(ocde,ncp=5,graph=F)\n",
    "res.hcpc=HCPC(acp,graph=F,nb.clust=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(acp, choix=\"ind\", habillage=\"ind\",col.hab=res.hcpc$data.clust$clust)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Que dire du choix de 4 classes? Revenir à 3.\n",
    "\n",
    "**Q** Quelle option de classification vous semble la plus pertinente sur ces données OCDE?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5 Classification de données qualitatives\n",
    "Plusieurs stratégies sont proposées pour classées des données qualitatives ou mélange de variables quantitatives et qualitatives. La plus simple consiste à recoder les variables quantitatives en qualitatives puis de calculer des scores à l'aide des composantes d'une [AFCM](http://wikistat.fr/pdf/st-m-explo-afcm.pdf).\n",
    "### 5.1 Données décrivant les races de chien\n",
    "\n",
    "Les données (race des chiens) sont celles illustrant l'utilisation de l'AFCM; R est nettement plus souple que SAS pour enchaîner de telles analyses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "chiens=read.csv(\"chiens.Data/csv\",row.names=1)\n",
    "summary(chiens)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.2 CAH sur composantes de l'AFCM\n",
    "Nous reprenons l'AFCM sur tableau disjonctif complet fournissant des scores ou variables quantitatives sur les individus, ici les races de chiens, il est alors facile de les classer. L'enchaînement par une CAH est la combinaison la plus simple à réaliser à l'aide de la librairie `FactoMineR`. Attention néanmoins de bien contrôler les options qui toutes sont prises par défaut. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(FactoMineR)\n",
    "#afcm avec la fonction en supplémentaire\n",
    "afcm=MCA(chiens,quali.sup=7,graph=F)\n",
    "plot(afcm, choix=\"ind\",habillage=\"quali\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`FactoMineR`propose des aides à l'interprétation des classees. La variable classe est croisée \n",
    "avec chacune des variables. Un test du chi2 apprécie la liaison."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res.hcpc=HCPC(afcm, graph=FALSE,order=FALSE,nb.clust=4)\n",
    "res.hcpc$desc.var$test.chi2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Avant de mettre en évidence les modalités les plus \"présentes\" dans chaque classe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res.hcpc$desc.var$category"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Justification a posteriori du choix du nombre de classes. Le choix par défaut (10) n'a pas d'intérêt."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=5, repr.plot.height=5)\n",
    "plot(res.hcpc,choice=\"tree\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(afcm, choix=\"ind\", invisible=\"var\",col.ind=res.hcpc$data.clust$clust)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "table(res.hcpc$data.clust$clust,res.hcpc$data.clust$fonction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "D'autres outils (arbre de discrimination) s'avéreront plus explicites pour interpréter les classes.\n",
    "\n",
    "**Q** Combien de dimensions sont retenues par défaut dans cette AFCM? \n",
    "\n",
    "**Q** Quelle distance est prise en compte entre les individus?\n",
    "\n",
    "**Q** Quelle est l'option par défaut pour la distance entre les groupes?\n",
    "\n",
    "D'autre par une option positionnée (`consol=TRUE`) fait enchaîner automatiquement un algorithme *k-means*.\n",
    "\n",
    "**Q** Que dire des options suivantes ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "afcm=MCA(chiens,quali.sup=7, ncp=10, graph=F)\n",
    "res.hcpc=HCPC(afcm,method=\"single\",graph=F,nb.clust=4)\n",
    "table(res.hcpc$data.clust$clust,res.hcpc$data.clust$fonction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "En fait les classes se distinguent assez facilement et restent stables sur ces données sauf en cas d'utilisation du saut minimum.\n",
    "\n",
    "Bien entendu, la fonction `HCPC` de `FactoMineR` pourrait être remplacée par tout autre algorithme de classification avec ses propres options: `hclust, kmeans, pam...`."
   ]
  }
 ],
 "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
}