Analyse multivariée: association, partitionnement et ordination

Pour cette section, vous aurez besoin que les modules suivants soient installés.


In [1]:
import numpy as np
import pandas as pd
import matplotlib.pylab as plt
%matplotlib inline

Espaces d'analyse

Abondance et occurence

L'abondance est le décompte d'espèces observées, tandis que l'occurence est la présence ou l'absence d'une espèce. Le tableau suivant contient des données d'abondance.


In [2]:
abundance = pd.DataFrame({'Bruant familier': [1, 0, 0, 3],
                          'Citelle à poitrine rousse': [1, 0, 0, 0],
                          'Colibri à gorge rubis': [0, 1, 0, 0],
                          'Geai bleu': [3, 2, 0, 0],
                          'Bruant chanteur': [1, 0, 5, 2],
                          'Chardonneret': [0, 9, 6, 0],
                          'Bruant à gorge blanche': [1, 0, 0, 0],
                          'Mésange à tête noire': [20, 1, 1, 0],
                          'Jaseur boréal': [66, 0, 0, 0]},
                         index = ['Site 1', 'Site 2', 'Site 3', 'Site 4'])
abundance


Out[2]:
Bruant chanteur Bruant familier Bruant à gorge blanche Chardonneret Citelle à poitrine rousse Colibri à gorge rubis Geai bleu Jaseur boréal Mésange à tête noire
Site 1 1 1 1 0 1 0 3 66 20
Site 2 0 0 0 9 0 1 2 0 1
Site 3 5 0 0 6 0 0 0 0 1
Site 4 2 3 0 0 0 0 0 0 0

Ce tableau peut être rapidement transformé en données d'occurence, qui ne comprennent que l'information booléenne de présence (noté 1) et d'absence (noté 0).


In [3]:
occurence = (abundance > 0).astype(int)
occurence


Out[3]:
Bruant chanteur Bruant familier Bruant à gorge blanche Chardonneret Citelle à poitrine rousse Colibri à gorge rubis Geai bleu Jaseur boréal Mésange à tête noire
Site 1 1 1 1 0 1 0 1 1 1
Site 2 0 0 0 1 0 1 1 0 1
Site 3 1 0 0 1 0 0 0 0 1
Site 4 1 1 0 0 0 0 0 0 0

L'espace des espèces (ou des variables ou descripteurs) est celui où les espèces forment les axes et où les sites sont positionnés dans cet espace. Il s'agit d'une perspective en mode R, qui permet principalement d'identifier quels espèces se retrouvent plus courrament ensemble


In [4]:
%pylab inline
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

species = ['Bruant familier', 'Bruant chanteur', 'Chardonneret']

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(abundance[species[0]],
           abundance[species[1]],
           abundance[species[2]], marker='o')
ax.set_xlabel(species[0])
ax.set_ylabel(species[1])
ax.set_zlabel(species[2])


Populating the interactive namespace from numpy and matplotlib
/home/essicolo/bin/anaconda3/lib/python3.6/site-packages/IPython/core/magics/pylab.py:160: UserWarning: pylab import has clobbered these variables: ['plt']
`%matplotlib` prevents importing * from pylab and numpy
  "\n`%matplotlib` prevents importing * from pylab and numpy"
Out[4]:
Text(0.5,0,'Chardonneret')

Dans l'espace des sites (ou les échantillons ou objets), on transpose la matrice d'abondance. On passe ici en mode Q, et l'on peut observer quels échantillons sont similaires.


In [5]:
sites = abundance.index[0:3]
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(abundance.T[sites[0]],
           abundance.T[sites[1]],
           abundance.T[sites[2]], marker='o')
ax.set_xlabel(sites[0])
ax.set_ylabel(sites[1])
ax.set_zlabel(sites[2])


Out[5]:
Text(0.5,0,'Site 3')

Environnement

L'espace de l'environnement comprend souvent un autre tableau contenant l'information sur l'environnement où se trouve les espèces: les coordonnées et l'élévation, la pente, le pH du sol, la pluviométrie, etc.

Analyse d'association

Nous utiliserons le terme association come une mesure pour quantifier la ressemblance ou la différence entre deux objets (échantillons) ou variables (descripteurs).

Alors que la corrélation et la covariance sont des mesures d'association entre des variables (analyse en mode R), la similarité et la distance sont deux types de une mesure d'association entre des objets (analyse en mode Q). Une distance de 0 est mesuré chez deux objets identiques. La distance augmente au fur et à mesure que les objets sont dissociés. Une similarité ayant une valeur de 0 indique aucune association, tandis qu'une valeur de 1 indique une association parfaite. À l'opposé, la dissimilarité est égale à 1-similarité.

La distance peut être liée à la similarité par la relation:

$$distance=\sqrt{1-similarité}$$

ou

$$distance=\sqrt{dissimilarité}$$

La racine carrée permet, pour certains indices de similarité, d'obtenir des propriétés euclédiennes. Pour plus de détails, voyez le tableau 7.2 de Legendre et Legendre (2012).

Les matrices d'association sont généralement présentées comme des matrices carrées, dont les dimensions sont égales au nombre d'objets (mode Q) ou de vrariables (mode R) dans le tableau. Chaque élément ("cellule") de la matrice est un indice d'association entre un objet (ou une variable) et un autre. Ainsi, la diagonale de la matrice est un vecteur nul (distance ou dissimilarité) ou unitaire (similarité), car elle correspond à l'association entre un objet et lui-même.

Puisque l'association entre A et B est la même qu'entre B et A, et puisque la diagonale retourne une valeur convenue, il est possible d'exprimer une matrice d'association en mode "compact", sous forme de vecteur. Le vecteur d'association entre des objets A, B et C contiendra toute l'information nécessaire en un vecteur de trois chiffres, [AB, AC, BC], plutôt qu'une matrice de dimension $3 \times 3$. L'impact sur la mémoire vive peut être considérable pour les calculs comprenant de nombreuses dimensions.

En Python, les calculs de similarité et de distances peuvent être effectués avec le module scipy. La fonction pdist permet de calculer les indices d'association en forme compact, que l'on peut étaler en forme carrée avec la fonction squareform.

Nous verons plus tard les méthodes de mesure de similarité et de distance plus loin. Pour l'instant, utilisons la méthode de Jaccard pour une démonstration sur des données d'occurence.


In [6]:
from scipy.spatial.distance import pdist, squareform
from sklearn import datasets

squareform(pdist(occurence, 'jaccard'))


Out[6]:
array([[ 0.        ,  0.77777778,  0.75      ,  0.71428571],
       [ 0.77777778,  0.        ,  0.6       ,  1.        ],
       [ 0.75      ,  0.6       ,  0.        ,  0.75      ],
       [ 0.71428571,  1.        ,  0.75      ,  0.        ]])

Remarquez que scipy retourne une matrice dont la diagonale est de 0. La diagonale est l'association d'un objet avec lui-même. Or la similarité d'un objet avec lui-même devrait être de 1! En fait, par convention scipy retourne des dissimilarités, non pas des similarités. La matrice de distance serait donc calculée en extrayant la racine carrée des éléments de la matrice de dissimilarité:


In [7]:
dissimilarity = squareform(pdist(occurence, 'jaccard'))
distance = np.sqrt(dissimilarity)
distance


Out[7]:
array([[ 0.        ,  0.8819171 ,  0.8660254 ,  0.84515425],
       [ 0.8819171 ,  0.        ,  0.77459667,  1.        ],
       [ 0.8660254 ,  0.77459667,  0.        ,  0.8660254 ],
       [ 0.84515425,  1.        ,  0.8660254 ,  0.        ]])

