2A.data - Classification, régression, anomalies - correction

Links: notebook, html, PDF, python, slides, GitHub

Le jeu de données Wine Quality Data Set contient 5000 vins décrits par leurs caractéristiques chimiques et évalués par un expert. Peut-on s’approcher de l’expert à l’aide d’un modèle de machine learning.

%matplotlib inline
import matplotlib.pyplot as plt
from jyquickhelper import add_notebook_menu
add_notebook_menu()

Les données

On peut les récupérer sur github…data_2a.

from ensae_teaching_cs.data import wines_quality
from pandas import read_csv
df = read_csv(wines_quality(local=True, filename=True))
df.head()
fixed_acidity volatile_acidity citric_acid residual_sugar chlorides free_sulfur_dioxide total_sulfur_dioxide density pH sulphates alcohol quality color
0 7.4 0.70 0.00 1.9 0.076 11.0 34.0 0.9978 3.51 0.56 9.4 5 red
1 7.8 0.88 0.00 2.6 0.098 25.0 67.0 0.9968 3.20 0.68 9.8 5 red
2 7.8 0.76 0.04 2.3 0.092 15.0 54.0 0.9970 3.26 0.65 9.8 5 red
3 11.2 0.28 0.56 1.9 0.075 17.0 60.0 0.9980 3.16 0.58 9.8 6 red
4 7.4 0.70 0.00 1.9 0.076 11.0 34.0 0.9978 3.51 0.56 9.4 5 red
ax = df['quality'].hist(bins=16)
ax.set_title("Distribution des notes");
../_images/td2a_correction_cl_reg_anomaly_5_0.png

Il y a peu de très mauvais ou très bons vins. On découpe en apprentissage / test ce qui va nécessairement rendre leur prédiction complexe : un modèle reproduit en quelque sorte ce qu’il voit.

from sklearn.model_selection import train_test_split
X = df.drop("quality", axis=1)
y = df["quality"]
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.5)
X_train.shape, X_test.shape
((3248, 12), (3249, 12))
from pandas import DataFrame

def distribution(y_train, y_test):
    df_train = DataFrame(dict(color=y_train))
    df_test = DataFrame(dict(color=y_test))
    df_train['ctrain'] = 1
    df_test['ctest'] = 1
    h_train = df_train.groupby('color').count()
    h_test  = df_test.groupby('color').count()
    merge = h_train.join(h_test, how='outer')
    merge["ratio"] = merge.ctest / merge.ctrain
    return merge

distribution(y_train, y_test)
ctrain ctest ratio
color
3 17 13 0.764706
4 97 119 1.226804
5 1075 1063 0.988837
6 1399 1437 1.027162
7 563 516 0.916519
8 95 98 1.031579
9 2 3 1.500000
ax = y_train.hist(bins=24, label="train", align="right")
y_test.hist(bins=24, label="test", ax=ax, align="left")
ax.set_title("Distribution des notes")
ax.legend();
../_images/td2a_correction_cl_reg_anomaly_9_0.png

Avec un peu de chance, les notes extrêmes sont présentes dans les bases d’apprentissages et de tests mais une note seule a peu d’influence sur un modèle. Pour s’assurer une meilleur répartition train / test, on peut s’assurer que chaque note est bien répartie de chaque côté. On se sert du paramètre stratify.

X_train, X_test, y_train, y_test = train_test_split(X, y, stratify=y, test_size=0.5)
X_train.shape, X_test.shape
ax = y_train.hist(bins=24, label="train", align="right")
y_test.hist(bins=24, label="test", ax=ax, align="left")
ax.set_title("Distribution des notes - statifiée")
ax.legend();
../_images/td2a_correction_cl_reg_anomaly_11_0.png
distribution(y_train, y_test)
ctrain ctest ratio
color
3 15 15 1.000000
4 108 108 1.000000
5 1069 1069 1.000000
6 1418 1418 1.000000
7 539 540 1.001855
8 96 97 1.010417
9 3 2 0.666667

La répartition des notes selon apprentissage / test est plus uniforme.

