Classification sous Python

Nous utilisons encore le module scikit-learn, dédié au Machine Learning.

import numpy
import pandas
import matplotlib.pyplot as plt
import seaborn
seaborn.set_style("white") # change le style par défaut des graphiques seaborn

%matplotlib inline

Pour tester les méthodes de classification, nous allons continuer l'étude de cas sur les données WGI.

# le dropna() permet de supprimer les pays pour lesquels il manque des informations
WGI_complet = pandas.read_csv("https://fxjollois.github.io/donnees/WGI/wgi2019.csv").dropna()
WGI_complet

	Country	Code	Voice and Accountability	Political Stability and Absence of Violence/Terrorism	Government Effectiveness	Regulatory Quality	Rule of Law	Control of Corruption
0	Aruba	ABW	1.294189	1.357372	1.029933	0.857360	1.263128	1.217238
1	Andorra	ADO	1.139154	1.615139	1.908749	1.228176	1.579939	1.234392
2	Afghanistan	AFG	-0.988032	-2.649407	-1.463875	-1.120555	-1.713527	-1.401076
3	Angola	AGO	-0.777283	-0.311101	-1.117144	-0.893871	-1.054343	-1.054683
5	Albania	ALB	0.151805	0.118570	-0.061331	0.274380	-0.411179	-0.528758
...	...	...	...	...	...	...	...	...
209	Serbia	SRB	0.026626	-0.091665	0.019079	0.113867	-0.119070	-0.445551
210	South Africa	ZAF	0.670388	-0.217931	0.367380	0.156172	-0.076408	0.084924
211	Congo, Dem. Rep.	ZAR	-1.365966	-1.808007	-1.627429	-1.509667	-1.786088	-1.538931
212	Zambia	ZMB	-0.286199	-0.102216	-0.675215	-0.554269	-0.462069	-0.640345
213	Zimbabwe	ZWE	-1.141875	-0.920179	-1.205337	-1.463199	-1.257009	-1.238796

202 rows × 8 columns

Classification

Nous allons ici réaliser une classification non supervisée avec 3 méthodes (CAH, $k$-means et DBSCAN). Bien évidemment, ici, nous pouvons peut-être supposé qu'une partition en 6 classes (1 par continent) est possible. Mais nous n'avons pas de réelle idée de ce nombre de classes. Pour l'estimer, nous allons utiliser 2 techniques.

Classification Ascendante Hiérarchique (CAH)

La première possibilité est d'utiliser la CAH, afin d'obtenir un dendrogramme représentant les regroupements 2 à 2 de toutes les classes.

Tout d'abord, nous importons la méthode, puis nous l'appliquons. Ici, on indique qu'on ne souhaite pas réaliser de découpage au préalable, et donc qu'on souhaite garder toute la hiérarchie. Voici pourquoi on choisit 0 pour le paramètre distance_threshold et None pour n_clusters. Le paramètre affinity permet de choisir la distance utilisée (euclidean par défaut - ce qui nous va bien). Le paramètre linkage permet lui de choisir le critère d'aggrégation entre deux classes (ward par défaut - ce qui nous va bien aussi).

from sklearn.cluster import AgglomerativeClustering

WGI_num = WGI_complet.drop(columns = ["Country", "Code"])
hac = AgglomerativeClustering(distance_threshold = 0, n_clusters = None)
hac.fit(WGI_num)

AgglomerativeClustering(distance_threshold=0, n_clusters=None)

Cette page décrit le programme permettant d'afficher le dendrogramme à partir du résultat de la fonction ci-dessous. Je recopie ici que la fonction (avec un peu d'adaptation).

from scipy.cluster.hierarchy import dendrogram

def plot_dendrogram(model, **kwargs):
    # Create linkage matrix and then plot the dendrogram

    # create the counts of samples under each node
    counts = numpy.zeros(model.children_.shape[0])
    n_samples = len(model.labels_)
    for i, merge in enumerate(model.children_):
        current_count = 0
        for child_idx in merge:
            if child_idx < n_samples:
                current_count += 1  # leaf node
            else:
                current_count += counts[child_idx - n_samples]
        counts[i] = current_count

    linkage_matrix = numpy.column_stack([model.children_, model.distances_, counts]).astype(float)

    # Plot the corresponding dendrogram
    dendrogram(linkage_matrix, **kwargs)

Nous utilisons donc cette méthode sur nos pays. On remarque un saut important entre 2 et 3 classes, puis entre 3 et 4 classes (et un peu moindre entre 4 et 5 classes). On imagine donc qu'on pourrait couper la hiérarchie à ces 2 niveaux (cf graphique). Ainsi, cela suggère qu'on pourrait avoir une partition soit en 2 classes, soit en 3 classes.

