Analyse de données sous Python¶
Nous utilisons dans ce TP le module scikit-learn, dédié au Machine Learning. Pour l'utiliser, nous importons les éléments un par un lorsqu'on en a besoin. Et nous allons aussi utiliser d'autres modules, que nous importons ici.
import numpy
import pandas
import matplotlib.pyplot as plt
import seaborn
seaborn.set_style("white")
# permet d'afficher les graphiques dans un notebook
%matplotlib inline
Pour appliquer les deux méthodes, nous allons réaliser une étude de cas, avec les données WGI. Celles-ci donnent certains indicateurs sur la gouvernance des pays (plus d'infos sur cette page). Voici donc le code pour les importer dans un data frame pandas, et leur aperçu. Pour information, plus la valeur est haute, meilleur est le pays sur cet indicateur (et inversement).
# le dropna() permet de supprimer les pays pour lesquels il manque des informations
WGI_complet = pandas.read_csv("WGI.csv").dropna()
WGI_complet
ACP¶
Dans le sous-module decomposition, nous allons importer la fonction PCA(). Le fonctionnement de celle-ci est similaire à toutes les autres dans scikit-learn.
- On créé d'abord un objet (nommé ici
pca) qui va contenir le résultat de l'ACP. Dans la fonctionPCA(), on pourra aussi indiquer les paramètres tels que le nombre de composants à retenir (n_components- ici on garde tout). - Ensuite, on ajuste (fit en anglais) sur des données avec la fonction
fit()de l'objet précédemment créé. Dans cette fonction, nous devons donc passer les données à utiliser (variables uniquement quantitatives).
Si l'on souhaite une ACP normée (ce qui n'est pas le cas ici), il nous faut standardiser les données en amont, avec par exemple la fonction scale() du sous-module preprocessing (importée aussi ici). Bien évidemment, il ne faut passer que des variables quantitatives, sans données manquantes.
from sklearn.decomposition import PCA
from sklearn.preprocessing import scale
# suppression des colonnes non numériques
WGI_num = WGI_complet.drop(columns = ["Country", "Continent"])
pca = PCA()
pca.fit(WGI_num)
Variance expliquée¶
L'objet pca comprend maintenant un certain nombre d'objets et de fonctions. Le premier objet est le tableau des variances expliquées (explained_variance_ - proportionnelles aux valeurs propres) par chaque dimension, et le ratio (proportion) de variance expliquée par dimension (explained_variance_ratio_).
print(pca.explained_variance_)
print(pca.explained_variance_ratio_)
Bien évidemment, il est possible (et préférable) de faire un tableau récapitulatif, avec les variances expliquées, les proportions de variance expliquée simples et cumulées. Voici un petit code permettant de faire ceci.
eig = pandas.DataFrame(
{
"Dimension" : ["Dim" + str(x + 1) for x in range(6)],
"Variance expliquée" : pca.explained_variance_,
"% variance expliquée" : numpy.round(pca.explained_variance_ratio_ * 100),
"% cum. var. expliquée" : numpy.round(numpy.cumsum(pca.explained_variance_ratio_) * 100)
}
)
eig
L'idée est de représenter graphiquement ces proportions de variances expliquées (qu'on passe en pourcentage par convenance).
eig.plot.bar(x = "Dimension", y = "% variance expliquée") # permet un diagramme en barres
plt.text(5, 18, "17%") # ajout de texte
plt.axhline(y = 17, linewidth = .5, color = "dimgray", linestyle = "--") # ligne 17 = 100 / 6 (nb dimensions)
plt.show()
On remarque ici qu'avec seulement deux dimensions suffisent à expliquer 91 % de la variance des données.
Représentation des individus¶
Nous allons maintenant calculer les coordonnées des iris sur les dimensions, avec la fonction transform() de l'objet pca.
WGI_pca = pca.transform(WGI_num)
Afin de manipuler plus facilement l'objet obtenu par la suite, nous allons créer un DataFrame pandas en ne prenant que les deux premières dimensions, ainsi que les espèces.
# Transformation en DataFrame pandas
WGI_pca_df = pandas.DataFrame({
"Dim1" : WGI_pca[:,0],
"Dim2" : WGI_pca[:,1],
"Country" : WGI_complet["Country"],
"Continent" : WGI_complet["Continent"]
})
# Résultat (premières lignes)
WGI_pca_df.head()
Il est maintenant possible de représenter les données sur le premier plan factoriel. Dans ce graphique, il est important de noter le pourcentage de variance expliquée.
WGI_pca_df.plot.scatter("Dim1", "Dim2") # nuage de points
plt.xlabel("Dimension 1 (83%)") # modification du nom de l'axe X
plt.ylabel("Dimension 2 (8 %)") # idem pour axe Y
plt.suptitle("Premier plan factoriel (91%)") # titre général
plt.show()
Dans ce même graphique, on veut repérer les individus aux extrêmes. Ci-dessous, nous affichons donc les noms des pays ayant des valeurs absolues élevées pour chaque composante.
# utilisation de subplots nécessaire car annotation du graphique
fig, ax = plt.subplots()
WGI_pca_df.plot.scatter("Dim1", "Dim2", ax = ax) # l'option ax permet de placer les points et le texte sur le même graphique
# boucle sur chaque pays
for k in WGI_pca_df.iterrows():
# annotation uniquement si valeur absolue sur une de 2 dimensions importantes (valeurs choisies empiriquement)
if (abs(k[1]['Dim1']) > 3.5) | (abs(k[1]['Dim2']) > 1.5):
ax.annotate(k[1]["Country"], (k[1]['Dim1'], k[1]['Dim2']), fontsize = 9)
plt.xlabel("Dimension 1 (83%)")
plt.ylabel("Dimension 2 (8 %)")
plt.suptitle("Premier plan factoriel (91%)")
plt.show()
On pourrait aussi se demander s'il y a un effet continent. Pour cela, nous allons colorer les points en fonction de leur continent donc.
# permet de créer une palette de couleurs, basée sur Color Brewer
palette = plt.get_cmap("Dark2")
# associe une couleur à chaque continent (cf ci-dessous)
couleurs = dict(zip(WGI_pca_df["Continent"].drop_duplicates(), palette(range(6))))
couleurs
Une fois qu'on a associé une couleur à chaque continent, on peut utiliser ce dictionnaire comme ci-dessous
couleurs["Asia"]
position = dict(zip(couleurs.keys(), range(6)))
position
# Affichage des points avec une liste de couleurs
WGI_pca_df.plot.scatter(x = "Dim1", y = "Dim2", c = [couleurs[p] for p in WGI_pca_df["Continent"]])
# boucle pour afficher la légende
for cont, coul in couleurs.items():
plt.scatter(3, position[cont] / 3 + 2.15, c = [coul], s = 20)
plt.text(3.2, position[cont] / 3 + 2, cont)
plt.xlabel("Dimension 1 (83%)")
plt.ylabel("Dimension 2 (8 %)")
plt.suptitle("Premier plan factoriel (91%)")
plt.show()
Il est aussi possible de différencier l'affichage de ce premier plan par continent.
fig = plt.figure()
fig.subplots_adjust(hspace = 0.4, wspace = 0.4)
for cont, coul in couleurs.items():
ax = fig.add_subplot(2, 3, position[cont] + 1)
temp = WGI_pca_df.query('Continent == "' + cont + '"')
temp.plot.scatter(x = "Dim1", y = "Dim2", c = [coul], ax = ax)
ax.set_title(cont)
ax.set_xlim(WGI_pca_df.Dim1.quantile([0, 1])) # permet d'avoir min et max
ax.set_ylim(WGI_pca_df.Dim2.quantile([0, 1])) # idem
ax.set_xlabel("")
ax.set_ylabel("")
plt.show()
Représentation des variables¶
Lorsqu'on réalise une ACP, surtout dans un contexte exploratoire (et non de compression), il est absolument nécessaire de décrire les axesà partir des variables, afin de comprendre ce qu'ils représentent. Pour cela, nous devons réaliser le cercle des corrélations. Malheureusement, il faut le réaliser à la main.
n = WGI_num.shape[0] # nb individus
p = WGI_num.shape[1] # nb variables
eigval = (n-1) / n * pca.explained_variance_ # valeurs propres
sqrt_eigval = numpy.sqrt(eigval) # racine carrée des valeurs propres
corvar = numpy.zeros((p,p)) # matrice vide pour avoir les coordonnées
for k in range(p):
corvar[:,k] = pca.components_[k,:] * sqrt_eigval[k]
# on modifie pour avoir un dataframe
coordvar = pandas.DataFrame({'id': WGI_num.columns, 'COR_1': corvar[:,0], 'COR_2': corvar[:,1]})
coordvar
Une fois ces coordonnées récupérées, on peut donc afficher ce cercle des corrélations comme ci-dessous.
# Création d'une figure vide (avec des axes entre -1 et 1 + le titre)
fig, axes = plt.subplots(figsize = (6,6))
fig.suptitle("Cercle des corrélations")
axes.set_xlim(-1, 1)
axes.set_ylim(-1, 1)
# Ajout des axes
axes.axvline(x = 0, color = 'lightgray', linestyle = '--', linewidth = 1)
axes.axhline(y = 0, color = 'lightgray', linestyle = '--', linewidth = 1)
# Ajout des noms des variables
for j in range(p):
axes.text(coordvar["COR_1"][j],coordvar["COR_2"][j], coordvar["id"][j])
# Ajout du cercle
plt.gca().add_artist(plt.Circle((0,0),1,color='blue',fill=False))
plt.show()
Ici, on va en déduire que l'axe 1 oppose à gauche les pays avec une forte valeur dans tous les indicateurs (ils sont tous à gauche), et à droite les pays à faible valeur pour ces indicateurs.
L'axe 2 semble opposer les pays avec une valeur positive pour les indicateurs sur la statibilité politique et la prise en compte de l'avis des citoyens en haut, avec les pays avec une valeur négative pour ces 2 indicateurs en bas.
Vous pouvez revenir sur le graphique des pays pour mieux comprendre.
Classification¶
Nous allons ici réaliser une classification non supervisée. 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 (avec le critère de Ward), 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.
from sklearn.cluster import AgglomerativeClustering
hac = AgglomerativeClustering(distance_threshold=0, n_clusters=None)
hac.fit(WGI_num)
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. 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.title("CAH (Ward) sur les pays")
# plot the top three levels of the dendrogram
plot_dendrogram(hac)
plt.axhline(y = 25, linewidth = .5, color = "dimgray", linestyle = "--")
plt.axhline(y = 12, linewidth = .5, color = "dimgray", linestyle = "--")
plt.show()
from sklearn.cluster import KMeans
kmeans2 = KMeans(n_clusters = 2)
kmeans2.fit(WGI_num)
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.Series(kmeans2.labels_).value_counts()
kmeans2.cluster_centers_
Cette présentation n'est pas totalement satisfaisante, pour une lecture et analyse des résultats. Nous allons remédier à cela.
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_k2 = WGI_num.assign(classe = kmeans2.labels_)
WGI_k2.groupby("classe").mean()
En combinant ces informations avec l'ACP, on peut représenter ces 2 classes sur le plan factoriel. On voit bien la séparation horizontale des 2 classes.
WGI_pca_k2 = WGI_pca_df.assign(classe = kmeans2.labels_)
WGI_pca_k2.plot.scatter(x = "Dim1", y = "Dim2", c = "classe", cmap = "Accent")
plt.show()
3 classes¶
kmeans3 = KMeans(n_clusters = 3)
kmeans3.fit(WGI_num)
La plus grande classe semble être coupée en 2.
pandas.Series(kmeans3.labels_).value_counts()
La classe 2 présente des valeurs très positives. La classe 0 semble être la classe des pays avec des valeurs proches de zéro. Et la classe 1, celle des pays avec des valeurs très faibles.
WGI_k3 = WGI_num.assign(classe = kmeans3.labels_)
WGI_k3.groupby("classe").mean()
Sur la plan factoriel avec les 3 classes, on retrouve la séparation horizontale des classes, sans grande surprise.
WGI_pca_k3 = WGI_pca_df.assign(classe = kmeans3.labels_)
WGI_pca_k3.plot.scatter(x = "Dim1", y = "Dim2", c = "classe", cmap = "Accent")
plt.show()
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.
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})
inertia.plot.line(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 allons travailler ici sur les données de température mondiale HadCRUT4, fournies par Climate Research Unit. Vous trouverez plus d’informations sur ces données sur ce lien. Nous avons ici plus exactement l'historique des anomalies moyennes mensuelles et annuelles depuis 1850, au niveau mondial, par rapport à la période de référence 1961-1990.
Le code ci-dessous télécharge directement les dernières données disponibles et les met dans un DataFrame dont vous avez un aperçu en dessous.
had = pandas.read_csv("https://crudata.uea.ac.uk/cru/data/temperature/HadCRUT4-gl.dat", 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"]
)
donnees.tail()
Vous devez donc réaliser les étapes suivantes, au propre dans un notebook :
- Décrire rapidement les données
- Calculer les statistiques de base sur chaque mois et sur l'année
- Réprésenter l'évolution des anomalies annuelles sur un graphique
- Réaliser une ACP sur les données mensuelles
- Produire les graphiques nécessaires à l’interprétation
- Identifier des années particulières
- 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 ?