Dans le chapitre sur l'analyse compositionnelle, nous avons abordé les significations différentes que peuvent prendre le zéro. L'information fournie par un zéro peut être différente selon les circonstances. Dans le cas d'une variable continue, un zéro signifie généralement une mesure sous le seuil de détection. Deux tissus dont la concentration en cuivre est nulle ont une afinité sous la perspective de la concentration en cuivre. Dans le cas de mesures d'abondance (décompte) ou d'occurence (présence-absence), on pourra décrire comme similaires deux niches écologiques où l'on retrouve une espèce en particulier. Mais deux sites où l'on de retouve pas d'ours polaires ne correspondent pas nécessairement à des niches similaires! En effet, il peut exister de nombreuses raisons écologiques et méthodologiques pour lesquelles l'espèces ou les espèces n'ont pas été observées. C'est le problème des double-zéros (espèces non observées à deux sites), problème qui est amplifié avec les grilles comprenant des espèces rares.

La ressemblance entre des objets comprenant des données continues devrait être calculée grâce à des indicateurs symétriques. Inversement, les affinités entre les objets décrits par des données d'abondance ou d'occurence susceptibles de générer des problèmes de double-zéros devraient être évaluées grâce à des indicateurs asymétriques. Un défi supplémentaire arrive lorsque les données sont de type mixte.

Nous utiliserons la convention de scipy et nous calculerons la dissimilarité, non pas la similarité. Les mesures de dissimilarité sont calculées sur des données d'abondance ou des données d'occurence. Notons qu'il existe beaucoup de confusion dans la littérature sur la manière de nommer les dissimilarités (ce qui n'est pas le cas des distances, dont les noms sont reconnus). Dans les sections suivantes, nous noterons la dissimilarité avec un $d$ minuscule et la distance avec un $D$ majuscule.

Association entre objets (mode Q)

Objets: Abondance

La dissimilarité de Bray-Curtis est asymétrique. Elle est aussi appelée l'indice de Steinhaus, de Czekanowski ou de Sørensen. Il est important de s'assurer de bien s'entendre la méthode à laquelle on fait référence. L'équation enlève toute ambiguité. La dissimilarité de Bray-Curtis entre les points A et B est calculée comme suit.

$$d_{AB} = \frac {\sum \left| A_{i} - B_{i} \right| }{\sum \left(A_{i}+B_{i}\right)}$$

Utilisons pdist pour générer les vecteurs d'association, qui peuvent être transformés en matrices avec squareform. Le format "dictionnaire" de Python est pratique pour enregistrer la collection de vecteurs d'association que nous allons créer dans cette section.


In [8]:
associations_abund = {}
associations_abund['BrayCurtis'] = pdist(abundance, 'braycurtis')
squareform(associations_abund['BrayCurtis'])


Out[8]:
array([[ 0.        ,  0.94339623,  0.96190476,  0.95918367],
       [ 0.94339623,  0.        ,  0.44      ,  1.        ],
       [ 0.96190476,  0.44      ,  0.        ,  0.76470588],
       [ 0.95918367,  1.        ,  0.76470588,  0.        ]])

La dissimilarité de Bray-Curtis est souvent utilisée dans la littérature. Toutefois, la version originale de Bray-Curtis n'est pas tout à fait métrique (semimétrique). Conséquemment, la dissimilarité de Ruzicka (une variante de la dissimilarité de Jaccard pour les données d'abondance) est métrique, et devrait probablement être préféré à Bary-Curtis (Oksanen, 2006).

$$d_{AB, Ruzicka} = \frac { 2 \times d_{AB, Bray-Curtis} }{1 + d_{AB, Bray-Curtis}}$$

In [9]:
associations_abund['Ruzicka'] = associations_abund['BrayCurtis'] * 2 / (1+associations_abund['BrayCurtis'])

La dissimilarité de Kulczynski (aussi écrit Kulsinski) est asymétrique et semimétrique, tout comme celle de Bray-Curtis. Elle est calculée comme suit.

$$d_{AB} = 1-\frac{1}{2} \times \left[ \frac{\sum min(A_i, B_i)}{\sum A_i} + \frac{\sum min(A_i, B_i)}{\sum B_i} \right]$$

La fonction kulsinski de la version de scipy 0.19.1 ne donne pas les résultats qu'elle devrait (rapport de bogue). Une fonction maison peut être conçue.


In [10]:
def kulczynski(X):
    d_AB = []
    for i in range(X.shape[0]):
        sum_i = X.iloc[i, :].sum()
        for j in range(X.shape[0]):
            if (i>j):
                df_ = X.iloc[[i, j], :]
                min_ = df_.apply(min, axis=0).sum()
                sum_j = X.iloc[j, :].sum()
                d_AB.append(1-0.5*(min_/sum_i + min_/sum_j))
    return(d_AB)

In [11]:
associations_abund['Kulczynski'] = kulczynski(abundance)

Une approche commune pour mesurer l'association entre sites décrits par des données d'abondance est la distance de Hellinger. Notez qu'il s'agit ici d'une distance, non pas d'une dissimilarité. La première étape est de diviser chaque donnée d'abondance par l'abondance totale pour chaque site, puis d'extraire la racine carrée de chaque élément. Enfin, on calcule la distance euclidienne entre chaque site. Pour rappel, une distance euclidienne est la généralisation en plusieurs dimensions du théorème de Pythagore, $c = \sqrt{a^2 + b^2}$.

$$D_{AB} = \sqrt {\sum \left( \frac{A_i}{\sum A_i} - \frac{B_i}{\sum B_i} \right)^2}$$

In [12]:
Hel_abundance = abundance.apply(lambda x: np.sqrt(x / np.sum(x)), axis=1) # prétraitement
associations_abund['Hellinger'] = pdist(Hel_abundance, 'euclidean')

Toute comme la distance d'Hellinger, la distance de chord est calculée par une distance euclidienne sur des données d'abondance transformées de sorte que chaque ligne ait une longueur (norme) de 1.


In [13]:
from sklearn.preprocessing import normalize
chord_abundance = normalize(abundance, axis=1)
associations_abund['chord'] = pdist(chord_abundance, 'euclidean')

La métrique du chi-carré, ou $\chi$-carré, ou chi-square, donne davantage de poids aux espèces rares qu'aux espèces communes. Son utilisation est recommandée lorsque les espèces rares sont de bons indicateurs de conditions écologiques particulières (Legendre et Legendre, 2012, p. 308).

$$ d_{AB} = \sqrt{\sum _j \frac{1}{\sum y_j} \left( \frac{A_j}{\sum A} - \frac{B_j}{\sum B} \right)^2 } $$

La métrique peut être transformée en distance en la multipliant par la racine carrée de la somme totale des espèces dans la matric d'abondance ($X$).

$$ D_{AB} = \sqrt{\sum X} \times d_{AB} $$

scipy n'offre pas les association du chi-carré dans sa trousse d'outil. Cette fonction fera le travail.


In [14]:
def chi_square(X, method='distance'):
    X = np.array(X)
    row_sums = X.sum(axis=1, keepdims=True)
    column_sums = X.sum(axis=0)
    scaled_X = X / (row_sums * np.sqrt(column_sums))
    if method == 'distance':
        D_AB = pdist(scaled_X, metric='euclidean') * np.sqrt(X.sum())
    elif method=='metric':
        D_AB = pdist(scaled_X, metric='euclidean')
    else:
        print('Wrong method')
    return(D_AB)

