2A.ml - Statistiques descriptives avec scikit-learn#

Links: notebook, html, python, slides, GitHub

ACP, CAH, régression lineaire.

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

Introduction#

Les statistiques descriptives sont abordées de la première année à l’ENSAE. Un des livres que je consulte souvent est celui de Gilles Saporta : Probabilités, analyse des données et statistique qui est en français.

Le module scikit-learn a largement contribué au succès de Python dans le domaine du machine learning. Ce module inclut de nombreuses techniques regroupées sous le terme statistiques descriptives. La correspondance anglais-français n’est pas toujours évidente. Voici quelques termes :

scikit-learn est orienté machine learning, les résultats qu’il produit sont un peu moins complets que statsmodels pour les modèles statistiques linéaires ou fastcluster pour la CAH. L’objectif de ces deux heures est d’utiliser ces modules pour étudier un jeu de données :

ACP (Analyse en Composantes Principales)#

Le site data.gouv.fr propose de nombreux jeux de données dont Séries chronologiques Education : les élèves du second degré. Ces données sont également accessibles comme ceci :

import pandas, numpy, pyensae.datasource
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt

fichier = pyensae.datasource.download_data("eleve_region.txt")
df = pandas.read_csv("eleve_region.txt", sep="\t", encoding="utf8", index_col=0)
print(df.shape)
df.head(n=5)
(27, 21)
1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 ... 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013
académie
Aix-Marseille 241357 242298 242096 242295 243660 244608 245536 247288 249331 250871 ... 250622 248208 245755 243832 242309 240664 240432 241336 239051 240115
Amiens 198281 196871 195709 194055 192893 191862 189636 185977 183357 180973 ... 175610 172110 168718 165295 163116 162548 163270 164422 165275 166345
Besançon 116373 115600 114282 113312 112076 110261 108106 105463 103336 102264 ... 100117 98611 97038 95779 95074 94501 94599 94745 94351 94613
Bordeaux 253551 252644 249658 247708 247499 245757 244992 243047 243592 245198 ... 244805 244343 242602 242933 243146 244336 246806 250626 252085 255761
Caen 145435 144369 141883 140658 139585 137704 135613 133255 131206 129271 ... 125552 123889 122550 121002 119857 119426 119184 119764 119010 119238

5 rows × 21 columns

On veut observer l’évolution du nombre d’élèves. On prend l’année 1993 comme base 100.

for c in df.columns:
    if c != "1993":
        df[c] /= df ["1993"]
df["1993"] /= df["1993"]
df.head()
1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 ... 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013
académie
Aix-Marseille 1.0 1.003899 1.003062 1.003886 1.009542 1.013470 1.017315 1.024574 1.033038 1.039419 ... 1.038387 1.028385 1.018222 1.010255 1.003944 0.997129 0.996168 0.999913 0.990446 0.994854
Amiens 1.0 0.992889 0.987029 0.978687 0.972826 0.967627 0.956400 0.937947 0.924733 0.912710 ... 0.885662 0.868011 0.850904 0.833640 0.822651 0.819786 0.823427 0.829237 0.833539 0.838936
Besançon 1.0 0.993358 0.982032 0.973697 0.963076 0.947479 0.928961 0.906250 0.887972 0.878761 ... 0.860311 0.847370 0.833853 0.823035 0.816976 0.812053 0.812895 0.814149 0.810764 0.813015
Bordeaux 1.0 0.996423 0.984646 0.976955 0.976131 0.969261 0.966243 0.958572 0.960722 0.967056 ... 0.965506 0.963684 0.956817 0.958123 0.958963 0.963656 0.973398 0.988464 0.994218 1.008716
Caen 1.0 0.992670 0.975577 0.967154 0.959776 0.946842 0.932465 0.916251 0.902162 0.888858 ... 0.863286 0.851851 0.842644 0.832001 0.824128 0.821164 0.819500 0.823488 0.818304 0.819871

5 rows × 21 columns

Il n’est pas évident d’analyser ce tableaux de chiffres. On utilise une ACP pour projeter les académies dans un plan.

pca = PCA(n_components=4)
print(pca.fit(df))
pca.explained_variance_ratio_
PCA(n_components=4)
array([9.63024476e-01, 3.37267932e-02, 1.93116666e-03, 6.61046399e-04])

Le premier axe explique l’essentiel de la variance. Les variables n’ont pas été normalisées car elles évoluent dans les mêmes ordres de grandeur.