Premier modèle

from sklearn.linear_model import LogisticRegression
logreg = LogisticRegression()
try:
    logreg.fit(X_train, y_train)
except Exception as e:
    print(e)
could not convert string to float: 'red'

Une colonne n’est pas numérique. On utilise un OneHotEncoder.

from sklearn.preprocessing import OneHotEncoder
one = OneHotEncoder()
one.fit(X_train[['color']])
tr = one.transform(X_test[["color"]])
tr
<3249x2 sparse matrix of type '<class 'numpy.float64'>'
    with 3249 stored elements in Compressed Sparse Row format>

La matrice est sparse ou creuse.

tr.todense()[:5]
matrix([[0., 1.],
        [0., 1.],
        [0., 1.],
        [0., 1.],
        [1., 0.]])

Ensuite il faut fusionner ces deux colonnes avec les données ou une seule puisqu’elles sont corrélées. Ou alors on écrit un pipeline…

from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer

numeric_features = [c for c in X_train if c != 'color']

pipe = Pipeline([
    ("prep", ColumnTransformer([
        ("color", Pipeline([
            ('one', OneHotEncoder()),
            ('select', ColumnTransformer([('sel1', 'passthrough', [0])]))
        ]), ['color']),
        ("others", "passthrough", numeric_features)
    ])),
])

pipe.fit(X_train)
pipe.transform(X_test)[:2]
array([[0.000e+00, 7.200e+00, 2.000e-01, 3.400e-01, 2.700e+00, 3.200e-02,
        4.900e+01, 1.510e+02, 9.900e-01, 3.160e+00, 3.900e-01, 1.270e+01],
       [0.000e+00, 6.700e+00, 2.600e-01, 2.600e-01, 4.100e+00, 7.300e-02,
        3.600e+01, 2.020e+02, 9.956e-01, 3.300e+00, 6.700e-01, 9.500e+00]])
from jyquickhelper import RenderJsDot
from mlinsights.plotting import pipeline2dot
dot = pipeline2dot(pipe, X_train)
RenderJsDot(dot)

Il reste quelques bugs. On ajoute un classifieur.

pipe = Pipeline([
    ("prep", ColumnTransformer([
        ("color", Pipeline([
            ('one', OneHotEncoder()),
            ('select', ColumnTransformer([('sel1', 'passthrough', [0])]))
        ]), ['color']),
        ("others", "passthrough", numeric_features)
    ])),
    ("lr", LogisticRegression(max_iter=1000)),
])

pipe.fit(X_train, y_train)
C:xavierdupre__home_github_forkscikit-learnsklearnlinear_modellogistic.py:932: ConvergenceWarning: lbfgs failed to converge (status=1): b'STOP: TOTAL NO. of ITERATIONS REACHED LIMIT'. Increase the number of iterations.
  n_iter_i = _check_optimize_result(solver, opt_res, max_iter)