In [15]:
associations_abund['Chi-square'] = chi_square(abundance)

Une mannière visuellement plus intéressante de présenter une matrice d'association est un graphique de type heatmap, que le module seaborn permet de générer rapidement. Il est possible de générer ces graphiques dans une boucle.


In [16]:
import seaborn as sns

plt.figure(figsize=(20, 10))
for i, key in enumerate(associations_abund):
    plt.subplot(2, 3, i+1)
    assoc_df = pd.DataFrame(squareform(associations_abund[key]), index=abundance.index, columns=abundance.index)
    sns.heatmap(assoc_df, annot=True, square=True, cmap=sns.light_palette("purple"))
    plt.title(key)


Peu importe le type d'association utilisée, les heatmaps montrent les mêmes tendances. Les assocaitions de dissimilarité (1ière ligne) s'étalent de 0 à 1, tandis que les distances (2ième ligne) partent de zéro, mais n'ont pas de limite supérieure. On note les plus grandes différences entre les sites 2 et 4, tandis que les sites 2 et 3 sont les plus semblables pour toutes les mesures d'association à l'exception de la dissimilarité de Kulczynski.

Objets: Occurence (présence-absence)

Des indices d'association différents devraient être utilisés lorsque des données sont compilées sous forme booléenne. En général, les tableaux de données d'occurence seront compilés avec des 1 (présence) et des 0 (absence).

La similarité de Jaccard entre le site A et le site B est la proportion de double 1 (présences de 1 dans A et B) parmi les espèces. La dissimilarié est la proportion complémentaire (comprenant [1, 0], [0, 1] et [0, 0]). La distance de Jaccard est la racine carrée de la dissimilarité.


In [17]:
associations_occ = {}
associations_occ['Jaccard'] = pdist(occurence, 'jaccard')

Les distances d'Hellinger, de chord et de chi-carré sont aussi appropriées pour les calculs de distances sur des tableaux d'occurence.


In [18]:
Hel_occurence = occurence.apply(lambda x: np.sqrt(x / np.sum(x)), axis=1)
associations_occ['Hellinger'] = pdist(Hel_occurence, 'euclidean')

Graphiquement,


In [19]:
plt.figure(figsize=(12, 5))
for i, key in enumerate(associations_occ):
    plt.subplot(1, 2, i+1)
    assoc_df = pd.DataFrame(squareform(associations_occ[key]), index=occurence.index, columns=occurence.index)
    sns.heatmap(assoc_df, annot=True, square=True, cmap=sns.light_palette("purple"))
    plt.title(key)


Il est attendu que les matrices d'association sur l'occurence sont semblables à celles sur l'abondance. Dans ce cas-ci, la distance d'Hellinger donne des résultats semblables à la dissimilarité de Jaccard.

Objets: Données quantitatives

Les données quantitative en écologie peuvent décrire l'état de l'environnement: le climat, l'hydrologie, l'hydrogéochimie, la pédologie, etc. En règle générale, les coordonnées des sites ne sot pas des variables environnementales, à que l'on soupçonne la coordonnée elle-même d'être responsable d'effets sur notre système: mais il s'agira la plupart du temps d'effets confondants (par exemple, on peut mesurer un effet de lattitude sur le rendement des agrumes, mais il s'agira probablement avant tout d'effets dus aux conditions climatiques, qui elles changent en fonction de la lattitude). D'autre types de données quantitative pouvant être appréhendées par des distances sont les traits phénologiques, les ionomes, les génomes, etc.

La distance euclidienne est la racine carrée de la somme des carrés des distances sur tous les axes. Il s'agit d'une application multidimensionnelle du théorème de Pythagore. La distance d'Aitchison, couverte dans le chapitre 6, est une distance euclidienne calculée sur des données compositionnelles préalablement transformées. La distance euclidienne est sensible aux unités utilisés: utiliser des milimètres plutôt que des mètres enflera la distance euclidienne. Il est recommandé de porter une attention particulière aux unités, et de standardiser les données au besoin (par exemple, en centrant la moyenne à zéro et en fixant l'écart-type à 1).

On pourrait, par exemple, mesurer la distance entre des observations des dimensions de différentes espèces d'iris. Ce tableau est inclu dans le module graphique seaborn.


In [20]:
import seaborn as sns
iris = sns.load_dataset("iris")
iris.sample(5)


Out[20]:
sepal_length sepal_width petal_length petal_width species
40 5.0 3.5 1.3 0.3 setosa
138 6.0 3.0 4.8 1.8 virginica
69 5.6 2.5 3.9 1.1 versicolor
146 6.3 2.5 5.0 1.9 virginica
54 6.5 2.8 4.6 1.5 versicolor

L'objet iris est un tableau pandas. Les mesures du tableau sont en centimètres. Pour éviter de donner davantage de poids aux longueur des sépales et en même temps de négliger la largeur des pétales, nous allons standardiser le tableau.


In [21]:
iris_features = iris.drop('species', axis=1)
iris_s = iris_features.apply(lambda x: (x-x.mean())/x.std(), axis=0)

La standadrisation semble correcte. Notez que la fonction apply_along_axis de numpy ressemble à fonction .axis de pandas.DataFrame, à l'exception de la manière de l'appeler, de l'ordre et du nom des arguments.


In [22]:
associations_cont = {}
associations_cont['euclidean'] = pdist(iris_s, 'euclidean')

La distance de Mahalanobis est semblable à la distance euclidienne, mais qui tient compte de la covariance de la matrice des objets. Cette covariance peut être utilisée pour décrire la structure d'un nuage de points. La figure suivante montre deux points verts qui se trouvent aux extrêmes d'un nuage de point. Ces points ont des distances euclidiennes par rapport au centre différentes: les lignes d'équidistance euclédienne sont tracées en rose. Toutefois, les deux points ont un distance de Mahalanobis égale à partir du centre.

