Mickaël Tits CETIC mickael.tits@cetic.be
Dans ce Chapitre, nous allons analyser un vrai Dataset de biens immobiliers: le "Ames Housing Dataset". A partir de ces données, nous allons développer un modèle prédictif permettant d'estimer le prix d'une maison à partir de nombreuses caractéristiques, telles que sa surface, le nombre de pièces, différents indices de qualité, etc.
Plus d'informations ici: https://www.kaggle.com/c/home-data-for-ml-course/overview/description
http://jse.amstat.org/v19n3/decock.pdf
Détails sur les variables du dataset: https://github.com/titsitits/Python_Data_Science/blob/master/Donn%C3%A9es/data_description.txt
import pandas as pd
import matplotlib.pyplot as plt
ames = pd.read_csv("https://raw.githubusercontent.com/titsitits/Python_Data_Science/master/Donn%C3%A9es/train.csv")
ames
Commençons par explorer brièvement la qualité des données.
On a 1460 observations (maisons), et 81 variables dont un Id et le prix de vente (SalePrice). 19 variables contiennent des données invalides ou manquantes.
#Quelques fonctions utiles pour l'exploration
#Par soucis de lisibilité, on affichera les series comme des dataframes d'une ligne
def display_series(series):
display(series.to_frame().transpose())
#Corrélation entre deux colonnes
def col_corr(df,col1,col2):
return df[[col1,col2]].corr().values[0,1]
#Pour analyser l'effet d'une variable continue sur une autre, on peut extraire deux groupes (chaque côté de la médiane), et afficher un boxplot par groupe
def mediangroups_boxplot_comparison(df, group_col, comparison_col):
df["above_median"] = df[group_col] > df[group_col].quantile(0.5)
df.boxplot(comparison_col, by = "above_median")
df.pop("above_median")
#Pour analyser l'effet d'une variable continue sur une autre, on peut extraire deux groupes (chaque côté de la médiane), et calculer la moyenne par groupe
def mediangroups_mean_comparison(df, group_col, comparison_col):
means = df.groupby(df[group_col] > df[group_col].quantile(0.5))[comparison_col].mean()
means.index = ['Below','Above']
return means
#Idem en séparant les groupes avec la moyenne
def meangroups_mean_comparison(df, group_col, comparison_col):
means = df.groupby(df[group_col] > df[group_col].mean())[comparison_col].mean()
means.index = ['Below','Above']
return means
print(ames.shape)
print(list(ames.columns))
#Colonnes incomplètes
counts = ames.count()
incomplete = counts[counts < len(ames)]
display_series(incomplete.sort_values())
len(incomplete.index)
ames.hist("SalePrice", bins = 30)
sur les 19 variables incomplètes, 3 sont des variables numériques (les variables de type "float64"), et les autres sont des variables catégorielles (les variables de type "object").
Après vérification de la description des variables, on remarque que les variables catégorielles contiennent des "NA" comme catégories, et qui sont automatiquement interprétées comme des NaN par pandas.
Pour les variables catégorielles, nous allons simplement remplacer les NaN par une catégorie "Nothing". Nous allons ensuite analyser et corriger chaque variable numérique.
#Pour vérifier le type d'une colonne: .dtypes()
display_series(ames.dtypes)
#Les variables catégorielles sont du type "object"
#Vérifions les données incomplètes uniquement:
incomplete_types = ames[incomplete.index].dtypes
display_series(incomplete_types)
#Trois variables de type "float64" contiennent des NaN
display_series(incomplete_types[incomplete_types == "float64"])
Pour les variables catégorielles, remplaçons les NaN par la catégorie "Nothing".
cat_columns = ames.columns[ames.dtypes == 'object']
for col in cat_columns:
ames[col] = ames[col].fillna('Nothing')
Analysons la qualité de chaque variable: quelle est la proportion de données manquantes ?
counts = ames.count()
incomplete = counts[counts < len(ames)]
display_series(incomplete.sort_values())
Est-ce que le fait que la variable soit manquante est en soi une information intéressante ? autrement dit, est-ce que ça a une influence sur le prix ?
#Mean price for NaNs and for non-NaNs
display(ames.groupby(ames.LotFrontage.isna())["SalePrice"].agg(['median','mean']))
display(ames.groupby(ames.GarageYrBlt.isna())["SalePrice"].agg(['median','mean']))
display(ames.groupby(ames.MasVnrArea.isna())["SalePrice"].agg(['median','mean']))
L'absence de LotFrontage ne semble pas avoir beaucoup d'influence sur le prix. Les deux autres variables ont une influence. GarageYrBlt = NaN signifie probablement qu'il n'y a pas de garage (ce qui diminue la valeur), et MasVnrArea = NaN semble augmenter la valeur en moyenne.
Pour une comparaison plus complète, nous pouvons visualiser les distributions des valeurs NaNs et non-NaN pour chacune de ces variables:
from matplotlib import pyplot as plt
fig, axes = plt.subplots(1,3)
for ax, feat in zip(axes, ['LotFrontage', 'GarageYrBlt','MasVnrArea']):
ax.violinplot(dataset = [ames[ames[feat].isna()]['SalePrice'].values, ames[ames.notna()]['SalePrice'].values] )
ax.set_xlabel(feat)
ax.set_xticks([1,2])
ax.set_xticklabels(['NaNs','valid'])
#Permet d'éviter l'overlap entre les subplots
plt.tight_layout()
#Analysons la relation entre "MasVnrArea" et le prix
print(col_corr(ames,"MasVnrArea","SalePrice"))
display_series(mediangroups_mean_comparison(ames,"MasVnrArea","SalePrice"))
mediangroups_boxplot_comparison(ames,"MasVnrArea","SalePrice")
ames.plot.scatter("MasVnrArea","SalePrice")
#Il semblerait que l'influence soit positive
#Une possibilité pour prendre en compte cette information est de remplacer les NaN la moyenne de la colonne (ce qui minimisera l'influence de la variable pour ces observations particulières), et d'ajouter une colonne "isnan_MasVnrArea" pour garder l'information
#La relation entre "MasVnrArea" et "SalePrice" ne devrait pas être trop faussée de cette manière, étant donné que le nombre de NaN est assez faible (81).
ames["isnan_MasVnrArea"] = ames.MasVnrArea.isna().astype(int) #.astype(int) permet de rendre la variable numérique. #.astype(str) permettrait d'en faire une variable catégorielle
ames["MasVnrArea"] = ames["MasVnrArea"].fillna(ames["MasVnrArea"].median())
#Analysons la relation entre "LotFrontage" et le prix
print(col_corr(ames,"LotFrontage","SalePrice"))
display_series(mediangroups_mean_comparison(ames,"LotFrontage","SalePrice"))
mediangroups_boxplot_comparison(ames,"LotFrontage","SalePrice")
ames.plot.scatter("LotFrontage","SalePrice", s = 1)
#Ajoutons une variable isnan_LotFrontage pour pouvoir utiliser cette information dans notre modèle plus tard:
ames["isnan_LotFrontage"] = ames.LotFrontage.isna().astype(int)
#Ensuite, étant donné le grand nombre de NaNs (~18%), il est préférable (dans une premier temps) de ne pas utiliser cette variable en remplaçant les NaN par la moyenne
#Une autre possibilité serait d'entraîner un modèle de régression pour prédire ces données manquantes à partir des autres variables.
out = ames.pop("LotFrontage")
On ne peut pas considérer Lotfrontage = NaN comme 0, auquel cas le prix moyen serait plus faible que pour les autres maisons. Il est préférable d'omettre la variable vu la quantité de NaN (plus de 10%). Remplacer NaN par une valeur apporterait des informations erronnées, et omettre les observations réduirait significativement la taille du dataset.
#Concernant l'année de construction du garage, une valeur nulle n'aurait aucun sens.
#On peut vérifier que la NaN correspond simplement à une absence de garage. En effet, leur surface est toujours nulle:
display_series(ames[ames.GarageYrBlt.isna()].GarageArea.describe())
#L'information sur l'absence de garage est donc déjà présente dans une autre variable.
ames.corrwith(ames["GarageYrBlt"]).sort_values(ascending = False).plot.bar(figsize = (10,4))
#On remarque aussi que l'année de construction du garage est logiquement corrélée avec l'année de construction du bien.
#Pour remplir les quelques NaN, on peut éventuellement "simuler" une date de construction du garage avec l'année de construction.
ames["GarageYrBlt"] = ames["GarageYrBlt"].fillna(ames["YearBuilt"])
out = ames.pop("Id")
Il est intéressant de redéfinir le pipeline de correction avec une fonction, pour pouvoir le réappliquer plus tard sur de nouvelles données (le set de test par exemple)
def clean_df(df):
#Fill categorical columns
cat_columns = df.columns[df.dtypes == 'object']
for col in cat_columns:
df[col] = df[col].fillna('Nothing')
#Clean MasVnrArea
df["isnan_MasVnrArea"] = df.MasVnrArea.isna().astype(int)
df["MasVnrArea"] = df["MasVnrArea"].fillna(df["MasVnrArea"].mean())
#Clean LotFrontage
df["isnan_LotFrontage"] = df.LotFrontage.isna().astype(int)
df = df.drop("LotFrontage", axis=1)
#Clean GarageYrBlt
df["GarageYrBlt"] = df["GarageYrBlt"].fillna(df["YearBuilt"])
#Suppression des variables inutiles
out = df.pop("Id")
return df
Remarque: les méthodes pandas.DataFrame.corr() et .corrwith() omettent automatiquement les variables non-numériques.
from matplotlib import pyplot as plt
#corrélations des variables avec le prix
corr_with_price = ames.corrwith(ames["SalePrice"])
#On trie les variables selon la valeur absolue de leur corrélation avec le prix
best_features = corr_with_price.abs().sort_values(ascending=False)
# On ne évidemment peut utiliser le label comme variable prédictive
best_features.pop("SalePrice")
best_features.plot.bar(figsize=(10,4))
best_features = best_features.index.to_list()
#Analyse de l'effet
continuous_col = ames.columns[(ames.dtypes == 'int64') | (ames.dtypes == 'float64')]
delta_mean = []
for col in list(continuous_col):
means = meangroups_mean_comparison(ames, col, "SalePrice")
delta_mean.append( means[1] - means[0] )
best_features_effects = pd.Series(delta_mean, continuous_col).abs().sort_values(ascending = False)
best_features_effects = best_features_effects.drop("SalePrice")
best_features_effects.plot.bar(figsize=(10,4))
best_features_effects = best_features_effects.index.to_list()
Les résultats sont assez similaires. En réalité, les différences principales se trouvent pour les variables dont la distribution est fortement asymmétrique (e.g.: PoolArea est la plupart du temps 0, tout comme isnan_MasVnrArea)
fig, axes = plt.subplots(1,4,figsize = (15,3))
ames.boxplot("OverallQual", ax = axes[0])
ames.boxplot("GrLivArea", ax = axes[1])
ames.boxplot("PoolArea", ax = axes[2])
ames.boxplot("isnan_MasVnrArea", ax = axes[3])
#Prenons les cinq meilleures variables prédictives (en sa basant sur leur corrélation avec le prix, ou le delta de la moyenne)
featurelist = best_features[:5]
#featurelist = best_features_effects[:5]
X = ames[featurelist]
y = ames["SalePrice"]
from sklearn.model_selection import train_test_split
import numpy as np
ids = ames.index
trainids, valids = train_test_split(ids, test_size = 0.4, random_state = 1)
#Pour avoir à chaque fois des ids différents, il suffit de ne pas fixer l'état pseudo alétoire "random_state":
#trainids, valids = train_test_split(ids, test_size = 0.4)
Xtrain, Xval, ytrain, yval = X.loc[trainids], X.loc[valids], y.loc[trainids], y.loc[valids]
#Remarque: ce code est équivalent, mais le but ici est de garder les mêmes ids pour comparer plusieurs modèles (ou pour calculer la moyenne des prédictions, voir fin du notebook)
#Xtrain, Xval, ytrain, yval = train_test_split(X, y, test_size = 0.4)
from sklearn.linear_model import LinearRegression
price_predictor = LinearRegression()
price_predictor.fit(Xtrain, ytrain)
import numpy as np
from sklearn.metrics import mean_absolute_error as MAE, mean_squared_error as MSE
print("R2 Scores (train, val):", price_predictor.score(Xtrain, ytrain), price_predictor.score(Xval, yval))
y_pred = price_predictor.predict(Xval)
print("biais:", np.mean(y_pred - yval.values) ) #si il est négatif: sous-évaluation (en moyenne), si il est positif: sur-évaluation
print("MAE:", MAE(yval.values, y_pred) )
print("RMSE:", np.sqrt(MSE(yval.values, y_pred)) )
#Tester en gardant ou en retirant les outliers; tester avec différentes ensembles de variables
Remarque: les résultats seront un peu différents chaque fois qu'on relancera le code, étant donné la sélection aléatoire des sets d'entraînement et de validation.
from sklearn.preprocessing import StandardScaler
maes = pd.DataFrame(columns = ['train','val'])
for n in range(1,len(best_features)):
#Choix des variables
featurelist = best_features[:n]
#featurelist = best_features_effects[:n]
#Commencer par les plus mauvaises
#featurelist = best_features[::-1][:n]
#featurelist = best_features_effects[::-1][:n]
X = ames[featurelist]
y = ames["SalePrice"]
Xtrain, Xval, ytrain, yval = X.loc[trainids], X.loc[valids], y.loc[trainids], y.loc[valids]
#scaler = StandardScaler()
#Xtrain = scaler.fit_transform(Xtrain)
#Xval = scaler.transform(Xval)
price_predictor.fit(Xtrain, ytrain)
trainpred = price_predictor.predict(Xtrain)
valpred = price_predictor.predict(Xval)
maes.loc[n] = MAE(ytrain, trainpred), MAE(yval, valpred)
maes.plot()
plt.title("Performances (MAE) pour différents nombres de variables prédictives")
plt.xlabel("Nombre de variables prédictives")
plt.ylabel("Maen Absolute Error")
maes['val'].min(), maes['val'].idxmin()
#Meilleur résultat en validation
featurelist = best_features[:28]
X = ames[featurelist]
y = ames["SalePrice"]
Xtrain, Xval, ytrain, yval = X.loc[trainids], X.loc[valids], y.loc[trainids], y.loc[valids]
price_predictor.fit(Xtrain, ytrain)
trainpred = price_predictor.predict(Xtrain)
valpred = price_predictor.predict(Xval)
def display_results(labels, predictions, title = "Model results", figsize = (10,6), **kwargs):
print("Mean average error: ", MAE(labels, predictions))
plt.scatter(labels, predictions, s=2, label = title, **kwargs)
plt.gcf().set_size_inches(figsize)
plt.title(title)
plt.xlabel("Labels")
plt.ylabel("Prédictions")
prev_legend = plt.legend()
plt.legend()
xmin, xmax = plt.gca().get_xlim()
plt.plot([xmin,xmax],[xmin,xmax], 'k')
print(MAE(ytrain, trainpred), MAE(yval, valpred))
display_results(yval, valpred, "Régression linéaire - validation set")
On constate que le nuage de doit semble former une légère courbe. Cela semble indiquer qu'il existe des relations non-linéaires, telles que des relations polynomiales entre les variables et le prix. La régression linéaire se base en effet sur l'hypothèse très simpliste que les relations seraient linéaires. Si on revisualise les nuages de points plus haut montrant l'interactions entre différentes variables (comme 'OverallQual') et le prix, on constate que la relation et plutôt quadratique que linéaire.
Une variation de la régression linéaire, appelée régression polynomiale, permet de facilement prendre en compte de potentielles interactions non-linéaires entre les variables, en créant de nouvelles variables à partir de polynômes des variables d'origine. Bien sûr, les possibilités ne se liminent pas aux polynômes.
def add_polynomials(X):
#(remove some unwanted warnings)
prev = pd.options.mode.chained_assignment
pd.options.mode.chained_assignment = None
for col in X.columns:
#prise en compte des carrés des variables
X[col+"_square"] = X[col]**2
#prise en compte des racines carrées
X[col+"_sqrt"] = X[col]**(1/2)
pd.options.mode.chained_assignment = prev
return X
featurelist = best_features[:28]
X = ames[featurelist]
X = add_polynomials(X)
Xtrain, Xval, ytrain, yval = X.loc[trainids], X.loc[valids], y.loc[trainids], y.loc[valids]
price_predictor.fit(Xtrain, ytrain)
trainpred, valpred = price_predictor.predict(Xtrain), price_predictor.predict(Xval)
print("MAE (train,val):", MAE(trainpred, ytrain), MAE(valpred, yval))
display_results(yval, valpred, "Régression polynomiale (sqrt et square) - validation set")
Les prédiction semblent suivre une ligne droite. On remarque par contre que certaines prédictions parfois sont très mauvaises, et parfois même négatives (et peuvent aller jusqu'à -5.000.000 selon l'état pseudo-aléatoire utilisé).
Un simple bornage des prédictions entre les prix minimum et maximum du set d'entraînement permet de limiter ces erreurs exceptionnelles.
Néanmoins, selon les objectifs, il convient de définir une métrique d'évaluation appropriée: si quelques rares erreurs très grandes sont tolérables, ou si il est préférable d'avoir des prédictions globalement un peu moins précises mais sans gros écarts. Dans tous les cas, on peut facilement améliorer les résultats en contraignant les prédictions entre des valeurs limites.
lowest_pred_index = yval.index[valpred.argmin()]
display_series(ames.loc[lowest_pred_index])
trainpred = trainpred.clip(ytrain.min(), ytrain.max())
valpred = valpred.clip(ytrain.min(), ytrain.max())
print("MAE (train,val):", MAE(trainpred, ytrain), MAE(valpred, yval))
display_results(yval, valpred, "Régression polynomiale (sqrt et square) corrigée - validation set")
maes2 = pd.DataFrame(columns = ['train','val'])
plot_limit = 0
leg = []
for nsamples in range(10,800,10):
trainset, valset = train_test_split(ames, test_size = 660, random_state=1)
#On extrait un sample sur le training set
trainset_sample = trainset.sample(nsamples, random_state=1)
Xtrain = trainset_sample[featurelist]
ytrain = trainset_sample["SalePrice"]
Xval = valset[featurelist]
yval = valset["SalePrice"]
price_predictor.fit(Xtrain, ytrain)
trainpred = price_predictor.predict(Xtrain)
valpred = price_predictor.predict(Xval)
trainpred = trainpred.clip(ytrain.min(), ytrain.max())
valpred = valpred.clip(ytrain.min(), ytrain.max())
maes2.loc[nsamples] = MAE(ytrain, trainpred), MAE(yval, valpred)
if nsamples in [10,50,100,790]:
display_results(yval, valpred, "Résultats pour nsamples = " + str(nsamples))
plt.title("Résultats pour différentes tailles de set d'entraînement")
maes2.plot()
plt.title("Performances (MAE) pour différentes tailles de set d'entraînement")
plt.xlabel("Nombre d'observations d'entraînement")
plt.ylabel("Mean Absolute Error")
On constate que les performances diminuent (l'erreur augmente) avec le nombre d'échantillons d'entraînement, et que les performances en validation sont au départ beaucoup plus faibles, ce qui est signe d'overfitting. En effet, en utilisant seulement 50 exemples (gauche du graphe), on remarque que le modèle compte presque autant de paramètres (un par variable prédictive). Dans ce cas, les paramètres peuvent être adaptés aux exemples spécifiques d'entraînement pour apprendre leur prédiction presque "par coeur". Le modèle n'est alors pas générique et fonctionne mal sur de nouvelles données (comme le montrent les faibles performances en validation). On constate ensuite que pour 400 échantillons d'entraînement, les performances en entraînement se stabilisent, et les performances en validation sont équivalentes, ce qui signifie que le nombre d'échantillons est suffisant (plus d'overfitting).
len(price_predictor.coef_)
Pour évaluer les performances réelles du modèles, nous devons utiliser un set de données qui n'a pas encore été utilisé, ni pour l'entraînement ni pour la validation. Cela permet de simuler une utilisation réelle du modèle, sur de nouvelles maisons (dont on voudrait estimer le prix). L'évaluation des performances à partir du set de validation seraient biaisée, car il a en effet permis de choisir les paramètres du modèles, en l'occurence le choix des variables prédictives (le modèle a donc été indirectement optimisé pour ces données spécifiques).
Sur la plateforme Kaggle, pour garantir l'absence de biais lors du design du modèle, un set de test est généralement fourni séparément, et les labels ne sont volontairement pas fournis. L'utilisateur de la plateforme doit alors soumettre à la plateforme les prédictions données par le modèle développé, et reçoit le résultat de l'évaluation effectuée par la plateforme. Cela permet notamment de classer différents compétiteurs lors d'un concours: https://www.kaggle.com/competitions
Concernant le dataset utilisé, on peut soumettre des prédictions ici: https://www.kaggle.com/c/home-data-for-ml-course/overview/evaluation
Remarque: lorsque le concours ne sera fini, il sera toujours possible de récupérer les données complète ici: https://www.kaggle.com/prevek18/ames-housing-dataset
# Vérifions déjà les performances sur le set de validation
#Régression linéaire
X = ames[featurelist]
y = ames["SalePrice"]
Xtrain, Xval, ytrain, yval = X.loc[trainids], X.loc[valids], y.loc[trainids], y.loc[valids]
price_predictor.fit(Xtrain, ytrain)
trainpreds = price_predictor.predict(Xtrain).clip(ytrain.min(), ytrain.max())
valpreds_lin = price_predictor.predict(Xval).clip(ytrain.min(), ytrain.max())
MAE(trainpreds,ytrain.values), MAE(valpreds_lin,yval.values)
#Régression polynomiale
X = ames[featurelist]
X = add_polynomials(X)
Xtrain, Xval, ytrain, yval = X.loc[trainids], X.loc[valids], y.loc[trainids], y.loc[valids]
poly_predictor = LinearRegression()
poly_predictor.fit(Xtrain, ytrain)
trainpreds = poly_predictor.predict(Xtrain).clip(ytrain.min(), ytrain.max())
valpreds_poly = poly_predictor.predict(Xval).clip(ytrain.min(), ytrain.max())
MAE(trainpreds,ytrain.values), MAE(valpreds_poly,yval.values)
display_results(yval.values, valpreds_lin, "Régression linéaire")
display_results(yval.values, valpreds_poly, "Régression polynomiale")
Remarquez que vous obtiendrez des résultats différents à chaque fois que vous relancez les cellules ci-dessus avec un random_state différent: cela est dû à la séparation aléatoire des sets d'entraînement et de validation.
Afin d'avoir un indicateur plus robuste du modèle, une possibilité serait de calculer une moyenne de ces résultats, ou d'utiliser un processus de k-fold cross-validation: on découpe l'ensemble des données en k échantillons, et pour chaque échantillon on entraînement un modèle sur les autres données, et on calcule les prédictions pour l'échantillon gardé de côté. On calcule de cette manière des prédictions pour chaque échantillon et on extrait une mesure (e.g., MAE) sur l'ensemble des prédictions.
testset = pd.read_csv("https://raw.githubusercontent.com/titsitits/Python_Data_Science/master/Donn%C3%A9es/test.csv")
ids = testset['Id']
print(len(testset))
#nettoyage du testset
testset = clean_df(testset)
#Colonnes incomplètes
counts = testset.count()
incomplete = counts[counts < len(testset)]
display_series(incomplete.sort_values())
len(incomplete.index)
Il reste deux observations contenant des variables invalides. On va simplement remplacer les valeurs manquantes par les moyennes. Si le nombre d'observations à nettoyer était plus conséquent, il serait pertinent d'analyser les variables à corriger.
for col in incomplete.index:
testset[col] = testset[col].fillna(ames[col].mean())
X = ames[featurelist]
y = ames["SalePrice"]
#On peut ré-entraîner le modèle sur tout le dataset (entraînement+validation), pour le rendre potentiellement plus robuste
price_predictor.fit(X, y)
#select features
Xtest = testset[featurelist]
testpreds_lin = price_predictor.predict(Xtest).clip(ytrain.min(), ytrain.max())
submission = ids.to_frame()
submission["SalePrice"] = testpreds_lin
submission.to_csv("mysubmission_linear_regression.csv", index = False)
Xtest2 = add_polynomials(Xtest)
X = ames[featurelist]
X = add_polynomials(X)
poly_predictor.fit(X, y)
testpreds_poly = poly_predictor.predict(Xtest2).clip(ytrain.min(), ytrain.max())
submission = ids.to_frame()
submission["SalePrice"] = testpreds_poly
submission.to_csv("mysubmission_polynomial_regression.csv", index = False)
#Analysons l'effet de chaque variable catégorielle
all_columns = ames.columns
cat_columns = ames.columns[ames.dtypes == 'object']
diffs_per_group = []
#Bonus: vous pouvez afficher des graphes pour chaque groupe
#ncols = 8
#ncats = len(cat_columns)
#fig, axes = plt.subplots(int(np.ceil(ncats/ncols)),ncols, figsize = (20,30))
i=0
for cat in cat_columns:
group_means = ames.groupby(cat)["SalePrice"].mean()
diffs_per_group.append(group_means.max() - group_means.min())
#Bonus: vous pouvez afficher des graphes pour chaque groupe
#plt.figure(figsize=(10,3))
#group_means.plot.bar()
#ax = axes[int(i/ncols), int(i%ncols)]
#ames.boxplot("SalePrice", by=cat, ax = ax)
i=i+1
best_cats = pd.Series(diffs_per_group, list(cat_columns)).abs().sort_values(ascending = False)
best_cats.plot.bar(figsize=(10,4))
best_cats = best_cats.index.to_list()
#Installer un nouveau package: on appelle une ligne de commande linux grâce au symbole "!". On utilise le gestionnaire de packages python pip pour installer un nouveau package
!pip install catboost
contfeatures = best_features[:28]
catfeatures = best_cats[:-5] #on retire les plus mauvaises
newfeaturelist = contfeatures + catfeatures
X = ames[newfeaturelist]
y = ames["SalePrice"]
Xtrain, Xval, ytrain, yval = X.loc[trainids], X.loc[valids], y.loc[trainids], y.loc[valids]
import numpy as np
from catboost import Pool, CatBoostRegressor
cat_feature_ids = [i for i in range(len(newfeaturelist)) if newfeaturelist[i] in catfeatures]
train_pool = Pool(Xtrain.values, ytrain.values, cat_features=cat_feature_ids)
val_pool = Pool(Xval, yval, cat_features=cat_feature_ids)
all_pool = Pool(X, y, cat_features=cat_feature_ids)
# Spécification des paramètres d'entraînement du modèle
model = CatBoostRegressor(iterations=300,
depth=3,
learning_rate=0.2,
loss_function='RMSE')
#Entraînement du modèle
model.fit(train_pool, silent=True)
# Réalisation de la prédiction en utilisant le modèle obtenu
trainpreds = model.predict(train_pool)
valpreds = model.predict(val_pool)
MAE(trainpreds,ytrain.values), MAE(valpreds,yval.values)
# Spécification des paramètres d'entraînement du modèle
model = CatBoostRegressor(iterations=300,
depth=3,
learning_rate=0.2,
loss_function='RMSE')
#Entraînement du modèle, avec un critère d'arrêt lorsque les performances en validation baissent: éviter le surentraînement (overfitting)
model.fit(train_pool,eval_set = val_pool, early_stopping_rounds = 100, silent = True)
# Réalisation de la prédiction en utilisant le modèle obtenu
trainpreds = model.predict(train_pool)
valpreds_catboost = model.predict(val_pool)
#Pour éviter de prédire des valeurs anormales, on limite les prédictions au range du set d'entraînement
trainpreds = trainpreds.clip(ytrain.min(), ytrain.max())
valpreds = valpreds.clip(ytrain.min(), ytrain.max())
MAE(trainpreds,ytrain.values), MAE(valpreds_catboost,yval.values)
display_results(yval.values, valpreds_lin, "Régression linéaire")
display_results(yval.values, valpreds_poly, "Régression polynomiale")
display_results(yval.values, valpreds_catboost, "CatBoost")
Xtest = testset[newfeaturelist]
test_pool = Pool(Xtest.values, cat_features=cat_feature_ids)
#Pour la soumission, on peut éventuellement réentraîner le modèle sur toutes les données d'entraînement et de validation (pour espérer avoir un modèle plus générique. Cependant: attention à l'overfitting sans critère d'arrêt)
model.fit(all_pool, silent = True)
#Prédictions
testpreds_catboost = model.predict(test_pool).clip(ytrain.min(), ytrain.max())
submission = ids.to_frame()
submission["SalePrice"] = testpreds_catboost
submission.to_csv("mysubmission_catboost.csv", index = False)
#Pour s'assurer de la plausibilité des résultats, on peut comparer leur distribution avec ceux du dataset d'entraînement
plt.subplot(1,2,1)
submission.SalePrice.hist()
plt.subplot(1,2,2)
ames.SalePrice.hist()
Les résultats avec le fichier de baseline (sample_submission.csv) ont été obtenus avec une régression linéaire sur l'année et le mosi de vente, la surface du lot et le nombre de chambres. Les résultats de la soumission sur Kaggle donnent:
MAE = 59346
Les trois soumissions proposées donnent:
Les algotihmes proposés prédisent donc bien mieux mieux les prix des maisons que la baseline.
Le résultat obtenu avec catboost amène à la place 544/13783 au classement de la compétition Kaggle. Il reste donc de la marge de manoeuvre pour améliorer le modèle. (Le meilleur score publique est actuellement à 11824).
Il existe de très nombreuses possibilités pour améliorer le modèle, dont par exemple:
Kaggle est une excellente source d'inspiration pour s'améliorer. On y trouve en effet de très nombreux notebooks (généralement en Python) proposant des idées intéressantes de feature engineering et de développement de modèles prédictifs. Par exemple, le notebook suivant est classé 26ème: https://www.kaggle.com/itslek/stack-blend-lrs-xgb-lgb-house-prices-k-v17 La méthode propose certaines nouvelles variables prédictives simplifiées, ainsi que la moyenne pondérée de différents modèles prédictifs.
w1, w2, w3 = 0.02, 0.08, 0.9
valpreds_blend = w1*valpreds_lin + w2*valpreds_poly + w3*valpreds_catboost
MAE(yval, valpreds_blend)
testpreds_blend = w1*testpreds_lin + w2*testpreds_poly + w3*testpreds_catboost
submission = ids.to_frame()
submission["SalePrice"] = testpreds_blend
submission.to_csv("mysubmission_blend.csv", index = False)
MAE = 14273
Nouveau classement: 516/13783