plt.figure(figsize = (16, 6))
plt.title("CAH (Ward) sur les pays")
# plot the top three levels of the dendrogram
plot_dendrogram(hac)
plt.axhline(y = 27, linewidth = .5, color = "dimgray", linestyle = "--")
plt.axhline(y = 15, linewidth = .5, color = "dimgray", linestyle = "--")
plt.axhline(y = 9, linewidth = .5, color = "dimgray", linestyle = "--")
plt.show()

Pour récupérer la partition, il faut relancer la méthode en changeant le paramètre n_clusters (et en laissant la valeur par défaut pour distance_threshold). On va tester 2 et 3 classes. Pour chaque partition obtenue, nous allons regarder la taille des classes et les valeurs moyennes de chaque variable.

Avec 2 classes

hac2 = AgglomerativeClustering(n_clusters = 2)
hac2.fit(WGI_num)

AgglomerativeClustering()

pandas.DataFrame(hac2.labels_, columns = ["Classe"]).assign(Effectif = 1).groupby("Classe").count()

	Effectif
Classe
0	129
1	73

WGI_num.assign(Classe = hac2.labels_).groupby("Classe").mean()

	Voice and Accountability	Political Stability and Absence of Violence/Terrorism	Government Effectiveness	Regulatory Quality	Rule of Law	Control of Corruption
Classe
0	-0.498394	-0.536321	-0.615752	-0.603835	-0.633716	-0.641430
1	0.854999	0.824537	1.006136	0.996369	1.032080	1.023469

On a ici une classe avec des valeurs plutôt négatives pour chaque variable et une autre avec des valeurs positives.

plt.figure(figsize = (16, 8))
df = pandas.melt(WGI_num.assign(Classe = hac2.labels_), id_vars = 'Classe')
seaborn.boxplot(data = df, y = "variable", x = "value", hue = "Classe")
plt.show()

Avec 3 classes

hac3 = AgglomerativeClustering(n_clusters = 3)
hac3.fit(WGI_num)

AgglomerativeClustering(n_clusters=3)

pandas.DataFrame(hac3.labels_, columns = ["Classe"]).assign(Effectif = 1).groupby("Classe").count()

	Effectif
Classe
0	73
1	66
2	63

WGI_num.assign(Classe = hac3.labels_).groupby("Classe").mean()

	Voice and Accountability	Political Stability and Absence of Violence/Terrorism	Government Effectiveness	Regulatory Quality	Rule of Law	Control of Corruption
Classe
0	0.854999	0.824537	1.006136	0.996369	1.032080	1.023469
1	-0.959242	-1.062068	-1.052668	-0.990906	-1.041426	-1.029118
2	-0.015602	0.014461	-0.158030	-0.198332	-0.206592	-0.235281

On a maintenant une classe avec des valeurs négatives (proches de -1), une autre avec des valeurs positives (proche de 1), et une dernière classe avec des valeurs autour de 0. C'est la classe avec les valeurs positives qui a été coupée en deux.

plt.figure(figsize = (16, 8))
df = pandas.melt(WGI_num.assign(Classe = hac3.labels_), id_vars = 'Classe')
seaborn.boxplot(data = df, y = "variable", x = "value", hue = "Classe")
plt.show()

$k$-means

L'autre méthode de classification classique est la méthode des $k$-means. Elle est très rapide, mais nécessite de déterminer à l'avance le nombre de classes. Nous allons l'appliquer avec 2 et 3 classes afin de voir quelle partition serait la plus intéressante.

2 classes

from sklearn.cluster import KMeans

kmeans2 = KMeans(n_clusters = 2)
kmeans2.fit(WGI_num)

KMeans(n_clusters=2)

On peut avoir ainsi les classes de chaque individus (qui nous servent ici à calculer la taille de chaque classe), ainsi que les centres des classes.

pandas.DataFrame(kmeans2.labels_, columns = ["Classe"]).assign(Effectif = 1).groupby("Classe").count()

	Effectif
Classe
0	82
1	120

kmeans2.cluster_centers_

array([[ 0.83555058,  0.81293581,  0.90273293,  0.87084453,  0.94783308,
         0.92581433],
       [-0.58660919, -0.63045795, -0.66673415, -0.63807482, -0.70108263,
        -0.6995669 ]])

La première classe concerne les pays avec des valeurs positives sur tous les indicateurs, la seconde classe étant ceux avec des valeurs négatives.

WGI_num.assign(Classe = kmeans2.labels_).groupby("Classe").mean()

	Voice and Accountability	Political Stability and Absence of Violence/Terrorism	Government Effectiveness	Regulatory Quality	Rule of Law	Control of Corruption
Classe
0	0.835551	0.812936	0.902733	0.870845	0.947833	0.925814
1	-0.586609	-0.630458	-0.666734	-0.638075	-0.701083	-0.699567

On retrouve ici les deux mêmes classes que la méthode CAH avec 2 classes.

plt.figure(figsize = (16, 8))
df = pandas.melt(WGI_num.assign(Classe = kmeans2.labels_), id_vars = 'Classe')
seaborn.boxplot(data = df, y = "variable", x = "value", hue = "Classe")
plt.show()

