2A.ml - Statistiques descriptives avec scikit-learn

Links: notebook, html, PDF, 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
La base de données est déjà renseignée.

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.

Régressions linéaires avec Statsmodels et Scikit-Learn 2A.ml - Statistiques descriptives avec scikit-learn - correction