plt.bar(numpy.arange(len(pca.explained_variance_ratio_)) + 0.5,
        pca.explained_variance_ratio_)
plt.title("Variance expliquée");
../_images/td2a_cenonce_session_3A_12_0.png

On affiche les académies dans le plan des deux premiers axes :

X_reduced = pca.transform(df)
plt.figure(figsize=(18, 6))
plt.scatter(X_reduced[:, 0], X_reduced[:, 1])

for label, x, y in zip(df.index, X_reduced[:, 0], X_reduced[:, 1]):
    plt.annotate(
        label,
        xy=(x, y), xytext=(-10, 10),
        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.title("ACP rapprochant des profils d'évolution similaires");
../_images/td2a_cenonce_session_3A_14_0.png

Puis on vérifie que deux villes proches ont le même profil d’évolution au cours des années :

sub = df.loc[["Paris", "Bordeaux", "Lyon", "Nice", "La Réunion", "Reims"], :]
ax = sub.transpose().plot(figsize=(10, 4))
ax.set_title("Evolution des effectifs au cours du temps");
../_images/td2a_cenonce_session_3A_16_0.png

L’ACP version statsmodels produit le même type de résultats. Un exemple est disponible ici : PCA and Biplot using Python.

Exercice 1 : CAH (classification ascendante hiérarchique)#

Le point commun de ces méthodes est qu’elles ne sont pas supervisées. L’objectif est de réduire la complexité des données. Réduire le nombre de dimensions pour l’ACP ou segmenter les observations pour les k-means et la CAH. On propose d’utiliser une CAH sur les mêmes données.

Le module scikit-learn.cluster ne propose pas de fonction pour dessiner le dendrogram. Il faudra utiliser celle-ci : dendrogram et sans doute s’inspirer du code suivant.

from sklearn.cluster import AgglomerativeClustering
from scipy.cluster.hierarchy import dendrogram

ward = AgglomerativeClustering(linkage='ward', compute_full_tree=True).fit(df)
dendro = [ ]
for a,b in ward.children_:
    dendro.append([a, b, float(len(dendro)+1), len(dendro)+1])

fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(1, 1, 1)
r = dendrogram(dendro, color_threshold=1, labels=list(df.index),
               show_leaf_counts=True, ax=ax, orientation="left")
../_images/td2a_cenonce_session_3A_19_0.png

Exercice 2 : régression#

Ce sont trois méthodes supervisées : on s’en sert pour expliquer prédire le lien entre deux variables X et Y (ou ensemble de variables) ou prédire Y en fonction de X. Pour cet exercice, on récupère des données relatives aux salaires Salaires et revenus d’activités (les chercher avec la requête insee données dads sur un moteur de recherche). La récupération des données est assez fastidieuse. La première étape consiste à télécharger les données depuis le site de l’INSEE. La seconde étape consiste à convertir les données au format sqlite3. Pour ce fichier, ce travail a déjà été effectué et peut être téléchargé depuis mon site. La base comprend 2 millions de lignes.

import pyensae.datasource
f = pyensae.datasource.download_data("dads2011_gf_salaries11_dbase.zip",
                          website="https://www.insee.fr/fr/statistiques/fichier/2011542/")
f
['.\salaries11.dbf', 'varlist_salaries11.dbf', '.\varmod_salaries11.dbf']
import pandas
try:
    from dbfread import DBF
    use_dbfread = True
except ImportError as e :
    use_dbfread = False

if use_dbfread:
    import os
    from pyensae.sql.database_exception import ExceptionSQL
    from pyensae.datasource import dBase2sqllite
    print("convert dbase into sqllite")
    try:
        dBase2sqllite("salaries2011.db3", "varlist_salaries11.dbf", overwrite_table="varlist")
        dBase2sqllite("salaries2011.db3", "varmod_salaries11.dbf", overwrite_table="varmod")
        dBase2sqllite("salaries2011.db3", 'salaries11.dbf',
                      overwrite_table="salaries", fLOG = print)
    except ExceptionSQL:
        print("La base de données est déjà renseignée.")
else :
    print("use of zipped version")
    import pyensae.datasource
    db3 = pyensae.datasource.download_data("salaries2011.zip")
    # pour aller plus vite, données à télécharger au
    # http://www.xavierdupre.fr/enseignement/complements/salaries2011.zip
convert dbase into sqllite
SQL  'drop table if exists salaries11'
SQL  'create table salaries (A6 TEXT, A17 TEXT, A38 TEXT, REGR TEXT, DEPR TEXT, REGT TEXT, DEPT TEXT, SEXE TEXT, PCS TEXT, CONT_TRAV TEXT, CONV_COLL TEXT, TYP_EMPLOI TEXT, DUREE REAL, DATDEB REAL, DATFIN REAL, CPFD TEXT, DOMEMPL TEXT, DOMEMPL_EM TEXT, FILT TEXT, AGE REAL, CS TEXT, NB_PER REAL, NB_PER_N REAL, NBHEUR REAL, NBHEUR_TOT REAL, TRALCHT TEXT, TREFF TEXT, TRNNETO TEXT, POND REAL)'
moving line  0  to table salaries
moving line  20000  to table salaries
moving line  40000  to table salaries
moving line  60000  to table salaries
moving line  80000  to table salaries
moving line  100000  to table salaries
moving line  120000  to table salaries
moving line  140000  to table salaries
moving line  160000  to table salaries
moving line  180000  to table salaries
moving line  200000  to table salaries
moving line  220000  to table salaries
moving line  240000  to table salaries
moving line  260000  to table salaries
moving line  280000  to table salaries
moving line  300000  to table salaries
moving line  320000  to table salaries
moving line  340000  to table salaries
moving line  360000  to table salaries
moving line  380000  to table salaries
moving line  400000  to table salaries
moving line  420000  to table salaries
moving line  440000  to table salaries
moving line  460000  to table salaries
moving line  480000  to table salaries
moving line  500000  to table salaries
moving line  520000  to table salaries
moving line  540000  to table salaries
moving line  560000  to table salaries
moving line  580000  to table salaries
moving line  600000  to table salaries
moving line  620000  to table salaries
moving line  640000  to table salaries
moving line  660000  to table salaries
moving line  680000  to table salaries
moving line  700000  to table salaries
moving line  720000  to table salaries
moving line  740000  to table salaries
moving line  760000  to table salaries
moving line  780000  to table salaries
moving line  800000  to table salaries
moving line  820000  to table salaries
moving line  840000  to table salaries
moving line  860000  to table salaries
moving line  880000  to table salaries
moving line  900000  to table salaries
moving line  920000  to table salaries
moving line  940000  to table salaries
moving line  960000  to table salaries
moving line  980000  to table salaries
moving line  1000000  to table salaries
moving line  1020000  to table salaries
moving line  1040000  to table salaries
moving line  1060000  to table salaries
moving line  1080000  to table salaries
moving line  1100000  to table salaries
moving line  1120000  to table salaries
moving line  1140000  to table salaries
moving line  1160000  to table salaries
moving line  1180000  to table salaries
moving line  1200000  to table salaries
moving line  1220000  to table salaries
moving line  1240000  to table salaries
moving line  1260000  to table salaries
moving line  1280000  to table salaries
moving line  1300000  to table salaries
moving line  1320000  to table salaries
moving line  1340000  to table salaries
moving line  1360000  to table salaries
moving line  1380000  to table salaries
moving line  1400000  to table salaries
moving line  1420000  to table salaries
moving line  1440000  to table salaries
moving line  1460000  to table salaries
moving line  1480000  to table salaries
moving line  1500000  to table salaries
moving line  1520000  to table salaries
moving line  1540000  to table salaries
moving line  1560000  to table salaries
moving line  1580000  to table salaries
moving line  1600000  to table salaries
moving line  1620000  to table salaries
moving line  1640000  to table salaries
moving line  1660000  to table salaries
moving line  1680000  to table salaries
moving line  1700000  to table salaries
moving line  1720000  to table salaries
moving line  1740000  to table salaries
moving line  1760000  to table salaries
moving line  1780000  to table salaries
moving line  1800000  to table salaries
moving line  1820000  to table salaries
moving line  1840000  to table salaries
moving line  1860000  to table salaries
moving line  1880000  to table salaries
moving line  1900000  to table salaries
moving line  1920000  to table salaries
moving line  1940000  to table salaries
moving line  1960000  to table salaries
moving line  1980000  to table salaries
moving line  2000000  to table salaries
moving line  2020000  to table salaries
moving line  2040000  to table salaries
moving line  2060000  to table salaries
moving line  2080000  to table salaries
moving line  2100000  to table salaries
moving line  2120000  to table salaries
moving line  2140000  to table salaries
moving line  2160000  to table salaries
moving line  2180000  to table salaries
moving line  2200000  to table salaries
moving line  2220000  to table salaries
moving line  2240000  to table salaries

Les données des salaires ne sont pas numériques, elles correspondent à des intervalles qu’on convertit en prenant le milieu de l’intervalle. Pour le dernier, on prend la borne supérieure.

import sqlite3, pandas
con = sqlite3.connect("salaries2011.db3")
df = pandas.io.sql.read_sql("select * from varmod", con)
con.close()

values = df[ df.VARIABLE == "TRNNETO"].copy()

def process_intervalle(s):
    # [14 000 ; 16 000[ euros
    acc = "0123456789;+"
    s0 = "".join(c for c in s if c in acc)
    spl = s0.split(';')
    if len(spl) != 2:
        raise ValueError("Unable to process '{0}'".format(s0))
    try:
        a = float(spl[0])
    except Exception as e:
        raise ValueError("Cannot interpret '{0}' - {1}".format(s, spl))
    b = float(spl[1]) if "+" not in spl[1] else None
    if b is None:
        return a
    else:
        return (a+b) / 2.0

values["montant"] = values.apply(lambda r : process_intervalle(r ["MODLIBELLE"]), axis = 1)
values.head()
VARIABLE MODALITE MODLIBELLE montant
8957 TRNNETO 00 [0 ; 200[ euros 100.0
8958 TRNNETO 01 [200 ; 500[ euros 350.0
8959 TRNNETO 02 [500 ; 1 000[ euros 750.0
8960 TRNNETO 03 [1 000 ; 1 500[ euros 1250.0
8961 TRNNETO 04 [1 500 ; 2 000[ euros 1750.0

On crée la base d’apprentissage :

import sqlite3, pandas
con = sqlite3.connect("salaries2011.db3")
data = pandas.io.sql.read_sql("select TRNNETO,AGE,SEXE from salaries", con)
con.close()

salaires = data.merge ( values, left_on="TRNNETO", right_on="MODALITE" )
salaires.head()
TRNNETO AGE SEXE VARIABLE MODALITE MODLIBELLE montant
0 02 49.0 1 TRNNETO 02 [500 ; 1 000[ euros 750.0
1 02 27.0 1 TRNNETO 02 [500 ; 1 000[ euros 750.0
2 02 22.0 1 TRNNETO 02 [500 ; 1 000[ euros 750.0
3 02 26.0 1 TRNNETO 02 [500 ; 1 000[ euros 750.0
4 02 29.0 2 TRNNETO 02 [500 ; 1 000[ euros 750.0

On récupère les variables utiles pour la régression.

salaires["M"] = salaires.apply( lambda r : 1 if r["SEXE"] == "1" else 0, axis=1)
salaires["F"] = 1 - salaires["M"]  # en supposant que le sexe est toujours renseigné
data = salaires[["AGE", "M", "F", "montant"]]
data = data[data.M + data.F > 0]
data.head()
AGE M F montant
0 49.0 1 0 750.0
1 27.0 1 0 750.0
2 22.0 1 0 750.0
3 26.0 1 0 750.0
4 29.0 0 1 750.0

Ce type d’écriture est plutôt lent car une fonction Python est exécutée à chaque itération. Il est préférable dès que c’est possible d’utiliser les expressions avec des indices sans passer par la fonction apply qui créé une copie de chaque ligne avant d’appliquer la fonction à appliquer à chacune d’entre elles.

salaires["M2"] = 0
salaires.loc[salaires["SEXE"] == "1", "M2"] = 1
salaires["F2"] = 1 - salaires["M2"]
salaires.head()
TRNNETO AGE SEXE VARIABLE MODALITE MODLIBELLE montant M F M2 F2
0 02 49.0 1 TRNNETO 02 [500 ; 1 000[ euros 750.0 1 0 1 0
1 02 27.0 1 TRNNETO 02 [500 ; 1 000[ euros 750.0 1 0 1 0
2 02 22.0 1 TRNNETO 02 [500 ; 1 000[ euros 750.0 1 0 1 0
3 02 26.0 1 TRNNETO 02 [500 ; 1 000[ euros 750.0 1 0 1 0
4 02 29.0 2 TRNNETO 02 [500 ; 1 000[ euros 750.0 0 1 0 1

Il ne reste plus qu’à faire la régression.

version scikit-learn#

Vous pouvez vous inspirer de cet exemple.

version statsmodels#

L’exemple avec ce module est ici.