En croisant les deux partitions, on voit bien que ce sont les mêmes, à quelques individus prêts.

pandas.crosstab(hac2.labels_, kmeans2.labels_)

col_0	0	1
row_0
0	9	120
1	73	0

3 classes

kmeans3 = KMeans(n_clusters = 3)
kmeans3.fit(WGI_num)

KMeans(n_clusters=3)

La plus grande classe semble être coupée en 2.

pandas.DataFrame(kmeans3.labels_, columns = ["Classe"]).assign(Effectif = 1).groupby("Classe").count()

	Effectif
Classe
0	66
1	49
2	87

On retrouve les 3 mêmes classes que pour la CAH :

• une présente des valeurs très positives ;
• une autre semble être la classe des pays avec des valeurs proches de zéro ;
• une dernère avec des valeurs très faibles.

WGI_num.assign(Classe = kmeans3.labels_).groupby("Classe").mean()

	Voice and Accountability	Political Stability and Absence of Violence/Terrorism	Government Effectiveness	Regulatory Quality	Rule of Law	Control of Corruption
Classe
0	-0.992739	-1.075007	-1.030913	-0.981857	-1.049120	-1.033342
1	1.020859	0.889478	1.344081	1.338705	1.327762	1.272946
2	0.156560	0.211170	-0.043721	-0.068433	-0.025587	-0.025345

plt.figure(figsize = (16, 8))
df = pandas.melt(WGI_num.assign(Classe = kmeans3.labels_), id_vars = 'Classe')
seaborn.boxplot(data = df, y = "variable", x = "value", hue = "Classe")
plt.show()

En croisant les deux partitions à 3 classes (CAH et $k$-means), on voit qu'on retrouve les mêmes.

pandas.crosstab(hac3.labels_, kmeans3.labels_)

col_0	0	1	2
row_0
0	0	49	24
1	64	0	2
2	2	0	61

Choix du nombre de classes avec $k$-means

L'algorithme $k$-means nous permet d'avoir à la fin l'inertie intra-classes, qui représente la disparité des individus à l'intérieur des classes. Plus cette valeur est proche de 0, meilleur est la partition. Malheureusement, la meilleure partition selon ce critère est donc celle avec autant de classes que d'individus (ce qui n'est pas très utile...).

On va donc chercher un point d'inflexion dans la courbe d'évolution de ce critère. Voici comment faire pour avoir ce graphique. Et ici, le point le plus marquant est celui à 2 classes. Ensuite, celui à 3 classes peut montrer aussi une certaine cassure dans l'évolution du critère.

plt.figure(figsize = (16, 6))
inertia = []
for k in range(1, 11):
    kmeans = KMeans(n_clusters = k, init = "random", n_init = 20).fit(WGI_num)
    inertia = inertia + [kmeans.inertia_]
inertia = pandas.DataFrame({"k": range(1, 11), "inertia": inertia})
seaborn.lineplot(data = inertia, x = "k", y = "inertia")
plt.scatter(2, inertia.query('k == 2')["inertia"], c = "red")
plt.scatter(3, inertia.query('k == 3')["inertia"], c = "red")
plt.show()

A faire

Températures mondiales (anomalies)

Nous continuons de travailler ici sur les données de température mondiale HadCRUT4, fournies par Climate Research Unit.

Pour rappel, voici le code pour l'importer avec la nouvelle version.

had = pandas.read_csv("https://crudata.uea.ac.uk/cru/data/temperature/HadCRUT5.0Analysis_gl.txt", header=None)
donnees = pandas.DataFrame(
    [list(map(lambda v: float(v), filter(lambda v: v!= "", h.split(" ")))) for h in had[0][::2]],
    columns = ["Year", "Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec", "Annual"]
).query("Year < 2022")
donnees.tail()

	Year	Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec	Annual
167	2017.0	0.952	1.067	1.065	0.846	0.780	0.658	0.805	0.811	0.729	0.809	0.806	0.815	0.845
168	2018.0	0.711	0.796	0.790	0.822	0.713	0.738	0.733	0.735	0.676	0.869	0.745	0.824	0.763
169	2019.0	0.800	0.844	1.076	0.939	0.778	0.809	0.857	0.858	0.803	0.956	0.937	1.037	0.891
170	2020.0	1.069	1.113	1.094	1.063	0.908	0.825	0.816	0.801	0.867	0.811	1.013	0.693	0.923
171	2021.0	0.701	0.565	0.726	0.760	0.706	0.713	0.792	0.799	0.867	0.907	0.854	0.751	0.762

Vous devez donc réaliser les étapes suivantes, au propre dans un notebook :

• Rechercher une partition intéressante des années
• Représenter cette partition sur le résultat de l'ACP
• Décrire les classes ainsi obtenues

Que peut-on dire globalement ?