Importation et préparation des données¶

In [ ]:
#dossier de travail
import os
os.chdir("C:/Users/ricco/Desktop/demo")

#chargement des données
import pandas
df = pandas.read_table("data.txt",sep='\t')
df.shape
Out[ ]:
(11069, 8316)
In [ ]:
#partition apprentissage-test
from sklearn.model_selection import train_test_split
dfTrain, dfTest =  train_test_split(df,train_size=1000,random_state=0,stratify=df.target)

#vérif.
#on est dans un exercice de style ici, d'où ces proportions app/test
print(dfTrain.shape)
print(dfTest.shape)

(1000, 8316)
(10069, 8316)
In [ ]:
#distribution des classes
dfTrain.target.value_counts()
Out[ ]:
neg    935
pos     65
Name: target, dtype: int64
In [ ]:
#fréquences relatives
#accuracy du modèle par défaut = 93.5%
dfTrain.target.value_counts(normalize=True)
Out[ ]:
neg    0.935
pos    0.065
Name: target, dtype: float64
In [ ]:
#recoder la cible pour la régression PLS
#recodage 0/1 simple parce qu'évaluation avec courbe ROC
#capacité du modèle à attribuer un score plus élévés aux individus de la classe cible
#pas de problème de seuil d'affectation
import numpy
yTrain = numpy.asarray(dfTrain.target == 'pos').astype(int)

#vérif.
yTrain[:20]
Out[ ]:
array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0])
In [ ]:
#verif. 2
numpy.sum(yTrain)
Out[ ]:
65
In [ ]:
#matrice des X
XTrain = dfTrain.iloc[:,1:]

#vérif.
XTrain.shape
Out[ ]:
(1000, 8315)

Régression PLS pour le scoring¶

Première version - 1 seule composante¶

In [ ]:
#régression PLS - solution initiale : 1 composante pour l'instant
from sklearn.cross_decomposition import PLSRegression
plsInit = PLSRegression(n_components = 1)
plsInit.fit(XTrain,yTrain)

#coefficients de la régression
plsInit.coef_

