Emissions de CO2 et de polluants des véhicules commercialisés en France¶
Source : data.gouv.fr¶
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
Définitions des colonnes du jeu de données¶
col_def = pd.read_excel("carlab-annuaire-variable.xlsx")
col_def
Chargement du jeu de données¶
df = pd.read_csv("mars-2014-complete.csv",encoding="latin",sep=";",decimal=",")
df = df.loc[:, ~df.columns.str.contains('^Unnamed')]
df.head()
Remarque : la colonne "lib_mrq" semble correspondre à la colonne définie comme "lib_marq_utac" (cf. définition des colonnes du jeu de données).
Résumé des données¶
df.describe()
df.describe(include=['O'])
Suppression de certaines colonnes¶
On décide de supprimer certaines colonnes du jeu de données dans la mesure où on estime qu'il n'est pas nécessaire d'en tenir compte dans le cadre de notre analyse. On supprime les colonnes correspondantes aux modèle du dossier, le modèle commercial, la désignation commerciale, le CNIT et le tvv(type variante version) car le numéro de ligne suffie à identifier un véhicule de manière unique et la marque permet de lui associé une "classe" d'appartenance. D'autre part, quand bien même ces colonnes apportent des informations supplémentaires (sur le nombre de chevaux par exemple), les colonnes numériques caractérisent déjà ces éléments (par exemple les colonnes relatives à la puissance du véhicule et aux résultats des tests). De même, on ne considérera pas la date de mise à jour (86% d'entre elles sont le 14 mars). On conserve la colonne du champ V9 car elle donne une information relative à la norme euro, elle-même reliée aux émissions de polluants (dont le CO2)
df = df.drop(['lib_mod_doss','lib_mod','dscom','cnit','tvv','date_maj'],axis=1)
df.head()
Nouvelle description¶
Le jeu de données à changé, on en propose donc une nouvelle description
df.describe()#description des variables numériques
df.describe(include=['O'])#description des variables catégorielles
Traitement des données manquantes (nan)¶
Il manque des données dans le dataset, ou elles ne sont pas "correctes" (pour diverses raisons). Quoiqu'il en soit, ce sont des éléments que l'on ne peut rentabiliser numériquement. Soit on en fait une estimation par diverses approches ent utilisant la méthode fillna. Ici, nous choisirons de simplement supprimer les lignes (et donc de perdre des données) contenant des valeurs manquantes pour ne pas "changer la réalité" : dans la mesure où chaque véhicule est différent faire une interpolation où propager une valeur d'un véhicule à l'autre ne semble pas approprié.
Avant tout, il est important de noter le nombre de valeurs manquantes, sur le total des valeurs du dataset.
val_manq = {}
val_manq["Valeurs manquantes"] = df.isna().sum().sum()
val_manq["Total de valeurs"] = df.shape[0]*df.shape[1]
vm = pd.DataFrame(val_manq,index=[""]).transpose()
vm
D'autre part, il est important d'observer le nombre de valeurs manquantes sur chacune des colonnes du dataframe:
df.isna().sum().plot.bar(figsize=(20,3),title="Répartition du nombre de valeur manquantes sur chaque colonne du dataset",rot=20).set_ylabel("Occurence")
Il y a 45271 valeurs manquantes sur la colonne 'hc' (les résultats d’essai HC en gramme par km) sur 55044 lignes dans le jeu de données. Autrement dit, si on supprime les lignes contenant les valeurs manquantes, notre jeu de données serait considérablement impacté en termes de pertes de données. Pour éviter ce coup de bélier, on supprime cette colonne du dataset avant de supprimer les lignes contenant des valeurs manquantes
df=df.drop(['hc'],axis=1)
df.dropna(inplace=True)
Le jeu de données étant modifié, il convient d'observer sa nouvelle description : À présent, il n'y a plus de valeurs manquantes :
df.describe()
df.describe(include=['O'])
À présent, il n'y a plus de valeurs manquantes :
pd.DataFrame(df.isna().sum()).transpose()
Visualisations¶
def visu(df,categorielles=False,numeriques=False,horizontal=True):
if numeriques==True:
for elt in df.columns:
if pd.api.types.is_numeric_dtype(df[elt]):
plt.figure()
df[elt].plot.hist(bins=50,density=True,
title="Distribution de la variable "+"'"+str(elt)+"'")
df[elt].plot.density().set_ylabel("Densité de probabilité / Fréquence")
plt.xlabel("Valeurs de la variable '"+elt+"' "+"("+str(col_def["unité"][col_def["nom-colonne"] == elt].values[0])+")"+"\n *"+str(col_def["légende"][col_def["nom-colonne"] == elt].values[0]))
if categorielles==True:
for elt in df.columns:
if pd.api.types.is_object_dtype(df[elt]):
if not(horizontal):
plt.figure()
df[elt].value_counts(sort=True).plot.bar(title="Distribution de la variable "+"'"+str(elt)+"'",log=True,zorder=3).grid()
else:
plt.figure()
df[elt].value_counts(sort=True,ascending=True).plot.barh(title="Distribution de la variable "+"'"+str(elt)+"'",log=True,zorder=3).grid()
plt.xlabel("Occurence (échelle logarithmique)")
visu(df,categorielles=True,numeriques=True,horizontal=False)
Concernant les variables quantitatives, on observe des similarités dans l'allure des distributions des variables 'conso_urb', 'conso_exurb', 'conso_mixte'et 'co2'. De même, il semble y avoir une ressemblance entre les distributions des variables 'nox' et 'hcnox'. Aussi, les distributions il semble y avoir un lien entre les distributions des variables 'puiss_admin_98' et 'puiss_max'. On pourrait faire la conjecture que ces variables sont corrélées. Cette conjecture est infirmée par la valeur des coefficients de corrélation de Pearson de la section Analyse des corrélations.
Par ailleurs, la marque Mercedes semble surreprésentée dans le jeu de données. Observons les distributions des variables pour cette marque:
for elt in df.columns:
if elt not in ['lib_mrq','cod_cbr','hybride','typ_boite_nb_rapp','champ_v9','Carrosserie','gamme']:
plt.figure()
df[df["lib_mrq"]=="MERCEDES"][elt].plot.hist(bins=50,density=True,
title="Distribution de la variable "+"'"+str(elt)+"' pour la marque Mercedes")
df[elt].plot.density().set_ylabel("Densité de probabilité / Fréquence")
plt.xlabel("Valeurs de la variable '"+elt+"' "+"("+str(col_def["unité"][col_def["nom-colonne"] == elt].values[0])+")"+"\n *"+str(col_def["légende"][col_def["nom-colonne"] == elt].values[0]))
for elt in ['cod_cbr','hybride','typ_boite_nb_rapp','champ_v9','Carrosserie','gamme']:
plt.figure()
df[df["lib_mrq"]=="MERCEDES"][elt].value_counts().plot.bar(title="Répartition des voitures Mercedes dans les catégories de la variable '"+str(elt)+"'",log=True,zorder=3).grid()
Par comparaison avec les répartitions dans les variables catégorielles dans tout le jeu de donné, on observe qu'il y a des similarités entre la répartition des Mercedes dans ces mêmes variables catégorielles. Or les voitures de la marque Mercedes sont surreprésentées. Le dataset comporte un biais de représentation. Il faudra en tenir compte dans la conception de l'apprentissage automatique.
Analyse des corrélations¶
df.corr(method="pearson").style.applymap(lambda x: "background-color: red" if round(abs(x),1) >= 0.7 else "background-color: white")
La matrice de corrélation nous permet de conclure quant aux potentielles connexions*** entre les variables quantitatives. On observe par exemple qu'il semble y avoir une forte corrélation entre les consommations urbaine, extra urbaine, mixte et les émissions de CO2. De même, il semblerait y avoir un lien entre les résultats d'essais NOx et HC+NOx. Aussi, dans une mesure plus faible, il y aurait corrélation entre es résultats d'essais NOx et les différents types de consommations et les émissions de CO2. Enfin, la masse en ordre de marche mini semblerait relativement "faiblement" corrélée avec les diffénts types de consommations, les émissions de CO2 tandis que la masse en ordre de marche maxi serait liée à la consommation extra urbaine.
L'analyse des corrélations pourrait présenter un intérêt par la suite, par exemple pour diminuer le nombre de variables, si nécessaire (par exemple pour rendre le jeu de données moins complexe), en perdant le moins d'information possible.
***Dans l'hypothèse où l'on considérerait que la corrélation (calculée statistiquement) serait significative si le coefficient de Pearson est supérieur ou égal à 0.7 en valeur absolue (la corrélation pourrait être "significative" pour des valeurs absolues plus faibles mais nous ne ferons pas d'hypothèse sur celles-ci)
En guise de résumé graphique, nous pouvons observer les nuages de points résultant du tracé des données d'une colonne en fonction de celles d'une autre colonne du dataset. C'est une façon d'apprécier visuellement le potentiel lien entre deux types de données quantitatives.
Tracés des relations par paires¶
import matplotlib as mpl
mpl.rcParams["axes.labelsize"] = 18
g = sns.PairGrid(df)
g.map_diag(sns.histplot)
g.map_offdiag(sns.scatterplot)
g.add_legend()
Remarque : Il semble y avoir un lien entre les masses en ordre de marche mini et maxi que l'on ne constate pas dans la matrice de corrélation (coefficient de pearson de l'ordre de 0.64). Cela pourrait s'expliquer par l'ensemble de points qui n'est pas sur la "diagonale" apparente, formée par les autres points du graphique. Ces derniers étant alignés horizontalement/verticalement (signe de l'indépendance de deux variables) seraient peut-être la cause d'un coefficient de pearson relativement plus faible rendant plus difficile la conclusion quant à la corrélation des variables.¶
Machine learning¶
On cherche à établir un lien entre les caractéristiques des véhicules et leurs émissions de CO2. On va donc dans les parties qui vont suivre, tester des approches machine learning (régression multilinéaire, réseaux de neurones) pour essayer d'établir ce lien. Une analyse de la performance de ces approches concluera cette étude¶
Préparation des données¶
On va préparer les données (i.e définir les données d'entrées et de sortie de notre modèle)
Encodage des variables qualitatives¶
On fait le choix d'utiliser un encodage one hot : chaque catégorie devient une colonne du tableau et si un véhicule appartient à celle-ci, la valeur 1 est insérée dans la cellule de cette colonne, et 0 sinon. Cela a pour effet d'augmenter le nombre de colonnes du tableau
df_enc = df.copy()
df_enc=pd.get_dummies(df_enc, columns=["lib_mrq","cod_cbr","hybride","typ_boite_nb_rapp","champ_v9","Carrosserie","gamme"])
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(df_enc.drop(['co2','co_typ_1','nox','hcnox','ptcl'],axis=1), df_enc[['co2','co_typ_1','nox','hcnox','ptcl']], train_size=0.8, random_state=42)
Modèles linéaires¶
Dans un premier temps, on propose d'utiliser un modèle linéaire multivarié. On fait l'hypothèse que chaque variable explicative contribue linéairement à la valeur des variables à expliquer : les émissions de CO2(co2), de monoxyde de carbone(co_typ_1) d'hydrocarbures imbrûlés(hc), de NOx et de particules (ptcl); toutes en grammes par kilomètres.
import statsmodels.api as sm
from sklearn.model_selection import train_test_split
# On définit un modèle linéaire par variable à expliquer. Une classe est définie à cette fin.
class multi_lm:
def __init__(self,dependent,independent,data_encoded,train_size=0.8): #Attention à encoder les variables catégorielles
self.dependent = dependent #liste des noms des colonnes des variables à expliquer
self.independent = independent #liste des noms des colonnes des variables explicatives
self.df = data_encoded
self.train_size = train_size
self.X_train,self.X_test,self.y_train,self.y_test = self.build_sets()
self.anovas = self.build_anovas() #liste de modèles linéaires
self.shape_tot = self.df.shape
self.r2 = self.r2()
def build_sets(self): #Subdivision du dataset en données d'entraînement, test et validation
X_train, X_test, y_train, y_test = train_test_split(self.df.drop(self.dependent,axis=1), self.df[self.dependent], train_size=self.train_size, random_state=42)
X_train = sm.add_constant(X_train)
X_test = sm.add_constant(X_test)
return X_train, X_test, y_train, y_test
def build_anovas(self): #Construction des modèles linéaires
l = []
for elt in self.dependent:
l.append(sm.OLS(self.y_train[[elt]], self.X_train).fit())
return l
def predict(self,x): # Prédiction sur un jeu de données nouveau (ou non)
l = []
x = sm.add_constant(x)
for elt in self.anovas:
l.append(elt.predict(x))
return pd.DataFrame(l,index=self.y_train.columns).transpose()
def resumes(self): #Résumés
print("Résumés des modèles linéaires \n")
for elt in self.anovas:
print("\n")
print(elt.summary())
print("\n")
def r2(self):
keys = self.y_train.columns
dico = {}
for elt,key in list(zip(self.anovas,keys)):
dico[key] = elt.rsquared
return pd.DataFrame(dico,index=[r'$R^2$'])
def coefs(self): #coefficients des variables explicatives
keys = self.y_train.columns
dico = {}
for elt,key in list(zip(self.anovas,keys)):
dico[key] = elt.params
return pd.DataFrame(dico)
def pvals(self):
keys = self.y_train.columns
dico = {}
for elt,key in list(zip(self.anovas,keys)):
dico[key] = elt.pvalues
return pd.DataFrame(dico)
def evaluate(self):
Y = self.predict(self.X_test)
YY = (Y - Y.mean()) / Y.std()
y = self.y_test
yy = (y - y.mean())/y.std()
return pd.DataFrame((((YY - yy)**2).mean()),columns=["mse"]).transpose()
multi_model = multi_lm(train_size=0.99,dependent=['co2','co_typ_1','nox','hcnox','ptcl'],independent=list(set(list(df_enc.columns)).difference(list(['co2','co_typ_1','nox','hcnox','ptcl']))),data_encoded=df_enc)
multi_model.resumes()
Remarques:
- On sait, d'après le tableau des corrélations, que plusieurs variables explicatives sont fortement corrélées entre elles(e.g les différents types de consommation, la puissance administrative et la puissance fiscale). Le modèle serait donc différent si l'on ne considérait qu'une seule des variables corrélées .
- Les modèles linéaires sont conçus sur l'ensemble du jeu de données mais on laisse la possibilité de les construire sur la base d'un jeu de données d'entraînement et de le tester sur le jeu de données de test (il suffit de changer le paramètre 'train_size' lors de l'instanciation des modèles)
- Parfois les p-values associées aux coefficient sont "conventionnellement" élevées (i.e supérieures à 0.05), signe que l'on prendrait "trop de risque" à considérer que les valeurs constatées (par le calcul) ne sont pas dues au hasard.
Performance du modèle linéaire¶
Les performances des modèles linéaires semblent dépendre de la variable que l'on veut expliquer. À ce propos, voici, les valeurs des R$^2$ pour chacun des modèles linéaires:
r2s = multi_model.r2
r2s
Nous pouvons noter que la qualité de la régression linéaire est meilleure lorsque la variable de sortie est l'émission de CO$_2$ (R$^2$ $\approx$ 0.999) que lorsqu'il s'agit des émissions de monoxyde de carbone (R$^2$ $\approx$ 0.63) ou de particules(R$^2$ $\approx$ 0.37). Or, nous pouvons remarquer dans le tableau des corrélations que la variable 'co2' est relativement fortement corrélée à plus certaines variables (relatives aux types de consommations) que la variable 'co_typ_1'(monoxyde de carbone) ou 'ptcl'(particules) dont la corrélation avec les autres variables semble très faible. Il serait donc cohérent, dans le cadre de l'utilisation d'un modèle linéaire, qu'il soit plus difficile de modéliser le liens entre des variables dont les corrélations sont plus faibles relativement. C'est d'ailleurs ce qui expliquerait ces difficultées.
Pour pallier ce problème, deux éléments semblent envisageable :
- Soit la corrélation, n'étant pas linéaire, pourrait être d'un autre type (quadratique, logarithmique, cubique...)
- Soit on teste une autre approche pour la prédiction des variables
La première approche semble compliquée : si on observe les Tracés des relations par paires, il semble difficile d'observer un type particulier de relations des variables 'co_typ_1' et 'ptcl' avec les variables restantes.
On propose donc d'aborder d'autres aspects.
Une autre évaluation¶
Nous avons déjà de bonnes indications concernant les performances du modèles.
Nous proposons d'ajouter une phase d'évaluation supplémentaire en testant le modèle (conçu sur la base des données d'entraînement) sur les données de test. Nous utiliserons l'erreur quadratique moyenne (mean squared error ou mse*) comme métrique. Nous diviserons donc le jeu de données initial en un jeu de données d'entraînement (80%) et un jeu de données de test(20%). Et les données (de test et prédites) seront standardisées (de sorte que les valeurs de chaque quantité soient échelonnées sur des plages de valeurs plus proches) en opérant la transformation suivante : $$y_i^{nouveau} = \frac{y_i - \bar{y_i}}{\sigma(y_i)}$$
*$mse = \frac{1}{N}\sum_{i=1}^{N}\left(\mathbf{y_i^{prédit}} - \mathbf{y_i^{réel}}\right)^{2}$
- $N$ : Taille de l'échantillon de test
- La mse est calculée en considérant les quantités standardisées (i.e. $y_i^{nouveau}$) pour les données prédites et de test
- Plus la mse est proche de zéro et plus cela signifie que la valeur prédite est proche de la vérité terrain
Voici les résultats que l'on obtient :
import warnings
warnings.filterwarnings('ignore')
multi_model_test = multi_lm(train_size=0.8,dependent=['co2','co_typ_1','nox','hcnox','ptcl'],
independent=list(set(list(df_enc.columns)).difference(list(['co2','co_typ_1','nox','hcnox','ptcl']))),
data_encoded=df_enc)
multi_model_test.evaluate()
Les tendances observées sur ce tableau semblent similaires à celles constatées sur le tableau des coefficients de pearson. Le modèle linéaire semble meilleur dans la prédiction des émissions de CO$_2$, moins bon dans la prédiction des résultats d'essais NO$_x$ et HC+NO$_x$. Enfin, il semble beaucoup moins performant dans la prédiction des résultats d'essais CO type 1 et d'émissions de particules.
Modélisation par réseaux de neurones¶
l'approche par modèles linéaires semblait performante dans la "prédiction" des émissions de CO$_2$, et, dans une mesure plus faible, des émissions de $NO_x$ et $HC+NO_x$ (hydrocarbures imbrûlés et $NO_x$). Par ailleurs, les modèles linéaires ne semblaient pas pertinent dans la "prédiction" de monoxyde de carbone et de particules. Ceci serait lié à une corrélation relativement faible aux variables explicatives.
Toujours dans l'optique d'évaluer des méthodes d'apprentissage automatique, on propose de tenter une approche par réseaux de neurones pour "prédire" les émissions de CO2 et polluants (particules, monoxyde de carbone, $NO_X$ et $HC+NO_x$).
X_train, X_test, y_train, y_test = train_test_split(df_enc.drop(['co2','co_typ_1','nox','hcnox','ptcl'],axis=1),
df_enc[['co2','co_typ_1','nox','hcnox','ptcl']], train_size=0.85)
#Pour générer un jeu de validation
X_train, X_valid, y_train, y_valid = train_test_split(X_train,
y_train, train_size=0.9)
# Création d'un dataframe pour récapituler les tailles de chaque jeux de données et leurs proportions respectives
dico={}
dico["Entraînement"] = y_train.shape[0]
dico["Validation"] = y_valid.shape[0]
dico["Test"] = y_test.shape[0]
data_nb = pd.DataFrame.from_dict(dico,orient='index',columns=['Nombre de données'])
data_nb['Proportions (%)'] = data_nb["Nombre de données"]*100 / data_nb["Nombre de données"].sum()
Le dataset est subdivisé en un jeu de données d'entraînement, de validation (pour garantir la qualité de l'apprentissage) et de test. En voici un récapitulatif :
data_nb
Environ 76% des données seront utilisées pour entraîner le modèle neuronal tandis que 9% serviront à observer la "qualité" de l'apprentissage et la capacité du modèle à généraliser. Enfin, on testera les performances du modèle sur les 15% des données restantes.
import tensorflow as tf
from tensorflow import keras
from tensorflow.keras import layers
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Dense
from tensorflow.keras.layers import Flatten
from tensorflow.keras import regularizers
model = keras.Sequential()
model.add(Flatten(input_shape=(X_train.shape[1],))) #entrées
model.add(Dense(X_train.shape[1],activation='relu',activity_regularizer=regularizers.L2(1e-03)))
model.add(Dense(100,activation='relu',activity_regularizer=regularizers.L2(1e-03)))
model.add(Dense(100,activation='relu',activity_regularizer=regularizers.L2(1e-03)))
model.add(Dense(50,activation='relu',activity_regularizer=regularizers.L2(1e-03)))
model.add(Dense(20,activation='relu',activity_regularizer=regularizers.L2(1e-03)))
model.add(Dense(y_train.shape[1],activation='linear'))#sorties
On choisit un réseau dense à 6 couches (sans compter la couche d'entrée) avec paramètre de régularisation de type l2 (nous pourrions ajouter un "dropout" pour éviter le surapprentissage mais choisissons de tester cette configuration dans un premier temps). En voici un résumé:
model.summary()
Cette architecture peut, bien entendu, être amenée à évoluer pour améliorer les performances du modèle et étudier d'autres aspects de l'approche par réseaux neuronaux mais l'objectif ici est de tester une première approche et de la comparer au modèle linéaire.
from tensorflow.python.ops import math_ops
from tensorflow.keras import backend as K
def rmse(y_true,y_pred):
return K.sqrt(K.mean(K.square(y_pred-y_true)))
model.compile(loss='mse',
optimizer = keras.optimizers.Adam(lr=0.0001),metrics=['mse','mae',tf.keras.metrics.RootMeanSquaredError()])
history = model.fit(X_train, y_train, batch_size=64, epochs=250,validation_data=(X_valid, y_valid))
model.save('model.h5')
plt.plot(history.history['loss'],label="Courbe d'apprentissage")
plt.plot(history.history['val_loss'],label="Loss sur les données de validation")
plt.title("Évolution de la fonction de perte (loss) durant l'entraînement")
plt.xlabel("Itération")
plt.ylabel("Valeur de la fonction de perte")
plt.legend()
plt.yscale('log')
def evaluate(x,y):
Y = pd.DataFrame(model.predict(x),columns=["co2","co_typ_1","nox","hcnox","ptcl"],index=y_test.index)
YY = (Y - Y.mean()) / Y.std()
yy = (y - y.mean())/y.std()
return pd.DataFrame((((YY - yy)**2).mean()),columns=["mse"]).transpose()
Performance du réseau de neurones¶
Nous proposons ici d'utiliser la même méthode d'évaluation que celle du modèle linéaire afin de comparer les deux méthodes. Voici donc un tableau récapitulatif de ce que l'on obtient:
evaluate(X_test,y_test)
Le réseau de neurones serait donc plus performant dans la prédiction des émissions de co2 que dans la prédiction des essais nox, hcnox et ptcl. On remarquera que le réseau de neurones est par ailleurs meilleur dans la prédiction des émissions de monoxyde de carbone (co_typ_1) que dans celles des nox et hcnox malgré une corrélation plus forte de ces variables avec les variables explicatives. Par ailleurs Les émissions de particules restent difficile à prédire avec ce dataset dans la mesure où elles sont relativement faiblement liées aux variables explicatives.
Plus généralement, le réseau de neurones actuellement conçu est moins performant que les modèles linéaires (cf. tableau de la section "Une autre évaluation"). Une possibilité de développement de cette activité serait donc de chercher si le modèle neuronal ne pourrait pas être amélioré (voir si des meilleurs hyperparamètres sont possibles ou si la méthodologie pourrait être améliorée; tester des modèles plus complexes ou d'autres approches...).