Analyse du revenu ou plutôt du dépassement d'un seuil de revenu sur des données issues d'un sondage aux USA. Modélisation et prévision du dépassement d'un seuil de revenu. Comparaison de la pertinence et de l'efficacité de différentes méthodes de modélisation ou apprentissage.
Des données publiques disponibles sur le site UCI repository sont extraites de la base de données issue du recensement réalisé aux Etats Unis en 1994. Ces données son largement utilisées et font référence comme outil de benchmark pour comparer les performances de méthodes d’apprentissage ou modélisation statistique. L’objectif est de prévoir la variable binaire « revenu annuel » (income) supérieur ou inférieur à 50k$. Il ne s’agit pas encore de données massives mais, en concaténant les fichiers fournis d'apprentissage et de test, 48841 individus sont décrits par les 14 variables du tableau ci-dessous. La phase préliminaire de préparation et exploration des données est décrite dans un scénario de la Saison 2: Exploration. Les codes de préparation des données sont rapidement repris afin de démarrer du même jeu de données.
| Num | Libellé | Ensemble de valeurs |
|---|---|---|
| 1 | Age |
real |
| 2 | workClass |
Private, Self-emp-not-inc, Self-emp-inc, Federal-gov, Local-gov, State-gov, Without-pay, Never-worked |
| 3 | fnlwgt |
real |
| 4 | education |
Bachelors, Some-college, 11th, HS-grad, Prof-school, Assoc-acdm, Assoc-voc, 9th, 7th-8th, 12th, Masters, 1st-4th, 10th, Doctorate, 5th-6th, Preschool |
| 5 | educNum |
integer |
| 6 | mariStat |
Married-civ-spouse, Divorced, Never-married, Separated, Widowed, Married-spouse-absent, Married-AF-spouse |
| 7 | occup |
Tech-support, Craft-repair, Other-service, Sales, Exec-managerial, Prof-specialty, Handlers-cleaners, Machine-op-inspct, Adm-clerical, Farming-fishing, Transport-moving, Priv-house-serv, Protective-serv, Armed-Forces |
| 8 | relationship |
Wife, Own-child, Husband, Not-in-family, Other-relative, Unmarried |
| 9 | origEthn |
White, Asian-Pac-Islander, Amer-Indian-Eskimo, Other, Black |
| 10 | sex |
Female, Male |
| 11 | capitalGain |
real |
| 12 | capitalLoss |
real |
| 13 | hoursWeek |
real |
| 14 | nativCountry |
United-States, Cambodia, England, Puerto-Rico, Canada, Germany, Outlying-US(Guam-USVI-etc), India, Japan, Greece, South, China, Cuba, Iran, Honduras, Philippines, Italy, Poland, Jamaica, Vietnam, Mexico, Portugal, Ireland, France, Dominican-Republic, Laos, Ecuador, Taiwan, Haiti, Columbia, Hungary, Guatemala, Nicaragua, Scotland, Thailand, Yugoslavia, El-Salvador, Trinadad&Tobago, Peru, Hong, Holand-Netherlands |
| 15 | income |
incHigh (>50K), incLow (<=50K) |
Une première étape permettant de vérifier, sélectionner, recoder les données, a permis de construire un fichier de type .csv qui est utilisé dans ce calepin.
Répondre aux questions en s'aidant des résultats des exécutions.
In [ ]:
%matplotlib inline
# Importations
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
adult=pd.read_csv('adultTrainTest.csv')
adult.head()
Un travail important de nettoyage des données doit être réalisé pour supprimer les données manquantes, éviter des erreurs de codage, regrouper des modalités bien trop nombreuses, éliminer des variables redondantes. Cette étape de préparation et d'exploration est réalisée dans la saison 2- exploration. Seuls sont reproduits ci-dessous les codes permmant la rpéapration du fichier.
In [ ]:
def create_categorical_data(df, column_name):
cat_columns = pd.Categorical(df[column_name], ordered=False)
return cat_columns
In [ ]:
print(np.sort(adult["mariStat"].unique()))
In [ ]:
# mariStat
cat_name_dic = {" Never-married": "Never-Married", " Married-AF-spouse": "Married",
" Married-civ-spouse": "Married", " Married-spouse-absent": "Not-Married",
" Separated": "Not-Married", " Divorced": "Not-Married", " Widowed": "Widowed"}
adult['mariStat'] = adult.mariStat.map(cat_name_dic)
In [ ]:
print(np.sort(adult["nativCountry"].unique()))
In [ ]:
# nativCountry
cat_country = {" Cambodia": "SE-Asia", " Canada": "British-Commonwealth", " China": "China", " Columbia": "South-America",
" Cuba": "Other", " Dominican-Republic": "Latin-America", " Ecuador": "South-America",
" El-Salvador": "South-America", " England": "British-Commonwealth", " France": "Euro_1",
" Germany": "Euro_1", " Greece": "Euro_2", " Guatemala": "Latin-America", " Haiti": "Latin-America",
" Holand-Netherlands": "Euro_1", " Honduras": "Latin-America", " Hong": "China", " Hungary": "Euro_2",
" India": "British-Commonwealth", " Iran": "Other", " Ireland": "British-Commonwealth", " Italy": "Euro_1",
" Jamaica": "Latin-America", " Japan": "Other", " Laos": "SE-Asia", " Mexico": "Latin-America",
" Nicaragua": "Latin-America", " Outlying-US(Guam-USVI-etc)": "Latin-America", " Peru": "South-America",
" Philippines": "SE-Asia", " Poland": "Euro_2", " Portugal": "Euro_2", " Puerto-Rico": "Latin-America",
" Scotland": "British-Commonwealth", " South": "Euro_2", " Taiwan": "China", " Thailand": "SE-Asia",
" Trinadad&Tobago": "Latin-America", " Vietnam": "SE-Asia", " United-States": "United-States",
" Yugoslavia": "Euro_2"}
adult["nativCountry"] = adult.nativCountry.map(cat_country)
In [ ]:
print(np.sort(adult["education"].unique()))
In [ ]:
# education
cat_educ = {" 10th": "Dropout", " 11th": "Dropout", " 12th": "Dropout", " 1st-4th": "Dropout", " 5th-6th": "Dropout",
" 7th-8th": "Dropout", " 9th": "Dropout", " Assoc-acdm": "Associates", " Assoc-voc": "Associates",
" Bachelors": "Bachelors", " Doctorate": "Doctorate", " HS-grad": "HS-grad", " Masters": "Masters",
" Preschool": "Dropout", " Prof-school": "Prof-School", " Some-college": "HS-Graduate"}
adult["education"] = adult.education.map(cat_educ)
In [ ]:
print(np.sort(adult["workClass"].unique()))
In [ ]:
# workClass
cat_work = {" Federal-gov": "Federal-Govt", " Local-gov": "Other-Govt", " State-gov": "Other-Govt", " Private": "Private",
" Self-emp-inc": "Self-Employed", " Self-emp-not-inc": "Self-Employed", " Without-pay": "Not-Working",
" Never-worked": "Not-Working"}
adult["workClass"] = adult.workClass.map(cat_work)
In [ ]:
print(np.sort(adult["occup"].unique()))
In [ ]:
# occup
cat_occup = {" Adm-clerical": "Admin", " Craft-repair": "Blue-Collar", " Exec-managerial": "White-Collar",
" Farming-fishing": "Blue-Collar", " Handlers-cleaners": "Blue-Collar", " Machine-op-inspct": "Blue-Collar",
" Other-service": "Service", " Priv-house-serv": "Service", " Prof-specialty": "Professional",
" Protective-serv": "Other-occups", " Sales": "Sales", " Tech-support": "Other-occups",
" Transport-moving": "Blue-Collar"}
adult["occup"] = adult.occup.map(cat_occup)
In [ ]:
print(np.sort(adult["origEthn"].unique()))
In [ ]:
# origEthn
cat_orig = {" White": "CaucYes", " Black": "CaucNo", " Amer-Indian-Eskimo": "CaucNo", " Asian-Pac-Islander": "CaucNo",
" Other": "CaucNo"}
adult["origEthn"] = adult.origEthn.map(cat_orig)
In [ ]:
adult["LcapitalGain"] = np.log(1 + adult["capitalGain"])
adult["LcapitalLoss"] = np.log(1 + adult["capitalLoss"])
# capital
def quantileCapitalGain(capital):
if type(capital) != int:
result = np.nan
elif capital <= 0:
result = "None"
elif capital <= np.median(adult[adult["capitalGain"] > 0]["capitalGain"]):
result = "cgLow"
else:
result = "cgHigh"
return result
adult["capitalGain"] = list(map(quantileCapitalGain, adult.capitalGain))
def quantileCapitalLoss(capital):
if type(capital) != int:
result = np.nan
elif capital <= 0:
result = "None"
elif capital <= np.median(adult[adult["capitalLoss"] > 0]["capitalLoss"]):
result = "clLow"
else:
result = "clHigh"
return result
adult["capitalLoss"] = list(map(quantileCapitalLoss, adult.capitalLoss))
In [ ]:
adult["ageQ"] = pd.qcut(adult.age, 5, labels=["Ag1", "Ag2", "Ag3", "Ag4", "Ag5"])
In [ ]:
adult["hoursWeekQ"] = pd.cut(adult.hoursWeek, bins=np.array([0, 39, 41, 100]), labels=["HW1", "HW2", "HW3"])
In [ ]:
def create_categorical_data_rename(df, column_name, cat_name_dic):
cat_columns = pd.Categorical(df[column_name], ordered=False)
new_categorie = [cat_name_dic[old_name] for old_name in cat_columns.categories]
return cat_columns.rename_categories(new_categorie)
In [ ]:
print(np.sort(adult["income"].unique()))
In [ ]:
adult["income"] = create_categorical_data_rename(adult, "income", {" <=50K": "incLow", " >50K": "incHigh"})
In [ ]:
for name in ["workClass", "education", "mariStat", "occup", "relationship", "origEthn", "sex", "capitalGain",
"capitalLoss", "nativCountry"]:
adult[name] = create_categorical_data(adult, name)
In [ ]:
adult = adult[np.logical_not(adult.isnull().any(axis=1))]
adult.head()
In [ ]:
adult = adult[(adult["sex"] != "Female") | (adult["relationship"] != "Husband")]
adult = adult[(adult["sex"] != "Male") | (adult["relationship"] != "Wife")]
In [ ]:
adult.describe()
Q Que dire de la distribution de la variable age, de celle income ?
In [ ]:
adult["age"].hist()
plt.show()
In [ ]:
adult["fnlwgt"].hist()
plt.show()
In [ ]:
adult["income"].value_counts()
In [ ]:
adult["relationship"].value_counts()
In [ ]:
adult.plot(kind="scatter",x="age",y="educNum")
plt.show()
Q Que dire des liaisons : age x hoursWeek, age x income, sex x income ?
In [ ]:
adult.plot(kind="scatter",x="hoursWeek",y="age")
plt.show()
In [ ]:
adult.boxplot(column="age",by="income")
plt.show()
Q La variable fnlwgt (Final sampling weight: Inverse of sampling fraction adjusted for non-response and over or under sampling of particular groups) est assez obscure. Que dire de sa liaison avec la variable cible? Elle est supprimée par la suite
In [ ]:
adult.boxplot(column="fnlwgt",by="income")
plt.show()
In [ ]:
# Mosaic plots
from statsmodels.graphics.mosaicplot import mosaic
mosaic(adult,["income","sex"])
plt.show()
In [ ]:
from statsmodels.graphics.mosaicplot import mosaic
mosaic(adult,["income","origEthn"])
plt.show()
Quelques modifications comlémentaires sont apportées de la base. Certaines variables en versions quantitatives et qualitatives comme le nombre d'heures par semaine, l'âge ou le niveau d'éducation sont conservées. Des variables sont supprimées afin de ne garder qu'une seule présence d'une information sensible: genre et origine ethnique.
fnlwgt qui n'a guère de signification pour cette analyse.Child: présence ou non d'enfants.relationship redondante avec le genre et le statut marital,nativCountry redondante avec l'origine ethnique.
In [ ]:
print(np.sort(adult["relationship"].unique()))
In [ ]:
cat_orig = {' Husband':"ChildNo",' Not-in-family':"ChildNo",' Other-relative':"ChildNo",' Own-child':"ChildYes",' Unmarried':"ChildNo",' Wife':"ChildNo"}
adult["child"] = adult.relationship.map(cat_orig)
In [ ]:
adult=adult.drop(["fnlwgt","nativCountry","relationship"],axis=1)
In [ ]:
adult.head()
In [ ]:
adultDum=pd.get_dummies(adult[["workClass","education","mariStat",
"occup","origEthn","sex","capitalGain","capitalLoss","ageQ","hoursWeekQ","child"]])
adultDum.head()
In [ ]:
adultJoin = adult[["age","educNum","hoursWeek","LcapitalGain","LcapitalLoss","income"]].join(adultDum)
In [ ]:
ind_ech = np.random.choice(adultJoin.index.values, 20000, replace=False)
adultEch=adultJoin.loc[ind_ech]
# Variable cible
Y=adultEch["income"]
# Variables prédictives
X=adultEch.drop(["income"],axis=1)
In [ ]:
from sklearn.model_selection import train_test_split
X_train,X_test,Y_train,Y_test=train_test_split(X,Y,test_size=5000,random_state=11)
Q Commenter les options de la commande GridSearchCV. A quoi sert param ?
In [ ]:
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import GridSearchCV
import time
tps0=time.perf_counter()
# Optimisation du paramètre de pénalisation
# grille de valeurs
param=[{"C":[0.8,0.9,1,1.1,1.2]}]
logit = GridSearchCV(LogisticRegression(penalty="l1",solver="liblinear"), param,cv=10,n_jobs=-1)
logitOpt=logit.fit(X_train, Y_train) # GridSearchCV est lui même un estimateur
# paramètre optimal
logitOpt.best_params_["C"]
tps1=(time.perf_counter()-tps0)
print("Temps logit = %f, Meilleur taux = %f, Meilleur paramètre = %s" % (tps1,
1.-logitOpt.best_score_,logitOpt.best_params_))
In [ ]:
# erreur sur l'échantillon test
1-logitOpt.score(X_test, Y_test)
In [ ]:
# Prévision
y_chap = logitOpt.predict(X_test)
# matrice de confusion
table=pd.crosstab(y_chap,Y_test)
print(table)
Q La matrice de confusion n’est pas symétrique. Quelle pourrait en être la raison ?
Q Quels algorithmes pourraient être exécutés en R pour la régression logistique? Quelles informations complémentaires en tirer?
In [ ]:
# Coefficients
LogisticRegression(penalty="l1",C=logitOpt.best_params_['C'],
solver='liblinear').fit(X_train, Y_train).coef_
In [ ]:
from sklearn.tree import DecisionTreeClassifier
tps0=time.perf_counter()
# Optimisation de la profondeur de l'arbre
param=[{"max_depth":list(range(2,10))}]
tree= GridSearchCV(DecisionTreeClassifier(),param,cv=10,n_jobs=-1)
treeOpt=tree.fit(X_train, Y_train)
# paramètre optimal
tps1=(time.perf_counter()-tps0)
print("Temps arbre = %f, Meilleur taux = %f, Meilleur paramètre = %s" % (tps1,
1. - treeOpt.best_score_,treeOpt.best_params_))
In [ ]:
# Estimation de l'erreur de prévision
1-treeOpt.score(X_test,Y_test)
In [ ]:
# prévision de l'échantillon test
y_chap = treeOpt.predict(X_test)
# matrice de confusion
table=pd.crosstab(y_chap,Y_test)
print(table)
In [ ]:
treeOpt.best_params_['max_depth']
In [ ]:
from sklearn.tree import export_graphviz
from sklearn.externals.six import StringIO
import pydotplus
treeG=DecisionTreeClassifier(max_depth=treeOpt.best_params_['max_depth'])
treeG.fit(X_train,Y_train)
dot_data = StringIO()
export_graphviz(treeG, out_file=dot_data)
graph=pydotplus.graph_from_dot_data(dot_data.getvalue())
graph.write_png("treeOpt.png")
In [ ]:
from IPython.display import Image
Image(filename='treeOpt.png')
Q Quelle est l’insuffisance de l’implémentation des arbres de décision dans Scikit-learn par rapport à celle de rpart de R ? Que dire de l’arbre?
In [ ]:
from sklearn.neural_network import MLPClassifier
from sklearn.preprocessing import StandardScaler
# L'algorithme ds réseaux de neurones nécessite éventuellement une normalisation
# des variables explicatives avec les commandes ci-dessous
scaler = StandardScaler()
scaler.fit(X_train.astype(float))
Xnet_train = scaler.transform(X_train.astype(float))
# Meme transformation sur le test
Xnet_test = scaler.transform(X_test.astype(float))
In [ ]:
tps0=time.perf_counter()
param_grid=[{"hidden_layer_sizes":list([(4,),(5,),(6,)])}]
nnet= GridSearchCV(MLPClassifier(max_iter=500),param_grid,cv=10,n_jobs=-1)
nnetOpt=nnet.fit(Xnet_train, Y_train)
# paramètre optimal
tps1=(time.perf_counter()-tps0)
print("Temps perceptron = %f, Meilleur taux = %f, Meilleur paramètre = %s" % (tps1,
1. - nnetOpt.best_score_,nnetOpt.best_params_))
Q Quelle stratégie d’optimisation est adoptée ? Quelle autre pourrait l’être? Quel réseau pourrait également être pris en compte? Quelles sont les fonctions d’activation des neurones?
In [ ]:
# Estimation de l'erreur de prévision sur le test
1-nnetOpt.score(Xnet_test,Y_test)
In [ ]:
# prévision de l'échantillon test
y_chap = nnetOpt.predict(Xnet_test)
# matrice de confusion
table=pd.crosstab(y_chap,Y_test)
print(table)
Q Commenter les choix de tous les paramètres.
In [ ]:
from sklearn.ensemble import RandomForestClassifier
# définition des paramètres
forest = RandomForestClassifier(n_estimators=500,
criterion='gini', max_depth=None,
min_samples_split=2, min_samples_leaf=1,
max_features='auto', max_leaf_nodes=None,
bootstrap=True, oob_score=True)
# apprentissage
rfFit = forest.fit(X_train,Y_train)
print(1-rfFit.oob_score_)
In [ ]:
# erreur de prévision sur le test
1-rfFit.score(X_test,Y_test)
In [ ]:
# optimisation du paramètre
tps0=time.perf_counter()
param=[{"max_features":list(range(2,10,1))}]
rf= GridSearchCV(RandomForestClassifier(n_estimators=100),param,cv=10,n_jobs=-1)
rfOpt=rf.fit(X_train, Y_train)
# paramètre optimal
tps1=(time.perf_counter()-tps0)
print("Temps r forest = %f, Meilleur taux = %f, Meilleur paramètre = %s" % (tps1,
1. - rfOpt.best_score_,rfOpt.best_params_))
In [ ]:
# erreur de prévision sur le test
1-rfOpt.score(X_test,Y_test)
In [ ]:
# prévision
y_chap = rfFit.predict(X_test)
# matrice de confusion
table=pd.crosstab(y_chap,Y_test)
print(table)
In [ ]:
rf= RandomForestClassifier(n_estimators=100,max_features=6)
rfFit=rf.fit(X_train, Y_train)
# Importance décroissante des variables
importances = rfFit.feature_importances_
indices = np.argsort(importances)[::-1]
for f in range(20):
print(X_train.columns[indices[f]], importances[indices[f]])
In [ ]:
# Graphe des importances
plt.figure()
plt.title("Importances des variables")
plt.bar(range(X_train.shape[1]), importances[indices])
plt.xticks(range(X_train.shape[1]), indices)
plt.xlim([-1, X_train.shape[1]])
plt.show()
Q Comment est obtenu le graphique ? Quelle importance ? Comment interpréter ces résultats ?
Q Pourquoi pas de paramètre njobs=-1?
Q En plus de celui optimiser, quels sont les 2 principaux paramètres cet algorithme laissés par défaut ?
In [ ]:
from sklearn.ensemble import GradientBoostingClassifier
tps0=time.perf_counter()
param=[{"n_estimators":[200, 250, 300]}]
gbm= GridSearchCV(GradientBoostingClassifier(),param,cv=10)
gbmOpt=gbm.fit(X_train, Y_train)
# paramètre optimal
tps1=(time.perf_counter()-tps0)
print("Temps boosting = %f, Meilleur taux = %f, Meilleur paramètre = %s" % (tps1,
1. - gbmOpt.best_score_,gbmOpt.best_params_))
In [ ]:
# erreur de prévision sur le test
1-gbmOpt.score(X_test,Y_test)
In [ ]:
# prévision de l'échantillon test
y_chap = gbmOpt.predict(X_test)
# matrice de confusion
table=pd.crosstab(y_chap,Y_test)
print(table)
Q En cohérence avec les résultats précédents, quelle est la courbe la plus au dessus des autres? Commenter ce graphique, que signifie AUC ?
In [ ]:
from sklearn.metrics import roc_curve
listMethod=[["GBM",gbmOpt],["RF",rfOpt],["NN",nnetOpt],["Tree",treeOpt],["Logit",logitOpt]]
In [ ]:
# Courbes ROC des précédents modèles optimaux
for method in enumerate(listMethod):
probas_ = method[1][1].fit(Xnet_train, Y_train).predict_proba(Xnet_test)
fpr, tpr, thresholds = roc_curve(Y_test,probas_[:,1], pos_label="incLow")
plt.plot(fpr, tpr, lw=1,label="%s"%method[1][0]),
plt.xlabel('Taux de faux positifs')
plt.ylabel('Taux de vrais positifs')
plt.legend(loc="best")
plt.show()
Q Quelle différence entre la validation croisée Monte Carlo et la V-fold cross validation?
Q Les variables sont « standardisées ». Pourquoi ? Est-ce important et pour quelles méthodes? Commenter les résultats. Quelle méthode choisir? Quel autre algorithme serait-il pertinent de tester?
In [ ]:
from sklearn.utils import check_random_state
tps0=time.perf_counter()
check_random_state(11)
# définition des estimateurs
logit= LogisticRegression(penalty="l1",solver="liblinear")
tree = DecisionTreeClassifier()
nnet = MLPClassifier(max_iter=400)
rf = RandomForestClassifier(n_estimators=200)
gbm = GradientBoostingClassifier()
# Nombre d'itérations
B=10 # pour utiliser le programme, mettre plutôt B=10
# définition des grilles de paramètres
listMethGrid=[[gbm,{"n_estimators":[200, 250, 300]}],
[rf,{"max_features":list(range(5,10))}],
[nnet,{"hidden_layer_sizes":list([(3,),(4,),(5,)])}],
[tree,{"max_depth":list(range(5,10))}],
[logit,{"C":[0.8,0.9,1,1.1,1.2]}]]
# Initialisation à 0
arrayErreur=np.empty((B,5))
for i in range(B): # itérations sur B échantillons test
# extraction apprentissage et test
X_train,X_test,Y_train,Y_test=train_test_split(X,Y,test_size=1000)
scaler = StandardScaler()
scaler.fit(X_train.astype(float))
Xnet_train = scaler.transform(X_train.astype(float))
# Meme transformation sur le test
Xnet_test = scaler.transform(X_test.astype(float))
# optimisation de chaque méthode et calcul de l'erreur sur le test
for j,(method, grid_list) in enumerate(listMethGrid):
methodGrid=GridSearchCV(method,grid_list,cv=10,n_jobs=-1).fit(X_train, Y_train)
methodOpt = methodGrid.best_estimator_
methFit=methodOpt.fit(Xnet_train, Y_train)
arrayErreur[i,j]=1-methFit.score(Xnet_test,Y_test)
tps1=time.perf_counter()
In [ ]:
dataframeErreur=pd.DataFrame(arrayErreur,columns=["GBM","RF","NN","Tree","Logit"])
print("Temps execution :",(tps1 - tps0))
In [ ]:
dataframeErreur[["GBM","RF","NN","Tree","Logit"]].boxplot(return_type='dict')
plt.show()
In [ ]:
# Moyennes
dataframeErreur.mean()
Q L'algoritme GBM n'a pas été complètement optimisé. Une optimisation fine de l'algorithme XGBoost permettra-t-elle une meilleure performance que celle de random forest?
Q Les SVM ne font pas partie de la comparaison à cause de temps rédhibitoires d’exécution. A quoi est-ce dû ? Commenter les temps d’exécution des différentes étapes.
Q Les temps d'exécution deviennent longs pour un ordinateur de bureau standard. Comment évaluer l'impact de la taille de l'échantillon sur le temps d'exécution et la précision des prévisions?
Q Est-il réellement nécessaire d'entraîner de type d'algorithme et pour ce type de problème sur de gros jeux de données?