c:\Users\ricco\anaconda3\envs\env_tutos\lib\site-packages\sklearn\cross_decomposition\_pls.py:503: FutureWarning: The attribute `coef_` will be transposed in version 1.3 to be consistent with other linear models in scikit-learn. Currently, `coef_` has a shape of (n_features, n_targets) and in the future it will have a shape of (n_targets, n_features).
  warnings.warn(
Out[ ]:
array([[ 0.        ],
       [ 0.        ],
       [ 0.        ],
       ...,
       [ 0.        ],
       [-0.00012556],
       [-0.0001751 ]])
In [ ]:
#constante de la régression
plsInit.intercept_
Out[ ]:
array([0.065])
In [ ]:
#valeur de l'AUC en resubstitution
from sklearn import metrics
metrics.roc_auc_score(yTrain, plsInit.predict(XTrain))
Out[ ]:
0.9958535582064993
In [ ]:
#définir un K-fold stratifié
from sklearn.model_selection import StratifiedKFold
skf = StratifiedKFold(n_splits=5,shuffle=True,random_state=0)
skf.get_n_splits(XTrain,yTrain)
Out[ ]:
5
In [ ]:
#AUC-s en validation croisée
from sklearn.model_selection import cross_val_score
auc_res = cross_val_score(plsInit, XTrain, yTrain, scoring='roc_auc',cv=skf)
auc_res
Out[ ]:
array([0.97614151, 0.96503497, 0.94282188, 0.98971617, 0.98354587])
In [ ]:
#AUC synthèse en validation croisée
numpy.mean(auc_res)
Out[ ]:
0.9714520773344304
In [ ]:
#AUC sur l'échantillon test
yTest = numpy.asarray(dfTest.target=='pos').astype(int)
XTest = dfTest.iloc[:,1:]

#
metrics.roc_auc_score(yTest, plsInit.predict(XTest))
Out[ ]:
0.9719458869255511

Recherche du nombre "optimal" de composantes¶

In [ ]:
#nombre max d'axes
max_components = 10

#auc en resubstitution
auc_resub = numpy.zeros(max_components)

#auc en validation croisée
auc_cv = numpy.zeros(max_components)

#itérer
for k in range(max_components):
    #régression PLS sur la totalité des individus
    pls = PLSRegression(n_components=k+1)
    pls.fit(XTrain,yTrain)
    #auc en resubstitution
    auc_resub[k] = metrics.roc_auc_score(yTrain, pls.predict(XTrain))
    #auc en validation croisée
    auc_cv[k] = numpy.mean(cross_val_score(pls, XTrain, yTrain, scoring='roc_auc',cv=skf))
    
#résultats
print(auc_resub)
print(auc_cv)

[0.99585356 0.9993912  0.99990128 1.         1.         1.
 1.         1.         1.         1.        ]
[0.97145208 0.98395722 0.98099548 0.97844508 0.97589469 0.97597696
 0.9747429  0.97359111 0.97202797 0.97112299]
In [ ]:
#graphique - s'intéresser au paramètre de régularisation est toujours important
import matplotlib.pyplot as plt
plt.plot(numpy.arange(1,max_components+1,1),auc_resub,marker='.',label='Resub.')
plt.plot(numpy.arange(1,max_components+1,1),auc_cv,marker='.',label='Cross-val.')
plt.xlabel("# components")
plt.ylabel('ROC AUC')
plt.legend()
plt.show()

Instanciation et évaluation de la solution optimisée¶

In [ ]:
#régression PLS - 2 composantes cette fois-ci
plsOpt = PLSRegression(n_components = 2)
plsOpt.fit(XTrain,yTrain)

#coefficients de la régression
print(plsOpt.coef_,'\n')

#constante
print(plsOpt.intercept_)

[[ 0.        ]
 [ 0.        ]
 [ 0.        ]
 ...
 [ 0.        ]
 [-0.00016245]
 [-0.00016625]] 

[0.065]

c:\Users\ricco\anaconda3\envs\env_tutos\lib\site-packages\sklearn\cross_decomposition\_pls.py:503: FutureWarning: The attribute `coef_` will be transposed in version 1.3 to be consistent with other linear models in scikit-learn. Currently, `coef_` has a shape of (n_features, n_targets) and in the future it will have a shape of (n_targets, n_features).
  warnings.warn(
In [ ]:
#performances en test
metrics.roc_auc_score(yTest, plsOpt.predict(XTest))
Out[ ]:
0.9774529942259497

Classement -- Analyse discriminante PLS¶

On aurait pu utiliser le modèle issu dans la régression PLS pour classer les individus supplémentaires, en comparant le score d'affectation avec la moyenne de la variable recodée. Mais puisque qu'avec la régression PLS on dispose d'un espace de représentation intermédiaire, nous pouvons en tirer profit.

Représentation factorielle des individus¶

In [ ]:
#PLS produit un syst. de représentation factorielle
#avec une propriété de réduction de dimension
#elle est construite en tenant compte de(s) variable(s) cible(s) y
ZTrain = plsOpt.transform(XTrain)
ZTrain.shape
Out[ ]:
(1000, 2)
In [ ]:
#projection des individus dans l'espace factoriel
import matplotlib.pyplot as plt
plt.scatter(ZTrain[:,0], ZTrain[:,1], marker='.')
plt.title("Plan factoriel PLS")
plt.xlabel('Composante 1')
plt.ylabel('Composante 2')
plt.show()
In [ ]:
#graphique - séparabilité des classes
plt.scatter(ZTrain[:,0], ZTrain[:,1], marker='.',c=numpy.array(['red','green'])[yTrain])
plt.title("Plan factoriel PLS + Classes d'appartenance")
plt.xlabel('Composante 1')
plt.ylabel('Composante 2')
plt.show()
In [ ]:
#les facteurs (définis à partir des X) sont orthogonaux
numpy.corrcoef(ZTrain,rowvar=False)
Out[ ]:
array([[ 1.0000000e+00, -4.9009449e-16],
       [-4.9009449e-16,  1.0000000e+00]])

On pourrait aussi construire un cercle des corrélations pour identifier les variables qui déterminent les facteurs, et par conséquent la séparabilité des classes.

Analyse dicriminante à partir de la représentation factorielle¶

Une ADL dans un espace à 8315 explicatives pour 1000 obs. aurait été très instable.

In [ ]:
#partir sur un schéma d'affectation cette fois-ci
#à partir de la représentation factorielle
#analyse discriminante PLS
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
adl = LinearDiscriminantAnalysis()
adl.fit(ZTrain,yTrain)

#coefficients
adl.coef_
Out[ ]:
array([[4.53363764, 1.9962874 ]])
In [ ]:
#constante
adl.intercept_
Out[ ]:
array([-35.50206944])
In [ ]:
#équation explicite de la frontière de séparation des classes
#dixit l'analyse discriminante linéaire
#pente et constante
a = -adl.coef_[0][0]/adl.coef_[0][1]
b = -adl.intercept_[0]/adl.coef_[0][1]

#
print(f"{a} * x1 + {b}")

-2.271034544237145 * x1 + 17.78404725158195
In [ ]:
#représentation de la frontière
#dans l'espace factoriel de la PLS
plt.scatter(ZTrain[:,0], ZTrain[:,1], marker='.',c=numpy.array(['red','green'])[yTrain])
plt.plot([0,30],[a*0+b,a*30+b])
plt.title("Plan factoriel - Frontière de séparation des classes")
plt.xlabel('Composante 1')
plt.ylabel('Composante 2')
plt.show()

Performances en test¶

On pourrait monter un "pipeline" pour enchaîner la PLS + l'analyse discriminante linéaire.

In [ ]:
#coord. des individus en test
#dans le plan factoriel
ZTest = plsOpt.transform(XTest)
ZTest.shape
Out[ ]:
(10069, 2)
In [ ]:
#représentation des points supplémentaires (de l'éch. test) dans la plan
plt.scatter(ZTest[:,0], ZTest[:,1], marker='.',s=10)
plt.title("Plan factoriel - Individus de l'ech. test")
plt.xlabel('Composante 1')
plt.ylabel('Composante 2')
plt.show()
In [ ]:
#prédiction +  évaluation = accuracy ici
adl.score(ZTest,yTest)
Out[ ]:
0.9637501241434104
In [ ]:
#détail - matrice de confusion
#visibilité sur la structure de l'erreur
metrics.ConfusionMatrixDisplay.from_predictions(yTest, adl.predict(ZTest))
Out[ ]:
<sklearn.metrics._plot.confusion_matrix.ConfusionMatrixDisplay at 0x16c73b2b010>
In [ ]:
#vérif. de l'accuracy
(9365+339)/10069
Out[ ]:
0.9637501241434104
In [ ]:
#dans la plan factoriel
#classements corrects et erronés
plt.scatter(ZTest[:,0], ZTest[:,1], marker='.',s=10,c=numpy.array(['red','green'])[numpy.array(dfTest.target=='pos').astype(int)])
plt.plot([0,20],[a*0+b,a*20+b])
plt.title("Plan factoriel - Frontière de séparation des classes")
plt.xlabel('Composante 1')
plt.ylabel('Composante 2')
plt.show()

Comparaison avec la régression pénalisée (Ridge)¶

Définie sur l'espace initial (les variables originelles) de représentation.

In [ ]:
#par comparaison une régression logistique
#par défaut dans sklearn = pénalisée Ridge
#https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html
from sklearn.linear_model import LogisticRegression
clf = LogisticRegression()
clf.fit(XTrain,yTrain)

#perfs
clf.score(XTest,yTest)
Out[ ]:
0.9433906048266958