Analyse de données sous Python - correction¶
import numpy
import pandas
import matplotlib.pyplot as plt
# permet d'afficher les graphiques dans un notebook
%matplotlib inline
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 ?
donnees.drop(columns = ("Year")).describe().round(2)
donnees.drop(columns = ("Year")).boxplot()
Evolution des anomalies¶
On peut aussi réprésenter l'évolution des anomalies au cours du temps, pour l'ensemble des mois. Tous semblent avoir le même comportement.
donnees.plot(x = "Year", figsize = (12,4))
plt.show()
En se concentrant sur la moyenne annuelle des anomalies par rapport à la période de référence (1961-1990 - représentée ci-dessous), on remarque une augmentation des anomalies à partir des années 80.
fig, ax = plt.subplots()
donnees.plot(x = "Year", y = "Annual", figsize = (12,4), ax = ax)
rect = plt.Rectangle((1961, -1), 30, 2,
facecolor = "black", alpha = 0.1)
ax.add_patch(rect)
ax.axhline(y = 0, linewidth = .5, color = "dimgray", linestyle = "--")
ax.text(1962, -0.55, "Période de référence")
plt.xlabel("")
plt.show()
ACP sur les données mensuelles¶
On peut appliquer l'ACP sur les données standardisées (avec la fonction scale()).
from sklearn.decomposition import PCA
from sklearn.preprocessing import scale
donnees_mens = donnees.drop(columns = ["Year", "Annual"])
pca = PCA()
pca.fit(scale(donnees_mens))
Le premier axe explique à lui seul 88 % de l'information.
eig = pandas.DataFrame(
{
"Dimension" : ["Dim" + str(x + 1) for x in range(12)],
"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
eig.plot.bar(x = "Dimension", y = "% variance expliquée") # permet un diagramme en barres
plt.show()
donnees_pca = pca.transform(donnees_mens)
donnees_pca_df = pandas.DataFrame({
"Dim1" : donnees_pca[:,0],
"Dim2" : donnees_pca[:,1],
"Year" : donnees["Year"]
})
donnees_pca_df.head()
En représentant les années sur le plan factoriel, et en faisant apparaître les points extrêmes, on remarque que les dernières années (2015, 2016, 2017 et 2019) se détache des autres. Ce sont les 4 années les plus chaudes depuis 1850.
# utilisation de subplots nécessaire car annotation du graphique
fig, ax = plt.subplots(figsize = (16,4))
donnees_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 donnees_pca_df.iterrows():
# annotation uniquement si valeur absolue sur une de 2 dimensions importantes (valeurs choisies empiriquement)
if (abs(k[1]['Dim1']) > 8) | (abs(k[1]['Dim2']) > 1.5):
ax.annotate(int(k[1]["Year"]), (k[1]['Dim1'], k[1]['Dim2']), fontsize = 9)
plt.xlabel("Dimension 1 (88%)")
plt.ylabel("Dimension 2 ( 4%)")
plt.suptitle("Premier plan factoriel (92%)")
plt.show()
En représentant les mois, on remarque qu'ils sont tous situés à droite du plan factoriel. Cela confirme bien que les dernières années sont celles avec des valeurs les plus élevées.
n = donnees_mens.shape[0] # nb individus
p = donnees_mens.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': donnees_mens.columns, 'COR_1': corvar[:,0], 'COR_2': corvar[:,1]})
fig, axes = plt.subplots(figsize = (4,4))
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()
Classification des années¶
On peut chercher maintenant une partition des années sur la base de ces anomalies.
from sklearn.cluster import AgglomerativeClustering
hac = AgglomerativeClustering(distance_threshold=0, n_clusters=None)
hac.fit(scale(donnees_mens))
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)
plt.title("CAH (Ward) sur les pays")
# plot the top three levels of the dendrogram
plot_dendrogram(hac)
plt.axhline(y = 37, linewidth = .5, color = "dimgray", linestyle = "--")
plt.axhline(y = 21, linewidth = .5, color = "dimgray", linestyle = "--")
plt.show()
Sur le graphique ci-dessous, on remarque qu'on peut choisir soit 2 classes, soit 3 classes. Ici, nous allons choisir 3 classes.
from sklearn.cluster import KMeans
kmeans3 = KMeans(n_clusters = 3)
kmeans3.fit(scale(donnees_mens))
Ces trois classes se répartissent comme suit :
- Classe 0 : années avec des valeurs globalement proche de 0
- Classe 1 : années chaudes (valeurs élévés)
- Classe 2 : années froides (valeurs faibles)
donnees_k3 = donnees_mens.assign(classe = kmeans3.labels_)
donnees_k3.groupby("classe").mean()
col = { 0: "gray", 1: "darkred", 2: "steelblue"}
donnees_pca_k3 = donnees_pca_df.assign(classe = kmeans3.labels_)
donnees_pca_k3.plot.scatter(x = "Dim1", y = "Dim2", c = [col[k] for k in kmeans3.labels_])
plt.show()
Conclusion¶
Il semble ainsi clair qu'il y a un échauffement global de la planète, en particulier sur les dernières années comme on peut le voir sur le graphique ci-dessous.
donnees_trans = donnees.drop(columns =["Year", "Annual"]).transpose()
fig, ax = plt.subplots(figsize = (12, 8))
im = ax.imshow(donnees_trans, cmap = "coolwarm", aspect = 4)
locs, labels = plt.xticks()
ax.set_xticklabels([int(l) for l in locs + 1850])
month = list(donnees.columns[1:13])
ax.set_yticks(numpy.arange(len(month)))
ax.set_yticklabels(month)
plt.setp(ax.get_xticklabels(), rotation = 45, ha = "right", rotation_mode = "anchor")
fig.colorbar(im, aspect = 1)
fig.tight_layout()
plt.show()