Pipeline(memory=None,
         steps=[('prep',
                 ColumnTransformer(n_jobs=None, remainder='drop',
                                   sparse_threshold=0.3,
                                   transformer_weights=None,
                                   transformers=[('color',
                                                  Pipeline(memory=None,
                                                           steps=[('one',
                                                                   OneHotEncoder(categories='auto',
                                                                                 drop=None,
                                                                                 dtype=<class 'numpy.float64'>,
                                                                                 handle_unknown='error',
                                                                                 sparse=True)),
                                                                  ('select',
                                                                   ColumnTransformer(n_jobs=None,
                                                                                     remainder='drop...
                                                   'free_sulfur_dioxide',
                                                   'total_sulfur_dioxide',
                                                   'density', 'pH', 'sulphates',
                                                   'alcohol'])],
                                   verbose=False)),
                ('lr',
                 LogisticRegression(C=1.0, class_weight=None, dual=False,
                                    fit_intercept=True, intercept_scaling=1,
                                    l1_ratio=None, max_iter=1000,
                                    multi_class='auto', n_jobs=None,
                                    penalty='l2', random_state=None,
                                    solver='lbfgs', tol=0.0001, verbose=0,
                                    warm_start=False))],
         verbose=False)
from sklearn.metrics import classification_report
print(classification_report(y_test, pipe.predict(X_test)))
              precision    recall  f1-score   support
           3       0.00      0.00      0.00        15
           4       0.00      0.00      0.00       108
           5       0.58      0.59      0.59      1069
           6       0.50      0.70      0.58      1418
           7       0.42      0.14      0.21       540
           8       0.00      0.00      0.00        97
           9       0.00      0.00      0.00         2
    accuracy                           0.52      3249
   macro avg       0.21      0.20      0.20      3249
weighted avg       0.48      0.52      0.48      3249
C:xavierdupre__home_github_forkscikit-learnsklearnmetricsclassification.py:1428: UndefinedMetricWarning: Precision and F-score are ill-defined and being set to 0.0 in labels with no predicted samples.
  'precision', 'predicted', average, warn_for)

Pas extraordinaire.

from sklearn.ensemble import RandomForestClassifier

pipe = Pipeline([
    ("prep", ColumnTransformer([
        ("color", Pipeline([
            ('one', OneHotEncoder()),
            ('select', ColumnTransformer([('sel1', 'passthrough', [0])]))
        ]), ['color']),
        ("others", "passthrough", numeric_features)
    ])),
    ("lr", RandomForestClassifier()),
])

pipe.fit(X_train, y_train)
print(classification_report(y_test, pipe.predict(X_test)))
              precision    recall  f1-score   support
           3       0.00      0.00      0.00        15
           4       0.45      0.09      0.15       108
           5       0.69      0.70      0.69      1069
           6       0.62      0.74      0.68      1418
           7       0.64      0.51      0.56       540
           8       0.84      0.33      0.47        97
           9       0.00      0.00      0.00         2
    accuracy                           0.65      3249
   macro avg       0.46      0.34      0.37      3249
weighted avg       0.64      0.65      0.64      3249
C:xavierdupre__home_github_forkscikit-learnsklearnmetricsclassification.py:1428: UndefinedMetricWarning: Precision and F-score are ill-defined and being set to 0.0 in labels with no predicted samples.
  'precision', 'predicted', average, warn_for)

Beaucoup mieux.

Courbe ROC pour chaque classe

from sklearn.metrics import roc_curve, auc

labels = pipe.steps[1][1].classes_
y_score = pipe.predict_proba(X_test)

fpr = dict()
tpr = dict()
roc_auc = dict()
for i, cl in enumerate(labels):
    fpr[cl], tpr[cl], _ = roc_curve(y_test == cl, y_score[:, i])
    roc_auc[cl] = auc(fpr[cl], tpr[cl])
fig, ax = plt.subplots(1, 1, figsize=(8,4))
for k in roc_auc:
    ax.plot(fpr[k], tpr[k], label="c%d = %1.2f" % (k, roc_auc[k]))
ax.legend();
../_images/td2a_correction_cl_reg_anomaly_31_0.png

Ces chiffres peuvent paraître élevés mais ce n’est pas formidable quand même.

from sklearn.metrics import confusion_matrix
confusion_matrix(y_test, pipe.predict(X_test), labels=labels)
array([[   0,    0,    6,    8,    1,    0,    0],
       [   1,   10,   62,   34,    1,    0,    0],
       [   0,   11,  746,  300,   12,    0,    0],
       [   0,    1,  260, 1045,  110,    2,    0],
       [   0,    0,   12,  250,  274,    4,    0],
       [   0,    0,    0,   34,   31,   32,    0],
       [   0,    0,    0,    1,    1,    0,    0]], dtype=int64)

Ce n’est pas très joli…

def confusion_matrix_df(y_test, y_true):
    conf = confusion_matrix(y_test, y_true)
    labels = list(sorted(set(y_test)))
    df = DataFrame(conf, columns=labels)
    df.set_index(labels)
    return df

confusion_matrix_df(y_test, pipe.predict(X_test))
3 4 5 6 7 8 9
0 0 0 6 8 1 0 0
1 1 10 62 34 1 0 0
2 0 11 746 300 12 0 0
3 0 1 260 1045 110 2 0
4 0 0 12 250 274 4 0
5 0 0 0 34 31 32 0
6 0 0 0 1 1 0 0

Mais cela veut dire que pour un score élevé, le taux de bonne classification s’améliore.

import numpy
ind = numpy.max(pipe.predict_proba(X_test), axis=1) >= 0.6
confusion_matrix_df(y_test[ind], pipe.predict(X_test)[ind])
3 4 5 6 7 8
0 0 0 1 2 0 0
1 0 5 28 4 0 0
2 0 0 528 76 0 0
3 0 0 90 521 13 0
4 0 0 3 60 163 0
5 0 0 0 6 7 27

Les petites classes ont disparu : le modèle n’est pas sûr du tout pour les classes 3, 4, 9. On voit aussi que le modèle se trompe souvent d’une note, il serait sans doute plus judicieux de passer à un modèle de régression plutôt que de classification. Cependant, un modèle de régression ne fournit pas de score de confiance. Sans doute serait-il possible d’en construire avec un modèle de détection d’anomalie…

Anomalies

Une anomalie est un point aberrant. Cela revient à dire que sa probabilité qu’un tel événement se reproduise est faible. Un modèle assez connu est EllipticEnvelope. On suppose que si le modèle détecte une anomalie, un modèle de prédiction aura plus de mal à prédire. On réutilise le pipeline précédent en changeant juste la dernière étape.

from sklearn.covariance import EllipticEnvelope

one = Pipeline([
    ("prep", ColumnTransformer([
        ("color", Pipeline([
            ('one', OneHotEncoder()),
            ('select', ColumnTransformer([('sel1', 'passthrough', [0])]))
        ]), ['color']),
        ("others", "passthrough", numeric_features)
    ])),
    ("lr", EllipticEnvelope()),
])

one.fit(X_train)
ano = one.predict(X_test)
ano
array([ 1,  1,  1, ..., -1, -1,  1])
from pandas import DataFrame
df = DataFrame(dict(note=y_test, ano=one.decision_function(X_test),
                    pred=pipe.predict(X_test),
                    errors=y_test == pipe.predict(X_test),
                    proba_max=numpy.max(pipe.predict_proba(X_test), axis=1),
                ))
df["anoclip"] = df.ano.apply(lambda x: max(x, -200))
df.head()
note ano pred errors proba_max anoclip
2172 7 80.129610 7 True 0.62 80.129610
2662 5 73.460879 5 True 0.61 73.460879
5093 6 85.916048 6 True 0.45 85.916048
2093 6 85.048457 6 True 0.52 85.048457
619 5 -10.361647 6 False 0.63 -10.361647
import seaborn
seaborn.lmplot(x="anoclip", y="proba_max", hue="errors",
               truncate=True, height=5, data=df,
               logx=True);
c:python372_x64libsite-packagesseabornregression.py:279: RuntimeWarning: invalid value encountered in log
  grid = np.c_[np.ones(len(grid)), np.log(grid)]
c:python372_x64libsite-packagesseabornregression.py:282: RuntimeWarning: invalid value encountered in log
  _x = np.c_[_x[:, 0], np.log(_x[:, 1])]
c:python372_x64libsite-packagesnumpylinalglinalg.py:1974: RuntimeWarning: invalid value encountered in greater
  large = s > cutoff
../_images/td2a_correction_cl_reg_anomaly_42_1.png
df.corr()
note ano pred errors proba_max anoclip
note 1.000000 0.083087 0.638612 -0.033520 -0.073341 0.152008
ano 0.083087 1.000000 0.113181 0.030611 -0.040923 0.572297
pred 0.638612 0.113181 1.000000 -0.019382 -0.107277 0.208969
errors -0.033520 0.030611 -0.019382 1.000000 0.370140 0.048471
proba_max -0.073341 -0.040923 -0.107277 0.370140 1.000000 0.008876
anoclip 0.152008 0.572297 0.208969 0.048471 0.008876 1.000000

Les résultats précédents ne sont pas probants. On peut changer de modèle de détection d’anomalies mais les conclusions restent les mêmes. Le score d’anomalie n’est pas relié au score de prédiction.

fig, ax = plt.subplots(1, 2, figsize=(14, 4))
df.proba_max.hist(ax=ax[0], bins=50)
df.anoclip.hist(ax=ax[1], bins=50)
ax[0].set_title("Distribution du score de classification")
ax[1].set_title("Distribution du score d'anomalie");
../_images/td2a_correction_cl_reg_anomaly_45_0.png

C’est joli mais ils n’ont rien à voir. Et c’était prévisible car le modèle de prédiction qu’on utilise est tout-à-fait capable de prédire ce qu’est une anomalie.

pipe_ano = Pipeline([
    ("prep", ColumnTransformer([
        ("color", Pipeline([
            ('one', OneHotEncoder()),
            ('select', ColumnTransformer([('sel1', 'passthrough', [0])]))
        ]), ['color']),
        ("others", "passthrough", numeric_features)
    ])),
    ("lr", RandomForestClassifier()),
])

pipe_ano.fit(X_train, one.predict(X_train))
confusion_matrix_df(one.predict(X_test), pipe_ano.predict(X_test))
-1 1
0 318 36
1 25 2870

Le modèle d’anomalie n’apporte donc aucune information nouvelle. Cela signifie que le modèle prédictif initial n’améliorerait pas sa prédiction en utilisant le score d’anomalie. Il n’y a donc aucune chance que les erreurs ou les score de prédiction soient corrélés au score d’anomalie d’une manière ou d’une autre.

Score de confiance pour une régression

from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import r2_score

pipe_reg = Pipeline([
    ("prep", ColumnTransformer([
        ("color", Pipeline([
            ('one', OneHotEncoder()),
            ('select', ColumnTransformer([('sel1', 'passthrough', [0])]))
        ]), ['color']),
        ("others", "passthrough", numeric_features)
    ])),
    ("lr", RandomForestRegressor()),
])

pipe_reg.fit(X_train, y_train)
r2_score(y_test, pipe_reg.predict(X_test))
0.4761770180844783

Pas super. Mais…

error = y_test - pipe_reg.predict(X_test)
score = numpy.max(pipe.predict_proba(X_test), axis=1)
fig, ax = plt.subplots(1, 2, figsize=(12, 4))
seaborn.kdeplot(score, error, ax=ax[1])
ax[1].set_ylim([-1.5, 1.5])
ax[1].set_title("Densité")
ax[0].plot(score, error, ".")
ax[0].set_xlabel("score de confiance du classifieur")
ax[0].set_ylabel("Erreur de prédiction")
ax[0].set_title("Lien entre classification et prédiction");
../_images/td2a_correction_cl_reg_anomaly_53_0.png

Comme prévu le modèle ne se trompe pas plus dans un sens que dans l’autre.

fig, ax = plt.subplots(1, 2, figsize=(12, 4))
seaborn.kdeplot(score, error.abs(), ax=ax[1])
ax[1].set_ylim([0, 1.5])
ax[1].set_title("Densité")
ax[0].plot(score, error.abs(), ".")
ax[0].set_xlabel("score de confiance du classifieur")
ax[0].set_ylabel("Erreur de prédiction")
ax[0].set_title("Lien entre classification et prédiction");
../_images/td2a_correction_cl_reg_anomaly_55_0.png

On vient de construire un indicateur de confiance pour un modèle de régression. On voit aussi que l’erreur de prédiction est concentrée autour de 0.5 lorsque le score est faible. C’est normal dans la mesure où la probabilité est faible lorsque le classifieur n’est pas sûr, c’est-à-dire que l’observation à prédire est proche de la frontière entre deux classes. Ces classes sont centrées autour des notes entières, les frontières sont au milieu, soit approximativement 3.5, 4.5, … Ce n’est pas une preuve mais seulement la vérifie que l’intervalle de confiance qu’on vient de fabriquer n’est pas complètement aberrant.