Analyses en composantes principales et une de ses alternatives, le tSNE
en langage python
L'essentiel de cette page
L'analyse en composantes principales (ACP ou PCA en anglais) permet de réduire le nombre de dimensions d'un jeu de données décrit par un grand nombre de variables. Cela permet une visualisation simplifiée et une accélération des calculs.
Le module scikit-learn permet de réaliser une telle réduction mais présente certaines faiblesses (écart-type utilisé non-approprié)
0- Dans cette page, nous partirons pour l'exemple du jeu de données Iris
import seaborn as sns # nécessite le module seaborn
import matplotlib.pyplot as plt
iris = sns.load_dataset('iris') # Charger le jeu de données iris
# Les 4 colonnes de mesures qui décrivent différents iris
X = iris.iloc[:,0:4]
colors = iris.species.astype('category')
# Les couleurs déduites de chaque espèce d'iris
y = colors.cat.codes
print(X) ; print(y)
from pandas.plotting import scatter_matrix
scatter_matrix(X,c=y) ; plt.show()
1- Réaliser une ACP avec scikit-learn (sklearn)
1.1- Centrer-réduire les valeurs (optionnel)
1) Méthode de base sklearn
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
scaler.fit(X)
X=scaler.transform(X)
Cette méthode prête à confusion : il semblerait que l'écart-type utilisé est celui d'une population et non d'un échantillon.
2) Autre méthode sklearn
from sklearn.preprocessing import scale
X = scale(X)
3) Autre méthode, pour les dataframe (pandas) - on obtient ainsi un X au format dataframe.
temp = X.sub(X.mean()) # Soustraire la moyenne de chaque colonne à chaque valeur : centrer les valeurs
X= temp.div(temp.std()) # Divisé les valeurs de la dataframe par l'écart-type de chaque colonne
1.2- Réaliser l'ACP
from sklearn.decomposition import PCA
mypca = PCA(n_components=3) # On paramètre ici pour ne garder que 3 composantes
# Modèle d'ACP
mypca.fit(X)
# Pourcentage de la variance expliquée par chacune des composantes sélectionnées.
print(mypca.singular_values_) # Valeurs de variance
print(mypca.explained_variance_ratio_) # Pourcentages
# Axes principaux dans l'espace des caractéristiques, représentant les directions de la variance maximale dans les données. Les composantes sont triées par variance expliquée.
print(mypca.components_) #
# Résultats de l'ACP
data_sortie= mypca.fit_transform(X)
# Bruit estimé lié à la covariance
print(mypca.noise_variance_)
Pour comprendre ces paramètres : lien externe.
1.3- Contrôler l'ACP par un diagramme de pareto
Un diagramme de pareto permet de s'assurer que les composantes retenues retiendront bien la majorité de la variabilité.
Un code rustique qui marche :
import numpy as np
y = list(mypca.explained_variance_ratio_)
x = range(len(y))
ycum = np.cumsum(y)
plt.bar(x,y)
plt.plot(x,ycum,"-r")
plt.show()
Remarque importante : pour qu'un pareto soit pertinent, il ne faut pas fixer le nombre de composante retenues sinon on ignore volontairement une partie de la variance perdu. Ici, si on avait eu plus de 4 variables, le choix de 3 composantes nous aurait conduit à négliger la variance des composantes au-delà de 3.
Un code plus élaboré et plus simple à utiliser :
pareto(mypca.explained_variance_ratio_)
Pour cela : copier-coller le code de la fonction ci-dessous
Ici je ne peux avoir que 3 composantes (car j'ai 4 variables à la base). Mes deux premières composantes absorbent bien la majorité de la variabilité.
Code de la fonction pareto() à copier-coller (cliquer sur le titre)
def pareto(data) :
from matplotlib.ticker import PercentFormatter
import numpy as np
y = list(data)
x = range(len(data))
ycum = np.cumsum(y)/sum(y)*100
fig, ax = plt.subplots()
ax.bar(x,y,color="yellow")
ax2 = ax.twinx()
ax2.plot(x,ycum, color="C1", marker="D", ms=7)
ax2.axhline(y=80,color="r")
ax2.yaxis.set_major_formatter(PercentFormatter())
ax.tick_params(axis="y", colors="C0")
ax2.tick_params(axis="y", colors="C1")
plt.ylim(0,110)
plt.show()
2- Visualiser les résultats de l'ACP
biplot - graphique qui représente la projection des individus sur deux composantes sélectionnées et la corrélation des variables avec ces composantes.
biplot(score=data_sortie[:,0:2],
coeff=np.transpose(mypca.components_[0:2, :]),cat=y,
cmap="viridis")
plt.show()
Pour avoir accès à la fonction biplot(), copier-coller le code ci-dessous à la fin de cette partie.
On peut obtenir directement un affichage avec le bon nom des variables si X est au format pandas :
biplot(mypca,x=X,cat=y,components=[0,1])
plt.show()
components permet ici d'indiquer quelles dimensions/composantes l'on souhaite afficher...
On peut aussi paramétrer la palette de couleurs et indiquer des labels pour les variables ou les individus :
cmap
coeff_labels
score_labels
Lorsqu'on travaille avec jeu de données trop grandes, on peut ne tracer que les nappes convexes des catégories afin de voir si elles ne se superposent pas
Paramètre : bigdata, au-delà de 100 valeurs, on passe en nappes convexes
Remarque : ici X est au format dataframe de pandas (et non numpy.array) : il a été centré-réduit selon la méthode 3.
biplot(score=data_sortie[:,0:2],
coeff=np.transpose(mypca.components_[0:2, :]),cat=y, bigdata=100,coeff_labels = list(X.columns),density=False)
plt.show()
On voit ainsi ici que 2 catégories se superposent.
Lorsque le paramètre BigData fait basculer en un affichage d'un grand jeu de valeurs, on peut aussi demander une graphique de densité
biplot(score=data_sortie[:,0:2],cat=y, bigdata=100,density=True)
plt.show()
On voit ainsi les 3 catégories à travers les 3 nappes convexes et l'on perçoit bien où se trouvent la majorité des points s'il y en avait trop pour tous les voir.
Code de la fonction biplot() à copier-coller (cliquer sur le titre)
#######################################
# biplot
# version 12/11/2021
#######################################
import matplotlib.pyplot as plt
from scipy.spatial import ConvexHull
import numpy as np
import pandas as pd
import matplotlib as mpl
import matplotlib.cm as cm
import seaborn as sns
from sklearn.decomposition import PCA
def biplot(pca=[],x=None,y=None,components=[0,1],score=None,coeff=None,coeff_labels=None,score_labels=None,circle='T',bigdata=1000,cat=None,cmap="viridis",density=True):
if isinstance(pca,PCA)==True :
coeff = np.transpose(pca.components_[components, :])
score= pca.fit_transform(x)[:,components]
if isinstance(x,pd.DataFrame)==True :
coeff_labels = list(x.columns)
if score is not None : x = score
if x.shape[1]>1 :
xs = x[:,0]
ys = x[:,1]
else :
xs = x
ys = y
if (len(xs) != len(ys)) : print("Warning ! x et y n'ont pas la même taille !")
scalex = 1.0/(xs.max() - xs.min())
scaley = 1.0/(ys.max() - ys.min())
#x_c = xs * scalex
#y_c = ys * scaley
temp = (xs - xs.min())
x_c = temp / temp.max() * 2 - 1
temp = (ys - ys.min())
y_c = temp / temp.max() * 2 - 1
data = pd.DataFrame({"x_c":x_c,"y_c":y_c})
print("Attention : pour des facilités d'affichage, les données sont centrées-réduites")
if cat is None : cat = [0]*len(xs)
elif len(pd.Series(cat)) == 1 : cat = list(pd.Series(cat))*len(xs)
elif len(pd.Series(cat)) != len(xs) : print("Warning ! Nombre anormal de catégories !")
cat = pd.Series(cat).astype("category")
fig = plt.figure(figsize=(6,6),facecolor='w')
ax = fig.add_subplot(111)
# Affichage des points
if (len(xs) < bigdata) :
ax.scatter(x_c,y_c, c = cat.cat.codes,cmap=cmap)
if density==True : print("Warning ! Le mode density actif n'apparait que si BigData est paramétré.")
# Affichage des nappes convexes (BigData)
else :
#color
norm = mpl.colors.Normalize(vmin=0, vmax=(len(np.unique(cat.cat.codes)))) #-(len(np.unique(c)))
cmap = cmap
m = cm.ScalarMappable(norm=norm, cmap=cmap)
if density==True :
sns.set_style("white")
sns.kdeplot(x="x_c",y="y_c",data=data)
if len(np.unique(cat)) <= 1 :
sns.kdeplot(x="x_c",y="y_c",data=data, cmap="Blues", shade=True, thresh= 0)
else :
for i in np.unique(cat) :
color_temp = m.to_rgba(i)
sns.kdeplot(x="x_c",y="y_c",data=data[cat==i], color=color_temp,
shade=True, thresh=0.25, alpha=0.25)
for cat_temp in cat.cat.codes.unique() :
x_c_temp = [x_c[i] for i in range(len(x_c)) if (cat.cat.codes[i] == cat_temp)]
y_c_temp = [y_c[i] for i in range(len(y_c)) if (cat.cat.codes[i] == cat_temp)]
points = [ [ None ] * len(x_c_temp) ] * 2
points = np.array(points)
points = points.reshape(len(x_c_temp),2)
points[:,0] = x_c_temp
points[:,1] = y_c_temp
hull = ConvexHull(points)
temp = 0
for simplex in hull.simplices:
color_temp = m.to_rgba(cat_temp)
plt.plot(points[simplex, 0], points[simplex, 1],color=color_temp)#, linestyle='dashed')#linewidth=2,color=cat)
if (temp == 0) :
plt.xlim(-1,1)
plt.ylim(-1,1)
temp = temp+1
if coeff is not None :
if (circle == 'T') :
x_circle = np.linspace(-1, 1, 100)
y_circle = np.linspace(-1, 1, 100)
X, Y = np.meshgrid(x_circle,y_circle)
F = X**2 + Y**2 - 1.0
#fig, ax = plt.subplots()
plt.contour(X,Y,F,[0])
n = coeff.shape[0]
for i in range(n):
plt.arrow(0, 0, coeff[i,0], coeff[i,1],color = 'r',alpha = 0.5,
head_width=0.05, head_length=0.05)
if coeff_labels is None:
plt.text(coeff[i,0]* 1.15, coeff[i,1] * 1.15, "Var"+str(i+1), color = 'g', ha = 'center', va = 'center')
else:
plt.text(coeff[i,0]* 1.15, coeff[i,1] * 1.15, coeff_labels[i], color = 'g', ha = 'center', va = 'center')
if score_labels is not None :
for i in range(len(score_labels)) :
temp_x = xs[i] * scalex
temp_y = ys[i] * scaley
plt.text(temp_x,temp_y,list(score_labels)[i])
plt.xlim(-1.2,1.2)
plt.ylim(-1.2,1.2)
plt.xlabel("PC{}".format(1))
plt.ylabel("PC{}".format(2))
plt.grid(linestyle='--')
plt.show()
3- Une alternative l'ACP, le plongement par tSNE
Le principe du tSNE est de réduire les dimensions en considérant que tous les points sont reliés par des petits ressorts. De façon stochastique, et par itération, le tSNE va essayer de proposer un modèle de réduction des dimensions correspondant à un modèle d'énergie minimisée.
Ce site n'explique pas le fonctionnement du tSNE mais comment en faire un en python avec sci-kit learn (module sklearn).
1) Mettre en place le tSNE
from sklearn.manifold import TSNE
# Instanciation
tsne = TSNE(n_components=2) # On peut choisir de réduire à 2 dimensions à plus
# Itérations
X_embedded = tsne.fit_transform(X)
# Résultat
print(X_embedded.shape)
2) Visualiser les résultats
# 2 composantes
from matplotlib import pyplot as plt
plt.scatter(X_embedded[:,0],X_embedded[:,1],c=y, cmap='Set1')
plt.show()
# Si on refait le tSNE avec 3 composantes au lieu de 2 : on voit que le passage de 2 à 3 composantes change la sortie du tSNE
from matplotlib import pyplot as plt
plt.subplot(121)
plt.scatter(X_embedded[:,0],X_embedded[:,1],c=y, cmap='Set1') # c2 = f(c1)
plt.subplot(122)
plt.scatter(X_embedded[:,0],X_embedded[:,2],c=y, cmap='Set1') # c3 = f(c2)
plt.show()
On peut aussi utiliser la fonction biplot() ci-dessus (partie ACP) pour visualiser rapidement les catégories représentées par des nappes convexes lorsqu'on dispose de trop de données :
biplot(x=X_embedded[:,0],y=X_embedded[:,1],c=y,bigdata=50)
plt.show()
Plutôt que les points sont représentés ici les contours des nuages des points
4- Adapter l'affichage des résultats de l'ACP ou du tSNE
4.1- Afficher les labels des points dans un graphique simple (non biplot)
Pour afficher les labels des points, le plus simple est de faire une boucle avec la fonction plt.text de matplotlib.pyplot :
from matplotlib import pyplot as plt
plt.scatter(data_sortie[:,0],data_sortie[:,1],c=col, cmap='Set1')
plt.show()
Sera à remplacer par :
Remarque : cette figure ne représente pas un résultat d'ACP mais un résultat d'une méthode alternative, le tSNE.
# listes des labels des individus
labels = list(range(1,80,1))
plt.scatter(data_sortie[:,0],data_sortie[:,1],c=col, cmap='Set1')
plt.grid(color='white', linestyle='solid')
plt.title("Scatter Plot", size=20)
for i in range(len(labels)) :
plt.text(data_sortie[i,0]data_sortie[i,1],labels[i])
plt.show()
On peut aussi adapter le code de la façon suivante :
labels = list(range(1,80,1))
fig, ax = plt.subplots()
plt.scatter(data_sortie[:,0],data_sortie[:,1],c=col, cmap='Set1')
ax.grid(color='white', linestyle='solid')
ax.set_title("Scatter Plot (with tooltips!)", size=20)
for label, x, y in zip(labels, data_sortie[:,0], data_sortie[:,1]):
plt.annotate(
label,
xy=(x, y), xytext=(-20, 20),
textcoords='offset points', ha='right', va='bottom',
bbox=dict(boxstyle='round,pad=0.5', fc='yellow', alpha=0.5),
arrowprops=dict(arrowstyle = '->', connectionstyle='arc3,rad=0'))
plt.show()
Remarque :
Il existe une alternative où le graphique sera construit avec le module mpld3 et affiché dans l'explorateur internet : lien externe. Il suffira alors de passer la souris au-dessus de chaque point pour avoir son nom.
4.2- Afficher 3 composantes avec un graphique 3D
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.scatter(data_sortie[:,0],data_sortie[:,1],
data_sortie[:,2],
zdir='z',s=40,depthshade=True,c=y)
plt.show()
Pour aller plus loin sur la réalisation de graphiques 3D : voir la page des graphiques 3D.
Remarque
Une fois la réduction de dimension faite, il peut être important de réaliser des classification : on peut alors faire appelle à différentes méthodes telles le k-means ou la classification hiérarchique ascendante (CAH).
Réaliser des réductions dimensionnelles sous R telles l'ACP, le cmdscale, la projection de Fischer...
