Important¶
Le code de ce notebook est montré par défaut pour une compréhension complète.
Si vous souhaitez parcourir ce rapport pour information sans porter attention aux détails techniques, vous pouvez cacher le code en cliquant sur le bouton ci-dessous.
from IPython.display import HTML
HTML('''<script>
code_show=true;
function code_toggle() {
if (code_show){
$('div.input').show();
} else {
$('div.input').hide();
}
code_show = !code_show
}
$( document ).ready(code_toggle);
</script>
<b><a href="javascript:code_toggle()">Cacher le code</a></b>.''')
Introduction¶
Ce notebook propose une analyse des accidents survenus dans le département de l'Essonne entre début 2015 et fin 2016.
Les informations à notre disposition sont les suivantes:
adr: Adresse postale : variable renseignée pour les accidents survenus en agglomérationagg: 1 – Hors agglomération, 2 – En agglomérationan: année de l'accidentatm: conditions atmosphériques : 1. Normale, 2. Pluie légère, 3. Pluie forte, 4. Neige - grêle, 5.Brouillard - fumée, 6. Vent fort - tempête, 7. Temps éblouissant, 8. Temps couvert, 9. Autrejour: Jour de l'accidentlat: latitudelong: longitudelum: conditions d’éclairage dans lesquelles l'accident s'est produit (1. Plein jour, 2. Crépuscule ou aube,3. Nuit sans éclairage public, 4. Nuit avec éclairage public non allumé, 5. Nuit avec éclairage public allumé) \n",mois: Mois de l'accidentNum_Acc: Numéro d'identifiant de l’accidenthrmn: heure et minutes (HHMM)gravite: gravité de l'accident
Notre client est particulièrement intéressé par les informations suivantes:
- combien y-a-t-il eu d’accidents en agglomération, par année ?
- quel mois est généralement le plus mortel ?
- quel mois est généralement le moins mortel ?
Nous répondrons à ces questions, avant d'appronfondir notre interrogation du dataset et de proposer un modèle.
Sommaire¶
Analyse¶
Exploration¶
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
df = pd.read_csv('dataset/accidents_essonne_2015_2016.csv')
df.info()
df.describe()
Le dataset comprend 1765 accidents uniques. Aucune valeur manquante ne semble à déplorer, et les différentes features ont déjà été converties dans le type idéal. Les valeurs minimales et maximales des différentes features concordent avec les indications du README.
agg étant une fonction pandas, il est possible que la manipulation de la colonne agg (indiquant si l'accident a eu lieu en agglomération ou hors agglomération) pose quelques problèmes syntactiques. Nous allons donc la renommer en agglo.
df = df.rename(columns={'agg':'agglo'})
df.columns
Il semble que le dataset ait déjà été nettoyé. Nous allons tout de même imprimer les valeurs uniques de différentes colonnes afin de nous assurer qu'aucune valeur étrange n'apparait dans le dataset.
unique_years = df['an'].unique()
np.sort(unique_years)
unique_months = df['mois'].unique()
np.sort(unique_months)
unique_days = df['jour'].unique()
np.sort(unique_days)
unique_lum = df['lum'].unique()
np.sort(unique_lum)
unique_agg = df['agglo'].unique()
np.sort(unique_agg)
unique_atm = df['atm'].unique()
np.sort(unique_atm)
Tout semble correct. Dans l'attente d'une confirmation de la part du client ou du collaborateur nous ayant fourni le dataset, nous procèderons à l'analyse directement.
Nous allons imprimer les 5 premiers et 5 derniers rangs du dataset pour référence, puis nous commencerons à répondre aux questions qui nous ont été posées.
df.head()
df.tail()
Il semble que la colonne hrmn (heure et minutes) apparaisse en style militaire (en milliers). Nous avons pris la décision de la laisser ainsi: elle demeure parfaitement compréhensible et exploitable. Si une conversion dans un format plus approprié s'avère nécessaire au cours de l'analyse, nous y procèderons à ce moment.
Questions¶
Combien y-a-t-il eu d’accidents en agglomération, par année ?¶
df[(df.an == 15) & (df.agglo == 2)].shape
264 accidents ont eu lieu en agglomération en 2015.
df[(df.an == 16) & (df.agglo == 2)].shape
598 accidents ont eu lieu en agglomération en 2016 (soit une augmentation de 126%).
Mortalité selon le mois¶
deathly = df[df.gravite == 'Mortel']
print(deathly.shape)
deathly = pd.DataFrame(deathly.groupby('mois')['gravite'].count())
deathly
There were 43 deathly accidents in 2015 and 2016.
months = ['Janvier', 'Février', 'Mars', 'Avril', 'Mai', 'Juin',
'Juillet', 'Août', 'Septembre', 'Octobre', 'Novembre', 'Décembre']
g = sns.barplot(deathly.index, deathly.gravite, palette="Blues")
g.set_xticklabels(months, rotation=45)
plt.title('Mortalité selon le mois')
plt.ylabel('nombre d\'accidents mortels')
plt.show()
Les données montrent que le mois le plus mortel dans l'Essonne est le mois de Décembre, suivi par le mois de Septembre.
À l'inverse, les mois de Janvier, Mars et Juin sont les moins mortels.
Il serait également intéressant d'observer l'impact des conditions atmosphériques et des conditions d'éclairage sur la mortalité.
Approfondissement¶
Combien y a-t-il eu d'accidents en 2015 et 2016?¶
df.groupby('an').size()
Le dataset compte 393 accidents en 2015 et 1372 en 2016. Nous nous gardons de conclure que le nombre d'accidents a plus que triplé en un an, et vérifirons que nous avons a disposition l'intégralité des données concernant 2015 et 2016 lors de notre prochain entretien avec le client.
Impact des conditions atmosphériques¶
deathly_atm = df[df.gravite == 'Mortel']
print(deathly_atm.shape, '\n')
deathly_atm = pd.DataFrame(deathly_atm.groupby('atm')['gravite'].count())
deathly_atm
conds = ['Normale', 'Pluie légère', 'Pluie forte', 'Brouillard — fumée', 'Temps éblouissant', 'Temps couvert']
g = sns.barplot(deathly_atm.index, deathly_atm.gravite, palette="Reds_d")
g.set_xticklabels(conds, rotation=45)
plt.title('Mortalité selon les conditions atmosphériques')
plt.ylabel('nombre d\'accidents mortels')
plt.xlabel('conditions atmosphériques')
plt.show()
Étonamment, il semble que l'écrasante majorité des accidents mortels se produit lorsque les conditions atmosphériques sont normales. D'autres facteurs peuvent entrer en jeu, comme la zone où l'accident se produit, ou la vitesse, qui n'apparaît pas dans ce dataset.
Toutefois, les conditions atmosphériques ont-elles une distribution différente, si l'on prend en compte tous les accidents?
Cette fois, nous allons regarder les proportions.
atm = pd.DataFrame(df.groupby('atm')['gravite'].count())
atm['gravite'] = atm['gravite'] / atm['gravite'].sum()
conds = ['Normale', 'Pluie légère', 'Pluie forte', 'Neige — grêle', 'Brouillard — fumée',
'Vent fort — tempête', 'Temps éblouissant', 'Temps couvert', 'Autre']
g = sns.barplot(atm.index, atm.gravite, palette="Oranges_d")
g.set_xticklabels(conds, rotation=45)
plt.title('Accidents selon les conditions atmosphériques')
plt.ylabel('proportion d\'accidents mortels')
plt.xlabel('conditions atmosphériques')
plt.show()
La distribution générale des accidents en fonction des conditions atmosphériques emble similaire: la majorité des accidents se produit sous des conditions atmosphériques normales.
Impact des conditions d'éclairage¶
deathly_lum = df[df.gravite == 'Mortel']
print(deathly_lum.shape, '\n')
deathly_lum = pd.DataFrame(deathly_lum.groupby('lum')['gravite'].count())
lum_conds = ['Plein jour', 'Crépuscule ou aube', 'Nuit sans éclairage public', 'Nuit avec éclairage public allumé']
g = sns.barplot(deathly_lum.index, deathly_lum.gravite, palette="Purples_d")
g.set_xticklabels(lum_conds, rotation=45)
plt.title('Accidents selon la luminosité')
plt.ylabel('nombre d\'accidents mortels')
plt.xlabel('conditions d\'éclairage')
plt.show()
La majorité des accidents mortels se produit de jour, mais un tiers d'entre eux se produit de nuit, sans éclairage public.
Cette distribution est-elle cohérente avec la totalité des accidents?
lum = pd.DataFrame(df.groupby('lum')['gravite'].count())
lum['gravite'] = lum['gravite'] / lum['gravite'].sum()
lum
lum_conds = ['Plein jour', 'Crépuscule ou aube', 'Nuit sans éclairage public', 'Nuit avec éclairage public non allumé', 'Nuit avec éclairage public allumé']
g = sns.barplot(lum.index, lum.gravite, palette="Purples_d")
g.set_xticklabels(lum_conds, rotation=30)
plt.title('Accidents selon la luminosité')
plt.ylabel('proportion d\'accidents mortels')
plt.xlabel('conditions d\'éclairage')
plt.show()
La majorité des accidents, tous types confondus, se produit de jour, et légèrement plus de 10% d'entre eux ont lieu la nuit sans éclairage public. Il semble donc que les conditions d'éclairage aient un impact sur la mortalité d'un accident; cependant, des test plus avancés sont nécessaires pour confirmer et quantifier cet impact.
Impact de l'horaire¶
mortel = df[df.gravite == 'Mortel']
grave = df[df.gravite == 'Grave']
leger = df[df.gravite == 'Leger/Indemnes']
mortel = df.loc[(df['gravite']== 'Mortel')]['hrmn']
grave = df.loc[(df['gravite']== 'Grave')]['hrmn']
leger = df.loc[(df['gravite']== 'Leger/Indemnes')]['hrmn']
hours = ('0:00', '1:00', '2:00', '3:00', '4:00', '5:00',
'6:00', '7:00', '8:00', '9:00', '10:00', '11:00',
'12:00', '13:00', '14:00', '15:00', '16:00', '17:00',
'18:00', '19:00', '20:00', '21:00', '22:00', '23:00')
bins = np.arange(0, 2400, 100)
g1 = sns.distplot(mortel, bins=bins, color='#ea0012', hist=False, label='mortel')
g2 = sns.distplot(grave, bins=bins, color='#ffb200', hist=False, label='grave')
g3 = sns.distplot(leger, bins=bins, color='#00af46', hist=False, label = 'léger / indemne')
plt.xlim(0, 2400)
plt.locator_params(axis='x', nbins=24)
g3.set_xticklabels(hours, rotation=45)
plt.title('Distribution des accidents en fonction de l\'heure')
plt.xlabel('heure')
plt.ylabel('fréquence')
plt.legend()
plt.show()
En observant ce graphique, on distingue une distribution binomiale pour les acidents légers et graves, et plus "normale" pour les accidents mortels.
L'on s'aperçoit que les horaires les plus accidentogènes sont le matin antre 7:00 et 9:00, et l'après-midi entre 16:00 et 19:00. Ces horaires correspondent aux heures de début et de fin des journées de travail traditionnelles: la circulation est plus dense, et le risque d'accidents augmente.
L'on peut constater que les courbes des accidents graves et des accidents légers se suivent de près: il serait intéressant de mener une étude plus approfondie afin de voir quels critères interviennent, et font basculer un accident potentiellement peu grave dans la catégorie des accidents graves.
Les accidents mortels ont tendance à survenir plus souvent l'après-midi. Cependant, le dataset ne comprend que 43 accidents mortels, sur une totalité de 1765 accidents. La distribution des accidents mortels est donc moins fiable que celle des accidents non-mortels, et il serait judicieux d'obtenir plus de données les concernant, par exemple en croisant ce dataset avec ceux d'autres départements limitrophes.
Visualisation¶
Puisque nous étudions des données géographiques, il va de soi de les visualiser sur une carte.
Nous allons d'abord les visualiser sous forme de nuage de points.
import bokeh.plotting as bp
bp.output_notebook()
import bokeh.charts
p = bokeh.charts.Scatter(df, x='long', y='lat', color='gravite',
tools="crosshair, hover, wheel_zoom, pan")
bp.show(p)
L'on peut distinguer 4 outliers ici: l'un au sud (entre Toulouse et Montpellier), deux à l'ouest (entre Strasbourg et Ulm, en Allemagne), et un au nord (proche de Bruxelles, au sud de Leuven), tous sans gravité, et tous éloignés du nuage de points principal (ces locations peuvent être vérifiées sur la carte ci-dessous).
Ces 4 points devraient être mentionnées lors de la prochaine réunion avec le client. Il peut s'agir de données exactes, ou d'erreurs; quoi qu'il en soit, il est important de comprendre comment ces données ont été générées et à quoi elles correspondent.
lat_mortel = df.loc[(df['gravite']== 'Mortel')]['lat']
lat_grave = df.loc[(df['gravite']== 'Grave')]['lat']
lat_leger = df.loc[(df['gravite']== 'Leger/Indemnes')]['lat']
long_mortel = df.loc[(df['gravite']== 'Mortel')]['long']
long_grave = df.loc[(df['gravite']== 'Grave')]['long']
long_leger = df.loc[(df['gravite']== 'Leger/Indemnes')]['long']
an_mortel = df.loc[(df['gravite']== 'Mortel')]['an']
import warnings
warnings.filterwarnings('ignore')
from bokeh.io import curdoc, output_file, show
from bokeh.layouts import row
from bokeh.models import ColumnDataSource, GMapOptions, HoverTool, Select
from bokeh.plotting import gmap, figure
map_options = GMapOptions(lat=48.56, lng=2.25, map_type="roadmap", zoom=9)
plot = figure()
p = gmap("AIzaSyAjei7RGGqFEKP4Qu0ZUzi8vU-wkjeMz_k", map_options, title='Carte des accidents en Essonne')
# Create the sources for all types of accidents
mortel_source = ColumnDataSource(
data=dict(lat=lat_mortel,
lon=long_mortel,
hrmn=df.hrmn,
lum=df.lum,
atm=df.atm,
agglo=df.agglo)
)
grave_source = ColumnDataSource(
data=dict(lat=lat_grave,
lon=long_grave,
hrmn=df.hrmn,
lum=df.lum,
atm=df.atm,
agglo=df.agglo)
)
leger_source = ColumnDataSource(
data=dict(lat=lat_leger,
lon=long_leger,
hrmn=df.hrmn,
lum=df.lum,
atm=df.atm,
agglo=df.agglo)
)
p.circle(x="lon", y="lat", size=15, color='#00af46', alpha=0.7, source=leger_source, legend='Léger/Indemnes')
p.circle(x="lon", y="lat", size=15, color='#ffb200', alpha=0.7, source=grave_source, legend='Grave')
p.circle(x="lon", y="lat", size=15, color='#ea0012', alpha=0.7, source=mortel_source, legend='Mortel')
# Add the HoverTool object to figure p
hover = HoverTool(tooltips=[ ('hrmn', '@hrmn'), ('lum', '@lum'), ('atm', '@atm'), ('agglo', '@agglo')])
p.add_tools(hover)
show(p)
La N20 près d'Égly peut constituer une zone accidentogène dangereuse: 4 accidents mortels s'y sont produits, mais étrangement, on y constate peu d'accidents d'un autre type. Les accidents sans gravité sont majoritaires, et se produitsent sur les axes principaux.
Cette carte peut être améliorée, en affichant les catégories sous forme alphabétique plutôt que numérique, ou encore en ajoutant des menus déroulants pour permettre une exploration par l'utilisateur.
Modèle¶
Est-il possible de prédire la gravité d'un accident en fonction des données à notre disposition?
Nous allons proposer un modèle de prédiction sous le forme d'un arbre de décision. L'arbre de décision a déjà fait ses preuves dans la prédiction de la gravité d'un accident (surpassant apparemment les réseaux neuraux).
Un arbre de décision est une simple représentation graphique des solution possibles à une décision en fonction de certaines conditions. C'est un modèle extrêmement efficace pour identifier des clusters dans les données, et de ce fait, très utile pour la segmentation. Cet algorithme est utilisé en industrie pour construire des moteurs de recommandation, tels celui de Netflix proposant des films en fonction des goûts et de l'historique de l'utilisateur. Les arbres de décision ont une fâcheuse tendance à overfitter, c'est-à-dire à apprendre le dataset par coeur en pénalisant la capacité à prédire sur des données inconnues. Nous allons préférer une forêt d'arbres décisionnels, permettant de palier cette faiblesse. Il est nécessaire de préciser qu'une forêt est en général utilisée pour une classification complexe. Ici, nous avons 3 catégories (Léger/Indemnes, Grave, et Mortel), par conséquent une forêt peut s'avérer excessive par rapport au besoin présent.
Optimisation des données pour la prédiction¶
Le dataset tel quel nécessite d'être optimisé. Avant tout, nous allons retirer la colonne d'identifiants uniques, qui n'apporte rien à l'apprentissage:
dfp = df.drop(['Num_Acc'], axis=1)
Les données concernant la latitude et la longitude du lieu de l'accident sont certainement très utiles pour prédire la gravité d'un accident (une zone accidentogène provoque nécessairement un impact). Toutefois, l'utilisation de ces données requiert une transformation afin d'être comprises de l'algorithme. Malheureusement, nous n'avons aucune expérience sur ce sujet. Dans l'absolu, nous aimerions consacrer du temps à l'étude du traitement de ce type de données, ou bénéficier des conseils d'un collègue plus expérimenté.
Après recherches, il semble qu'une transformation trigonométrique soit possible, afin d'obtenir des coordonées en trois dimensions. Une autre solution consiste à définir des zones. L'agorithme AZP (Automated Zoning Procedures) pourrait alors être utilisé. Nous pourrions également essayer d'obtenir les codes postaux. Nous pourrions peut-être utiliser les adresses déjà fournies: explorons cette possibilité.
df_adr = pd.DataFrame(df.groupby('adr')['gravite'].count())
df_adr.sort_values(by=['gravite'], ascending=False).head()
Il semble que certaines routes soient plus accidentogènes que d'autres (notamment les autoroutes), nous allons donc, avec précaution, tenter de prendre en compte ces catégories dans notre algorithme de prédiction. Nous conserverons donc la colonne contenant les adresses, et supprimerons celle contenant les latitudes et longitudes.
dfp = dfp.drop(['lat', 'long'], axis=1)
De manière similaire, les données concernant les années, les mois, les jours et le moment de la journée constituent des données cycliques. Elles nécessitent également une transformation, afin que ce caractère cyclique soit pris en compte. Ici encore, une transformation trigonométrique semble être la solution appropriée. Là encore, nous aimerions étudier sérieusement la question de la transformation de ces données cycliques, ou bénéficier des conseils d'un collègue plus expérimenté.
Dans le cas présent, nous allons nous amputer, à regret, des données touchant au temps. Il va de soi que cela impactera considérablement notre modèle, comme il a été démontré que le mois avait un impact sur les accidents et leur gravité, de même que le moment de la journée.
dfp = dfp.drop(['an', 'mois', 'jour', 'hrmn'], axis=1)
Quelles colonnes sont donc encore présentes dans le dataset?
dfp.columns
Toutes ces colonnes sont catégoriques, et doivent donc être interprétées comme telles.
dfp.lum = dfp.lum.astype('category')
dfp.agglo = dfp.agglo.astype('category')
dfp.atm = dfp.atm.astype('category')
dfp.gravite = dfp.gravite.astype('category')
features_columns = (pd.get_dummies(data=dfp, columns=['adr'])).drop(['gravite'], axis=1)
features_columns.head()
target = pd.get_dummies(data=dfp.gravite, columns=['gravite'])
Forêt¶
Nous allons maintenant implémenter la forêt d'arbre décisionnels.
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split
X = features_columns
y = target
# Split data
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42, stratify=y)
## Instantiate the classifier
rf = RandomForestClassifier()
# Fit on the training data
rf.fit(X_train, y_train)
# Predict on the testing data
y_pred = rf.predict(X_test)
# Check the accuracy
rf.score(X_test, y_test)
from sklearn.model_selection import cross_val_score
# Instantiate the model
rf = RandomForestClassifier()
# Compute cross-validation results
cv_results = cross_val_score(rf, X, y, cv=10)
# Check the scores obtained
print(cv_results)
# Compute statistic of interest
print('Mean cross-validated score:', np.mean(cv_results))
Nous obtenons un score de 71.8%, ce qui signifie que notre algorithme a su prédire la gravité de l'accident dans environ 71 cas sur 100.
### ATTENTION — Cette optimisation du modèle est particulièrement longue
from sklearn.model_selection import GridSearchCV
# Set up the hyperparameter grid to test
param_grid = {'n_estimators': np.arange(1, 15),
'criterion': ['gini', 'entropy'],
'max_depth': np.arange(1, 5),
'min_samples_split': [2, 5, 10],
'min_samples_leaf': [1, 2, 4]}
# Instantiate the classifier
rf = RandomForestClassifier()
# Tune the hyperparameters
rf_cv = GridSearchCV(rf, param_grid, cv=10)
# Fit the GridSearchCV object to the data
rf_cv.fit(X, y)
# Get the best parameters
rf_cv.best_params_
# Get the best mean cross validation score over the fold that performs the best
rf_cv.best_score_
En optimisant les paramètres, nous obtenons un score de 75.8%, soit quatre points de mieux qu'auparavant.
Avec les différents axes d'amélioration énoncés précédemment, il devrait être possible d'obtenir un résultat encore plus fiable. Les données concernant le moment et le lieu précis de l'accident permettraient certainement de rendre ce modèle plus performant. D'autres données dont nous ne disposons pas encore, telles que la vitesse, l'âge ou l'expérience des conducteurs, le type de véhicules, la consommation ou non de substances altérant les réflexes, pourraient également se révéler utile: il peut s'avérer utile de le mentionner au client afin de voir si ces données pourraient être enregistrées à l'avenir.
Par ailleurs, une analyse géographique plus poussée pourrait permettre d'identifier des zones particulièrement accidentogène, et d'intervenir en conséquence (pose de radar, limitation de vitesse, rénovation ou modification de la chaussée...).