Source: [Parent et al. (2012)](https://www.intechopen.com/books/soil-fertility/nutrient-balance-as-paradigm-of-plant-and-soil-chemometricsnutrient-balance-as-paradigm-of-soil-and-).

La distance de Mahalanobis permet de représenter des distances dans un espace fortement corrélé. Elle est courramment utilisée pour détecter les valeurs aberrantes selon des critères de distance à partir du centre d'un jeu de données multivariées.


In [23]:
associations_cont['Mahalanobis'] = pdist(iris_features, 'mahalanobis')

La distance de Manhattan porte aussi le nom de distance de cityblock ou de taxi. C'est la distance que vous devrez parcourir pour vous rendre du point A au point B à Manhattan, c'est-à-dire selon une séquence de tronçons perpendiculaires.

$$ D_{AB} = \sum _i \left| A_i - B_i \right| $$

La distance de Manhattan est appropriée lorsque les gradients (changements d'un état à l'autre ou d'une région à l'autre) ne permettent pas des changements simultanés. Mieux vaut standardiser les variables pour éviter qu'une dimension soit prépondérante.


In [24]:
associations_cont['Manhattan'] = pdist(iris_s, 'cityblock')

Lors de la création du graphique, l'argument annot=True demande beaucoup de temps de calcul. La valeur par défaut étant False, on peut seulement retirer l'argument ou le spécifier explicitement.


In [25]:
plt.figure(figsize=(20, 6))
for i, key in enumerate(associations_cont):
    plt.subplot(1, 3, i+1)
    assoc_df = pd.DataFrame(squareform(associations_cont[key]))
    sns.heatmap(assoc_df, annot=False, square=True, cmap=sns.light_palette("purple"))
    plt.title(key)


Le tableau iris est ordonné par espèce. Les distances euclidienne et deManhattan permettent aisément de distinguer les espèces selon les dimensions des pétales et des sépales. Toutefois, l'utilsation de la covariance avec la distance de Mahalanobis crée des distinction moins tranchées.

Objets: Données mixtes

Les données catégorielles ordinales peuvent être transformées en données continues par gradations linéaires ou quadratiques. Les données catégorielles nominales, quant à elles, peuvent être dummyfiées en données similaires à des occurences. Attention toutefois: contrairement à la régression linéaire qui demande d'exclure une catégorie, la dummyfication doit inclure toutes les catégories. Le comportement par défaut de la fonction pandas.get_dummies est de garder toutes les catégories. La similarité de Gower a été développée pour mesurer des associations entre des objets dont les données sont mixtes: booléennes, catégorielles et continues. La similarité de Gower est calculée en additionnant les distances calculées par colonne, individuellement. Si la colonne est booléenne, on utilise les distances de Jaccard (qui exclue les double-zéro) de manière univariée: une variable à la fois. Pour les variables continues, on utilise la distance de Manhattan divisée par la plage de valeurs de la variable (pour fin de standardisation). Puisqu'elle hérite de la particularité de la distance de Manhattan et de la similarité de Jaccard univariée, la similarité de Gower reste une combinaison linéaire de distances univariées.


In [26]:
# Notes
## possible avec ecopy, mais erreur. pull-request: https://github.com/Auerilas/ecopy/pulls
## possible avec scikit-learn: https://github.com/matchado/Misc/blob/master/gower_dist.py, mais jaccard non pas dice
## Attention: ça ne fonctionne pas bien. Changer la référence de dummification change la distance.
## Attendre la version de scikit-learn en développement: https://github.com/scikit-learn/scikit-learn/pull/9555

from sklearn.neighbors import DistanceMetric
def gower_distance(X):
    individual_variable_distances = []
    for i in range(X.shape[1]):
        feature = X.iloc[:,[i]]
        if feature.dtypes[0] == np.object:
            feature_dist = DistanceMetric.get_metric('jaccard').pairwise(feature)
        else:
            feature_dist = DistanceMetric.get_metric('manhattan').pairwise(feature) / np.ptp(feature.values)

        individual_variable_distances.append(feature_dist)

    return np.array(individual_variable_distances).mean(0)

Générons un jeu de données.


In [27]:
X = pd.DataFrame({'age':[21,21,19,30,21,21,19,30],
                'gender':['M','M','N','M','F','F','F','F'],
                'civil_status':['MARRIED','SINGLE','SINGLE','SINGLE','MARRIED','SINGLE','WIDOW','DIVORCED'],
                'salary':[3000.0,1200.0 ,32000.0,1800.0 ,2900.0 ,1100.0 ,10000.0,1500.0],
                'children':[True,False,True,True,True,False,False,True],
                'available_credit':[2200,100,22000,1100,2000,100,6000,2200]})

Il faut préalablement dummifier les variables catégorielles nominales.


In [28]:
X_dum = pd.get_dummies(X)
X_dum


Out[28]:
age available_credit children salary civil_status_DIVORCED civil_status_MARRIED civil_status_SINGLE civil_status_WIDOW gender_F gender_M gender_N
0 21 2200 True 3000.0 0 1 0 0 0 1 0
1 21 100 False 1200.0 0 0 1 0 0 1 0
2 19 22000 True 32000.0 0 0 1 0 0 0 1
3 30 1100 True 1800.0 0 0 1 0 0 1 0
4 21 2000 True 2900.0 0 1 0 0 1 0 0
5 21 100 False 1100.0 0 0 1 0 1 0 0
6 19 6000 False 10000.0 0 0 0 1 1 0 0
7 30 2200 True 1500.0 1 0 0 0 1 0 0

Calculons la dissimilarité de Gower.


In [29]:
d_gow = gower_distance(X_dum)

plt.figure(figsize=(6, 6))
sns.heatmap(d_gow, square=True, annot=True, cmap=sns.light_palette("purple"))


Out[29]:
<matplotlib.axes._subplots.AxesSubplot at 0x7feaf0085748>

Une avenue possible en mode multivarié est d'effectuer une moyenne (pondérée par une décision arbitraire) d'associations spécifiques.

Associations entre variables (mode R)

Il existe de nombreuses approches pour mesurer les associations entre variables. La plus connue est la corrélation. Mais les données d'abondance et d'occurence demandent des approches différentes.

Variables: Abondance

La distance du chi-carré est suggérée par Borcard et al. (2011).


In [30]:
D_chisq_R = squareform(chi_square(abundance.T, method='distance'))
D_chisq_R_df = pd.DataFrame(D_chisq_R, 
                         columns=abundance.columns, 
                         index=abundance.columns)

plt.figure(figsize=(6, 6))
sns.heatmap(D_chisq_R_df, square=True, annot=True, cmap=sns.light_palette("purple"))


Out[30]:
<matplotlib.axes._subplots.AxesSubplot at 0x7feaf002f470>

Variables: Occurence

La dissimilarité de Jaccard peut être utilisée.


In [31]:
D_jacc_R = squareform(pdist(abundance.T, metric='jaccard'))
D_jacc_R_df = pd.DataFrame(D_jacc_R, 
                         columns=abundance.columns, 
                         index=abundance.columns)

plt.figure(figsize=(6, 6))
sns.heatmap(D_jacc_R_df, square=True, annot=True, cmap=sns.light_palette("purple"))


Out[31]:
<matplotlib.axes._subplots.AxesSubplot at 0x7feadbea92e8>

Variables: Quantités

La matrice des corrélations de Pearson peut être utilisée pour les données continues. Quant aux variables ordinales, elles devraient idéalement être liées linéairement ou quadratiquement. Si ce n'est pas le cas, c'est-à-dire que les catégories sont ordonnées par rang seulement, vous pourrez avoir recours aux coefficients de corrélation de Spearman ou de Kendall.


In [32]:
iris_cor = pd.DataFrame(np.corrcoef(iris_features.T),
                        columns = iris_features.columns, 
                        index = iris_features.columns)

plt.figure(figsize=(6, 6))
sns.heatmap(iris_cor, annot=True, square=True, cmap=sns.light_palette("purple"))


Out[32]:
<matplotlib.axes._subplots.AxesSubplot at 0x7feaf28a5080>

Conclusion sur les associations

Il n'existe pas de règle claire pour déterminer quelle technique d'association utiliser. Cela dépend en premier lieu de vos données. Vous sélectionnerez votre méthode d'association selon le type de données que vous abordez, la question à laquelle vous désirez répondre ainsi l'expérience dans la littérature comme celle de vos collègues scientifiques. S'il n'existe pas de règle clair, c'est qu'il existe des dizaines de méthodes différentes, et la plupart d'entre elles vous donneront une perspective juste et valide. Il faut néanmoins faire attention pour éviter de sélectionner les méthodes qui ne sont pas appropriées.

Partitionnement

Les données suivantes ont été générées par Leland McInnes (Tutte institute of mathematics, Ottawa). Êtes-vous en mesure d'identifier des groupes? Combien en trouvez-vous?


In [33]:
data_mcinnes = np.load('data/clusterable_data.npy')
df_mcinnes = pd.DataFrame({'x': data_mcinnes[:, 0],
                          'y': data_mcinnes[:, 1]})
plt.figure(figsize=(8, 8))
sns.lmplot('x', 'y', data=df_mcinnes, fit_reg=False, markers='.', size=6, aspect=1)


Out[33]:
<seaborn.axisgrid.FacetGrid at 0x7feadbf892e8>
<matplotlib.figure.Figure at 0x7feaf3315400>

En 2D, l'oeil humain peut facilement détecter les groupes. En 3D, c'est toujours possible, mais au-delà de 3D, le partitionnement cognitive devient rapidement maladroite. Les algorithmes sont alors d'une aide précieuse. Mais ils transportent en pratique tout un baggage de limitations. Quel est le critère d'association entre les groupes? Combien de groupe devrions-nous créer? Comment distinguer une donnée trop bruitée pour être classifiée?

Le partitionnement de données (clustering en anglais), et inversement leur regroupement, permet de créer des ensembles selon des critères d'association. On suppose donc que Le partitionnement permet de créer des groupes selon l'information que l'on fait émerger des données. Il est conséquemment entendu que les données ne sont pas catégorisées à priori: il ne s'agit pas de prédire la catégorie d'un objet, mais bien de créer des catégories à partir des objets par exemple selon leurs dimensions, leurs couleurs, leurs signature chimique, leurs comportements, leurs gènes, etc.

Plusieurs méthodes sont aujourd'hui offertes aux analystes pour partitionner leurs données. Dans le cadre de ce manuel, nous couvrirons ici deux grandes tendances dans les algorithmes.

  1. Méthodes hiérarchique et non hiérarchiques. Dans un partitionnement hiérarchique, l'ensemble des objets forme un groupe, comprenant des sous-regroupements, des sous-sous-regroupements, etc., dont les objets forment l'ultime partitionnement. On pourra alors identifier comment se décline un partitionnement. À l'inverse, un partitionnement non-hiérarchique des algorhitmes permettent de créer les groupes non hiérarchisés les plus différents que possible.

  2. Membership exclusif ou flou. Certaines techniques attribuent à chaque une classe unique: l'appartenance sera indiquée par un 1 et la non appartenance par un 0. D'autres techniques vont attribuer un membership flou où le degré d'appartenance est une variable continue de 0 à 1. Parmi les méthodes floues, on retrouve les méthodes probabilistes.

Évaluation d'un partitionnement

Le choix d'une technique de partitionnement parmi de nombreuses disponibles, ainsi que le choix des paramètres gouvernant chacune d'entre elles, est avant tout basé sur ce que l'on désire définir comme étant un groupe, ainsi que la manière d'interpréter les groupes. En outre, le nombre de groupe à départager est toujours une décision de l'analyste. Néanmoins, on peut se fier des indicateurs de performance de partitionnement. Parmis ceux-ci, retenons le score silouhette ainsi que l'indice de Calinski-Harabaz.

Score silouhette

En anglais, le h dans silouhette se trouve après le l: on parle donc de silhouette coefficient pour désigner le score de chacun des objets dans le partitionnement. Pour chaque objet, on calcule la distance moyenne qui le sépare des autres points de son groupe ($a$) ainsi que la distance moyenne qui le sépare des points du groupe le plus rapproché.

$$s = \frac{b-a}{max \left(a, b \right)}$$

Un coefficient de -1 indique le pire classement, tandis qu'un coefficient de 1 indique le meilleur classement. La moyenne des coefficients silouhette est le score silouhette.

Indice de Calinski-Harabaz

L'indice de Calinski-Harabaz est proportionnel au ratio des dispersions intra-groupe et la moyenne des dispersions inter-groupes. Plus l'indice est élevé, mieux les groupes sont définis. La mathématique est décrite dans la documentation de scikit-learn.

Note. Les coefficients silouhette et l'indice de Calinski-Harabaz sont plus appropriés pour les formes de groupes convexes (cercles, sphères, hypersphères) que pour les formes irrégulières (notamment celles obtenues par la DBSCAN, discutée ci-desssous).

Partitionnement non hiérarchique

Il peut arriver que vous n'ayez pas besoin de comprendre la structure d'agglomération des objets (ou variables). Plusieurs techniques de partitionnement non hiérarchique sont disponibles dans le module scikit-learn. On s'intéressera en particulier à celles-ci.

Kmeans (sklearn.cluster.Kmeans). L'objectif des kmeans est de minimiser la distance euclédienne entre un nombre prédéfini de k groupes exclusifs.

  1. L'algorhitme commence par placer une nombre k de centroides au hasard dans l'espace d'un nombre p de variables (vous devez fixer k, et p est le nombre de colonnes de vos données).
  2. Ensuite, chaque objet est étiquetté comme appartenant au groupe du centroid le plus près.
  3. La position du centroide est déplacée à la moyenne de chaque groupe.
  4. Recommencer à partir de l'étape 2 jusqu'à ce que l'assignation des objets aux groupes ne change plus.

Source: [David Sheehan](https://dashee87.github.io/data%20science/general/Clustering-with-Scikit-with-GIFs/)

La technique des kmeans suppose que les groupes ont des distributions multinormales - représentées par des cercles en 2D, des sphères en 3D, des hypersphères en plus de 3D. Cette limitation est problématique lorsque les groupes se présentent sous des formes irrégulières, comme celles du nuage de points de Leland McInnes, présenté plus haut. De plus, la technique classique des kmeans est basée sur des distances euclidiennes: l'utilisation des kmeans n'est appropriée pour les données comprenant beaucoup de zéros, comme les données d'abondance, qui devraient préalablement être transformées en variables centrées et réduites (Legendre et Legendre, 2012). La technique des mixtures gaussiennes (gaussian mixtures, sklearn.mixture.GaussianMixture) est une généralisation des kmeans permettant d'intégrer la covariance des groupes. Les groupes ne sont plus des hyper-sphères, mais des hyper-ellipsoïdes.

DBSCAN (sklearn.cluster.DBSCAN). La technique DBSCAN ( Density-Based Spatial Clustering of Applications with Noise) sousentend que les groupes sont composés de zones où l'on retrouve plus de points (zones denses) séparées par des zones de faible densité. Pour lancer l'algorithme, nous devons spécifier une mesure d'association critique (distance ou dissimilarité) d ainsi qu'un nombre de point critique k dans le voisinage de cette distance.

  1. L'algorithme comme étiqueter chaque point selon l'une de ces catégories:

    • Noyau: le point a au moins k points dans son voisinage, c'est-à-dire à une distance inférieure ou égale à d.
    • Bordure: le point a moins de k points dans son voisinage, mais l'un de des points voisins est un noyau.
    • Bruit: le cas échéant. Ces points sont considérés comme des outliers.

  2. Les noyaux distancés de d ou moins sont connectés entre eux en englobant les bordures.

Le nombre de groupes est prescrit par l'algorithme DBSCAN, qui permet du coup de détecter des données trop bruitées pour être classées.

Damiani et al. (2014) a développé une approche utilisant la technique DBSCAN pour partitionner des zones d'escale pour les flux de populations migratoires.

Application

Dans scikit-learn, on définit d'abord le modèle (par exemple Kmeans(...)), puis on l'applique à nos données (fit(...)), enfin on applique le modèle sur des données (predict(...)). Certaines fonctions utilisent toutefois le raccourcis fit_predict. Chaque algorithme doit être ajusté avec les paramètres qui convient. De nombreux paramètres par défaut sont utilisés dans les exécutions ci-dessous. Lors de travaux de recherche, l'utilsation d'un argument ou d'un autre dans une fonction doit être justifié: qu'un paramètre soit utilisé par défaut dans une fonction n'est a priori pas une justification convainquante.

Nous allons classifier les données de Leland McInnes, présentée ci-dessus. Importons les fonctions nécessaires.


In [34]:
from sklearn.cluster import KMeans, DBSCAN
from sklearn.mixture import GaussianMixture
from sklearn import metrics

Pour les kmeans, on doit fixer le nombre de groupes. Le graphique des données de Leland McInnes montrent 6 groupes. Toutefois, il est rare que l'on puisse visualiser de telles démarquations, qui plus est dans des cas où l'on doit traiter de plus de deux dimensions. Je vais donc lancer le partitionnement en boucle pour plusieurs nombres de groupes, de 3 à 10 et pour chaque groupe, évaluer le score silouhette et de Calinski-Habaraz. J'utilise un argument random_state pour m'assurer que les groupes seront les mêmes à chaque fois que la cellule sera lancée.


In [35]:
silouhette_score = []
calinski_index = []
n_groups = list(range(3, 10))
for i in n_groups:
    mcinnes_kmeans_labels = KMeans(n_clusters=i, random_state=5437925).fit_predict(data_mcinnes)
    silouhette_score.append(metrics.silhouette_score(data_mcinnes, mcinnes_kmeans_labels, 'euclidean'))
    calinski_index.append(metrics.calinski_harabaz_score(data_mcinnes, mcinnes_kmeans_labels))

plt.figure(figsize=(20, 6))
plt.subplot(1, 2, 1)
plt.plot(n_groups, silouhette_score, '-o')
plt.xlabel('Number of groups')
plt.ylabel('Silouhette score')
plt.subplot(1, 2, 2)
plt.plot(n_groups, calinski_index, '-o')
plt.xlabel('Number of groups')
plt.ylabel('Calinski-Harabaz index')


Out[35]:
Text(0,0.5,'Calinski-Harabaz index')

La forme des groupes n'étant pas convexe, il fallait s'attendre à ce que indicateurs maximaux pour les deux indicateurs soient différents. C'est d'ailleurs souvent le cas. Cet exemple supporte que le choix du nombre de groupe à départager repose sur l'analyste, non pas uniquement sur les indicateurs de performance. Choisissons 6 groupes, puisque que c'est visuellement ce que l'on devrait chercher pour ce cas d'étude.


In [36]:
from scikitplot import clustering_factory

mcinnes_kmeans = clustering_factory(KMeans(n_clusters=6, random_state=5437925))
df_mcinnes['group_kmeans'] = mcinnes_kmeans.fit_predict(data_mcinnes)
print("Score silouhette:",
      metrics.silhouette_score(data_mcinnes, df_mcinnes['group_kmeans'], 'euclidean'))
print("Indice de Calinski-Habaraz:",
      metrics.calinski_harabaz_score(data_mcinnes, df_mcinnes['group_kmeans']))


/home/essicolo/bin/anaconda3/lib/python3.6/site-packages/sklearn/utils/deprecation.py:70: DeprecationWarning: Function clustering_factory is deprecated; This will be removed in v0.4.0. The Factory API has been deprecated. Please migrate existing code into the various new modules of the Functions API. Please note that the interface of those functions will likely be different from that of the Factory API.
  warnings.warn(msg, category=DeprecationWarning)
Score silouhette: 0.404590768399
Indice de Calinski-Habaraz: 2146.05914177

Remarquez que la ligne mcinnes_kmeans = ... comprend la fonction clustering_factory. Cette fonction issue du module scikit-plot implémente Kmeans en lui dotant la possibilité de générer le graphique silouhette, qui présente les coefficients silouhette pour chacun des objets.


In [37]:
mcinnes_kmeans.plot_silhouette(data_mcinnes, figsize=(10, 6))


/home/essicolo/bin/anaconda3/lib/python3.6/site-packages/sklearn/utils/deprecation.py:70: DeprecationWarning: Function plot_silhouette is deprecated; This will be removed in v0.4.0. Please use scikitplot.metrics.plot_silhouette instead.
  warnings.warn(msg, category=DeprecationWarning)
Out[37]:
<matplotlib.axes._subplots.AxesSubplot at 0x7feaf3124ba8>

La technique DBSCAN n'est pas basée sur le nombre de groupe, mais sur la densité des points. Il est possible de spécifier le type d'association à utiliser dans l'argument metric, mais dans ce cas on doit se limiter aux métriques offertes par scikit-learn (qui sont celles de scipy). Si la métrique souhaitée n'est pas disponible parmi celles-là, la matrice d'association peut être spécifiée dans l'argument X plutôt que les données directe quand l'argument metric est égal à "precomputed". L'argument min_samples spécifie le nombre minimal de points qui l'on doit retrouver à une distance critique d pour la formation des noyaux et la propagation des groupes, spécifiée dans l'argument eps. La distance d peut être estimée en prenant une fraction de la moyenne, mais on aura volontiers recours à sont bon jugement.


In [38]:
mcinnes_dbscan_distmat = squareform(pdist(data_mcinnes, 'euclidean'))
mcinnes_dbscan = DBSCAN(metric='precomputed', min_samples=6, eps=0.026)
df_mcinnes['group_dbscan'] = mcinnes_dbscan.fit_predict(mcinnes_dbscan_distmat)
df_mcinnes['group_dbscan'].unique()


Out[38]:
array([ 0, -1,  1,  2,  3,  4,  5,  6,  7])

Les paramètres spécifiés donnent 8 groupes (0, 1, 2, ..., 7) et des points trop bruités pour être classifiés (étiquetés -1). Voyons comment les groupes ont été formés.


In [39]:
df_mcinnes_tidy = pd.melt(df_mcinnes, id_vars=['x', 'y'])
sns.lmplot('x', 'y', data=df_mcinnes_tidy, 
           fit_reg=False, 
           hue='value', markers='.', size=6, 
           aspect=1,
           col='variable')


Out[39]:
<seaborn.axisgrid.FacetGrid at 0x7feaf2c53518>

L'algorithme kmeans est loin d'être statisfaisant. Cela est attendu, puisque les kmeans recherchent des distribution gaussiennes sur des groupes vraisemblablement non-gaussiens. Quant à DBSCAN, le partitionnement semble plus conforme à ce que l'on recherche. Néanmoins, DBSCAN cré quelques petits groupes indésirables ainsi que des grands groupes (violet et vert) qui auraient lieu d'être partitionnés.

Partitionnement hiérarchique

Les techniques de partitionnement hiérarchique sont basées sur les matrices d'association. La technique pour mesurer l'association (entre objets ou variables) déterminera en grande partie le paritionnement des données. Les partitionnements hiérarchiques ont l'avantage de pouvoir être représentés sous forme de dendrogramme (ou arbre) de partition. Un tel dendrogramme présente des sous-groupes qui se joignent en groupes jusqu'à former un seul ensemble.

Le partitionnement hiérarchique est abondamment utilisé en phylogénie, pour étudier les relations de parenté entre organismes vivants, populations d'organismes et espèces. La phénétique, branche empirique de la phylogénèse interspécifique, fait usage du partitionnement hiérarchique à partir d'associations génétiques entre unités taxonomiques. On retrouve de nombreuses ressources académiques en phylogénétique ainsi que des outils pour Python et R. Toutefois, la phylogénétique en particulier ne fait pas partie de la présente ittération de ce manuel.

Techniques de partitionnement hiérarchique

Le partitionnement hiérarchique est typiquement effectué avec une des quatres méthodes suivantes, dont chacune possède ses particularités, mais sont toutes agglomératives: à chaque étape d'agglomération, on fusionne les deux groupes ayant le plus d'affinité sur la base des deux sous-groupes les plus rapprochés.

Single link (single). Les groupes sont agglomérés sur la base des deux points parmi les groupes, qui sont les plus proches.

Complete link (complete). À la différence de la méthode single, on considère comme critère d'agglomération les éléments les plus éloignés de chaque groupe.

Agglomération centrale. Il s'agit d'une fammille de méthode basées sur les différences entre les tendances centrales des objets ou des groupes.

  • Average (average). Appelée UPGMA (Unweighted Pair-Group Method unsing Average), les groupes sont agglomérés selon un centre calculés par la moyenne et le nombre d'objet pondère l'agglomération (le poids des groupes est retiré). Cette technique est historiquement utilisée en bioinformatique pour partitionner des groupes phylogénétiques (Sneath et Sokal, 1973).
  • Weighted (weighted). La version de average, mais non pondérée (WPGMA).
  • Centroid (centroid). Tout comme average, mais le centroïde (centre géométrique) est utilisé au lieu de la moyenne. Accronyme: UPGMC.
  • Median (median). Appelée WPGMC. Devinez! ;)

Ward (ward). L'optimisation vise à minimiser les sommes des carrés par regroupement.

Quel outil de partitionnement hiérarchique utiliser?

Alors que le choix de la matrice d'association dépend des données et de leur contexte, la technique de partitionnement hiérarchique peut, quant à elle, être basée sur un critère numérique. Il en existe plusieurs, mais le critère recommandé pour le choix d'une technique de partitionnement hiérarchique est la corrélation cophénétique. La distance cophénétique est la distance à laquelle deux objets ou deux sous-groupes deviennent membres d'un même groupe. La corrélation cophénétique est la corrélation de Pearson entre le vecteur d'association des objets et le vecteur de distances cophénétiques.

Application

Les techniques de partitionnement hiérarchique présentées ci-dessus, ainsi que la fonction cophenet pour calculer la corrélation cophénétique, sont disponibles dans scipy.cluster.hierarchy. Nous allons classifier les dimensions des iris grâce à la distance de Manhattan.


In [40]:
mcinnes_hclust_distmat = pdist(data_mcinnes, 'euclidean')

In [41]:
from scipy.cluster.hierarchy import *

clustering_methods = ['single', 'complete', 'average', 'weighted', 'centroid', 'median', 'ward']
clust_l = []
coph_corr_l = []
for i in range(len(clustering_methods)):
    clust_m = linkage(y=mcinnes_hclust_distmat, method=clustering_methods[i])
    clust_l.append(clust_m)
    coph_corr_l.append(cophenet(clust_m, mcinnes_hclust_distmat)[0])

# créer un dictionnaire avec les listes
clust_d = dict(zip(clustering_methods, clust_l))
coph_corr_d = dict(zip(clustering_methods, coph_corr_l))

# créer un tableau grâce au dictionnaire
coph_corr_df = pd.DataFrame(coph_corr_d, index=[0]).melt()
coph_corr_df


Out[41]:
variable value
0 average 0.697121
1 centroid 0.713944
2 complete 0.534390
3 median 0.516069
4 single 0.486935
5 ward 0.584880
6 weighted 0.557286

La méthode centroid retourne la corrélation la plus élevée. Pour cet exemple, nous utiliserons average afin de faciliter l'intégration avec scikit-learn, qui ne comprend pas dans sa version actuelle le partitionnement par centroid dans ses fonctions. Pour plus de flexibilité, enchâssons le nom de la méthode dans une variable. Ainsi, en chageant le nom de cette variable, le reste du code sera conséquent.


In [42]:
best_method = 'average'

Le partitionnement hiérarchique peut être visualisé par un dendrogramme.


In [43]:
plt.figure(figsize=(20, 6))
dnd = dendrogram(clust_d[best_method], color_threshold=0, leaf_font_size=8)
plt.xlabel('Sample')
plt.ylabel('Height')


Out[43]:
Text(0,0.5,'Height')

Combien de groupes utiliser?

La longueur des lignes verticales est la distance séparant les groupes enfants. Bien que la sélection du nombre de groupe soit avant tout basée sur les besoins du problème, nous pouvons nous appuyer sur certains outils. La hauteur totale peut servir de critère pour définir un nombre de groupes adéquat. On pourra sélectionner le nombre de groupe où la hauteur se stabilise en fonction du nombre de groupe. On pourra aussi utiliser le graphique silhouette, comprenant une collection de largeurs de silouhette, représentant le degré d'appartenance à son groupe. La fonction sklearn.metrics.silhouette_score, du module scikit-learn, s'en occupe.


In [44]:
nodes_f = fcluster(Z=clust_d[best_method], t=2, criterion='maxclust')
np.unique(nodes_f)


Out[44]:
array([1, 2], dtype=int32)

In [45]:
from sklearn import metrics

# hauteur
heights = clust_d[best_method][:, 2]

# silouhette
silouhette = []
max_clus = 15
n_clust = list(range(3, max_clus))
for i in n_clust:
    nodes_f = fcluster(Z=clust_d[best_method], t=i, criterion='maxclust')
    silouhette.append(metrics.silhouette_score(data_mcinnes, nodes_f, metric='euclidean'))

# nombre de groupes correspondant à la largeur de silhouette maximale
n_clust_smax = np.array(n_clust)[silouhette == np.max(silouhette)]
silouhette_smax = np.array(silouhette)[silouhette == np.max(silouhette)]

In [46]:
n_clust_smax


Out[46]:
array([14])

La largeur de silhouette maximale est atteinte à seulement deux groupes. Sur graphique,


In [47]:
f, ax = plt.subplots(2, sharex=True, figsize=(12, 6))

ax[0].plot(list(range(len(heights)))[::-1], heights, '-o')
ax[0].set_ylabel('Node height')
ax[0].axvline(x=n_clust_smax, c='k', linewidth=1)

ax[1].plot(n_clust, silouhette, '-o')
ax[1].axvline(x=n_clust_smax, c='k', linewidth=1)
ax[1].plot(n_clust_smax, silouhette_smax, 'o')
ax[1].set_xlabel('Number of clusters')
ax[1].set_ylabel('Silouhette score')
ax[1].set_xlim([0, max_clus])


Out[47]:
(0, 15)

Un nombre élevé de groupe pourrait être créés. Mais puisque nous pouvons discerner 6 groupes dans les données, retenons 6 groupes, coupés à une hauteur entre 0.35 et 0.40 sur le dendrogramme.


In [48]:
heights_test = np.arange(0.35, 0.40, 0.005)
n_groups = []
for i in heights_test:
    nodes_f = fcluster(Z=clust_d[best_method], t=i, criterion='distance')
    n_groups.append(len(np.unique(nodes_f)))
plt.plot(n_groups, heights_test, '-o')
plt.xlabel('Nombre de groupes')
plt.ylabel('Hauteur')


Out[48]:
Text(0,0.5,'Hauteur')

Il n'y a qu'une seule valeur de distance calculé pouvant générer 6 groupes.


In [49]:
heights_test[np.array(n_groups)==6][0]


Out[49]:
0.375

Coupons le dendrorgamme à cette hauteur.


In [50]:
cut_off = heights_test[np.array(n_groups)==6][0]
plt.figure(figsize=(15, 5))
dendrogram(clust_d[best_method], leaf_rotation=90, color_threshold=cut_off,
          above_threshold_color='black')
plt.axhline(y=cut_off, c='grey', linestyle='dashed')
plt.ylabel('Height')


Out[50]:
Text(0,0.5,'Height')

On voit clairement qu'un le cluster comporte peu de données. On fera intervenir scikit-learn pour créer un graphique silouhette.


In [51]:
from sklearn.cluster import AgglomerativeClustering
mcinnes_hclust_sk = clustering_factory(AgglomerativeClustering(n_clusters=6, affinity='euclidean',
                                                               linkage='average').fit(data_mcinnes))
mcinnes_hclust_sk.plot_silhouette(data_mcinnes, figsize=(10, 6))


/home/essicolo/bin/anaconda3/lib/python3.6/site-packages/sklearn/utils/deprecation.py:70: DeprecationWarning: Function clustering_factory is deprecated; This will be removed in v0.4.0. The Factory API has been deprecated. Please migrate existing code into the various new modules of the Functions API. Please note that the interface of those functions will likely be different from that of the Factory API.
  warnings.warn(msg, category=DeprecationWarning)
/home/essicolo/bin/anaconda3/lib/python3.6/site-packages/sklearn/utils/deprecation.py:70: DeprecationWarning: Function plot_silhouette is deprecated; This will be removed in v0.4.0. Please use scikitplot.metrics.plot_silhouette instead.
  warnings.warn(msg, category=DeprecationWarning)
Out[51]:
<matplotlib.axes._subplots.AxesSubplot at 0x7feacfe26ef0>

Les groupes 3 et 4 comportent peut d'objets. Les groupes 0 et 1 comportent un bon nombre d'objets apparamment mal classés.

Une autre forme de représentation de le partitionnement hiérarchique est le clustermap, qui présente le partitionnement ainsi que la matrice de données.


In [52]:
sns.clustermap(data_mcinnes, metric='euclidean', method=best_method, cmap=sns.light_palette("purple"))


Out[52]:
<seaborn.matrix.ClusterGrid at 0x7feacfddbd30>

Intégrons nos résultats dans le tableau.


In [53]:
df_mcinnes['group_hclust'] = fcluster(Z=clust_d[best_method], t=6, criterion='maxclust')

Partitionnement hiérarchique basée sur la densité des points

La tecchinque HDBSCAN, dont l'algorithme est relativement récent (Campello et al., 2013), permet une partitionnement hiérarchique sur le même principe des zones de densité de la technique DBSCAN. Le HDBSCAN a été utilisée pour partitionner les lieux d'escale d'oiseaux migrateurs en Chine (Xu et al., 2013).

Avec DBSCAN, un rayon est fixé dans une métrique appropriée. Pour chaque point, on compte le nombre de point voisins, c'est à dire le nombre de point se situant à une distance (ou une dissimilarité) égale ou inférieure au rayon fixé. Avec HDBSCAN, on spécifie le nombre de points devant être recouverts et on calcule le rayon nécessaire pour les recouvrir. Ainsi, chaque point est associé à un rayon critique que l'on nommera $d_{noyau}$. La métrique initiale est ensuite altérée: on remplace les associations entre deux objets A et B par la valeur maximale entre cette association, le rayon critique de A et le rayon critique de B. Cette nouvelle distance est appelée la distance d'atteinte mutuelle: elle accentue les distances pour les points se trouvant dans des zones peu denses. On applique par la suite un algorithme semblable à la partition hiérarchique single link: En s'élargissant, les rayons se superposent, chaque superposition de rayon forment graduellement des groupes qui s'agglomèrent ainsi de manière hiérarchique. Au lieu d'effectuer une tranche à une hauteur donnée dans un dendrogramme de partitionnement, la technique HDBSCAN se base sur un dendrogramme condensé qui discarte les sous-groupes comprenant moins de n objets ($n_{gr min}$). Dans nouveau dendrogramme, on recherche des groupes qui occupent bien l'espace d'analyse. Pour ce faitre, on utilise l'inverse de la distance pour créer un indicateur de persistance (semblable à la similarité), $\lambda$. Pour chaque groupe hiérarchique dans le dendrogramme condensé, on peut calculer la persistance où le groupe prend naissance. De plus, pour chaque objet d'un groupe, on peut aussi calculer une distance à laquelle il quitte le groupe. La stabilité d'un groupe est la domme des différences de persistance entre la persistance à la naissance et les persistances des objets. On descend dans le dendrogramme. Si la somme des stabilité des groupes enfants est plus grande que la stabilité du groupe parent, on accepte la division. Sinon, le parent forme le groupe. La documentation du module hdbscan pour Python offre une description intuitive et plus exhaustive des principes et algorithme de HDBSCAN.

Paramètres

Outre la métrique d'association dont nous avons discuté, HDBSCAN demande d'être nourri avec quelques paramètres importants. En particulier, le nombre minimum d'objets par groupe, $n_{gr min}$ (ou min_cluster_size dans le module hdbscan) dépend de la quantité de données que vous avez à votre disposition, ainsi que de la quantité d'objets que vous jugez suffisante pour créer des groupes. Ce paramètre a évedemment une influence sur le nombre de groupes formés, mais aussi sur la quantité de valeurs aberrantes détectées. Un autre paramètre important est le nombre de points devant être recouverts, $k$ (min_samples dans hdbscan). Plus $k$ est élevé, plus la densité nécessaire pour créer un groupe sera élevée et plus HDBSCAN désignera de points somme du bruit. Parmi les algorithmes de sélection des groupes du module hdbscan, on retrouve l'excess of mass (cluster_selection_method='eom', par défaut) qui a tendance à créer des groupes de tailles différentes, et leaf (cluster_selection_method='leaf') qui permet de créer des groupes comprenant des quantités d'objets plus homogènes.

Application

Utilisons le module hdbscan ainsi que les données de Leland McInnes. La métrique utilisée sera euclidienne. Les métriques disponibles sont celles de scikit-learn. Il est néamoins possible d'utiliser une métrique personalisée (Gower, par exemple) en utilisant la matrice d'association en guise de données et en spécifiant metric='precomputed'.


In [54]:
import hdbscan
hdbscan_clust = hdbscan.HDBSCAN(min_cluster_size=20, # le nombre de point minimum par groupe
                                min_samples=5, # le nombre de points dans le voisinage pour calculer le "rayon noyau"
                                metric='euclidean', 
                                gen_min_span_tree=True)
hdbscan_clust.fit(data_mcinnes)
df_mcinnes['group_hdbscan'] = hdbscan_clust.labels_

Demandons les valeurs uniques de la nouvelle colonne group_hdbscan de notre tableau pour obtenir le nombre de groupes calculés.


In [55]:
df_mcinnes['group_hdbscan'].unique()


Out[55]:
array([ 1, -1,  2,  0,  3,  5,  4])

Nous avons 6 groupes, numérotés de 0 à 5, ainsi que des étiquettes identifiant des objets désignés comme étant du bruit de fond, numéroté -1. Le dendrogramme non condensé peu être produit.


In [56]:
plt.figure(figsize=(15, 5))
hdbscan_clust.single_linkage_tree_.plot(cmap='viridis', colorbar=True)


Out[56]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fead0138278>