Liens de téléchargement des données nettoyées :
Données de mortalité due aux cancers
Séparation des bases entre base totale, base hommes et base femmes
Traitement des données manquantes
Statistiques descriptives sur les bases
Représentation des données et détection des outliers
Méthode des K-Means appliquée à l'ACP
Quelques représentations graphiques finales...
Séparation entre base d'apprentissage et base de test
Premier modèle : Régression multilinéaire
Régression multilinéaire sur la base à 13 variables
Deuxième modèle : Random Forest
C'est au moment où il fallait choisir un sujet pour notre projet de machine learning que nous avons reçu un email nous proposant de participer à un challenge open data sur le thème du cancer (le Challenge4Cancer). Notre engagement associatif (nous faisons tous les deux partie de Cheer Up!, association étudiante dont les membres visitent des jeunes atteints de cancer et les aident à réaliser des projets personnels) et notre conception de la data science comme une discipline qui peut avoir un impact positif sur la société nous a conduit à choisir ce challenge comme cadre pour notre projet.
C'est dans le cadre proposé par les organisateurs du challenge - la Paillasse et le laboratoire Roche - que nous avons étudiés les risques associés au cancer. Plus précisément, nous avons choisi d'axer nos recherches sur l'impact de la pollution et de certains facteurs comportementaux et sociaux sur la mortalité liée au cancer.
from IPython.display import Image
from IPython.core.display import HTML
Image(url= "http://wiki.epidemium.cc/images/thumb/2/20/EPDM_visuel_Partenaires.jpeg/800px-EPDM_visuel_Partenaires.jpeg")
Base totale : https://XChangeFile.ensae.fr/download.php?file=845916
Base réduite à 39 variables : https://XChangeFile.ensae.fr/download.php?file=703674
Base réduite à 13 variables : https://XChangeFile.ensae.fr/download.php?file=114045
Notre projet a la caractéristique d'avoir requis un temps de préparation considérable. Effectivement, nous avons dû tout d'abord effectuer un travail de recherche de bases de données pertinentes pour notre sujet (en tant que participants à un challenge Open Data, nous avions simplement décidé quel thème nous allions aborder avant d'avoir véritablement toutes les bases à notre disposition). Pour ce faire, nous avons utilisé des données récupérées sur différents sites Internet mais surtout celles proposées par le portail mis à notre disposition dans le cadre du challenge (http://data.epidemium.cc/dataset, plus de 20 000 data sets).
Plus précisément, voici les liens vers les datasets utilisés pour notre projet :
- Mortalité due au cancer en France : http://www.invs.sante.fr/applications/cancers/cancer_1985_2009/default.htm
- Données mensuelles météo scrapées au format pdf : https://donneespubliques.meteofrance.fr/?fond=produit&id_produit=129&id_rubrique=52
Abandonnées car trop peu pertinentes.
- Données météo saisonnières (abandonnées également) : http://data.epidemium.cc/organization/534fff8ba3a7292c64a77ed4
- Données sur le potentiel radon des communes françaises : http://data.epidemium.cc/dataset/53834c53a3a72906c7ec5c4c
- Données sur la pollution (EIDER): http://data.epidemium.cc/dataset/5369945ba3a729239d204413
- Données sur la pollution des eaux : http://data.epidemium.cc/dataset/53699e7ba3a729239d205eb0
- Données sur les OGM en 2007 : http://data.epidemium.cc/dataset/53699e7ba3a729239d205eb0
- Données sur les déchets dangereux : http://www.sinoe.org/contrib/module/catalogue/flex/exploitCatalogue.php?C_THEME_INDIC=1&C_SOUS_THEME_INDIC=2&ACCES=AT
- Données sur le chômage : http://data.epidemium.cc/dataset/5369a184a3a729239d206612
- Données démographiques générales : http://www.ecosante.fr/index2.php?base=DEPA&langs=FRA&langh=FRA&source=800341
- Données sur le vieillissement de la population : http://www.ecosante.fr/index2.php?base=DEPA&langs=FRA&langh=FRA&source=800190
- Données sur les médecins français : http://www.has-sante.fr/portail/jcms/c_1752855/fr/donnees-publiques-open-data
Nous avons des données sur la mortalité liée au cancer pour les deux sexes (plus précisément, le nombre de morts sur cinq ans dans chaque département français pour différents types de cancer, entre 1985 et 2009), c'est pourquoi nous avons créé trois bases : une pour les femmes, une pour les hommes et une pour les deux sexes (sachant que, bien entendu, chaque sexe est affecté par des cancers "généraux" et des cancers spécifiques, comme le cancer du sein). Ces trois bases et ces différents types de cancer pourront être utilisés par la suite pour affiner les modèles.
Du fait de la spécificité de notre sujet, nous avons formulé l'hypothèse que les facteurs de pollution affectaient la mortalité dûe au cancer mais sur le long terme. C'est pourquoi nous construisons un modèle de prédiction de la mortalité du cancer dans les départements français entre 2000 et 2009 à partir de données sur ces années ou sur les années précédentes.
Les données liées à la mortalité sont réparties dans des tableaux html sur plusieurs pages web du site de l'INVS. Les pages sont cependant indisponibles depuis quelques semaines. Voici malgré tout le programme qui a permis d'extraire les tableaux.
from bs4 import BeautifulSoup
import csv
from urllib.request import urlopen
basicurl = 'http://webcancer.invs.sante.fr/mortalite_8408/print.php?return_url=%2Fmortalite_8408%2Fdonnees_localisation%2Fevolution%2Flocal.php%3Flocalisation%3DL{location}&data=tab_CI_02_GR_L{location}_S{gender}'
locationlist = ["0" + str(i) for i in range(1,10)] + [str(i) for i in range(10,27)] + ["99"]
validchars = ' abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789'
def souptocsv(basicurl, validchars, location, gender):
try:
soup = BeautifulSoup(urlopen(basicurl.format(location=location,gender=gender)))
tables = soup.findAll('table') #cherche les tableaux dans la page
table1 = tables[0]
table2 = tables[1]
rows = []
for row in table1.find_all('tr'): #tr et td marquent les séparations
rows.append([val.text for val in row.find_all('td')])
title = gender + "".join([c for c in rows[1][0] if c in validchars])
headers = [header.text for header in table2.find_all('th')]
rows = []
for row in table2.find_all('tr'):
rows.append([val.text.encode('utf8') for val in row.find_all('td')])
folder = 'Mortalite/'
with open(folder + title + '.csv', 'w', newline='') as f:
writer = csv.writer(f)
writer.writerow(headers)
writer.writerows(row for row in rows if row)
except:
print(location + gender)
for i in locationlist:
souptocsv(basicurl, validchars, i, "1")
souptocsv(basicurl, validchars, i, "2")
Les données météorologiques de Météo France ne sont accessibles que dans les bulletins météorologiques publiés en pdf. Les données mensuels correspondent aux années 1999 à 2015.
#programme de téléchargement des données
def download_file(date):
import urllib
urllib.urlretrieve('https://donneespubliques.meteofrance.fr/donnees_libres/bulletins/BCM/{0}.pdf'.format(date), 'meteo/BCM{0}.pdf'.format(date))
print("Completed")
#Nous bouclons sur les dates de publication des bulletins météo
date = 199901
while date < 201511:
download_file(date)
if date%100==12:
date+=89
else:
date+=1
Les données se trouvent sur deux pages avec pour titre : "Résumé mensuel". Nous cherchons à savoir le numéro de ces pages pour chaque pdf.
def fnPDF_FindText(xFile, xString):
# xfile : chemin vers le pdf
# xString : string à chercher
import pyPdf, re
pagenumber = set()
pdfDoc = pyPdf.PdfFileReader(file(xFile, "rb"))
for i in range(2, pdfDoc.getNumPages()):
content = ""
content += pdfDoc.getPage(i).extractText() + "\n"
content1 = content.encode('ascii', 'ignore')
ResSearch = re.search(xString, content1)
if ResSearch is not None:
pagenumber.add(i)
return pagenumber
date = 199901
pagenumbers={}
while date < 201511:
if date%100 == 1:
print(str(date) + "**********************************************")
try:
pagenumbers[date]=fnPDF_FindText("meteo/BCM{0}.pdf".format(date), 'Rsum mensuel')
except ValueError: #Des erreurs de lecture sont possibles
print(date)
if date%100==12:
date+=89
else:
date+=1
Voici les deux principales fonctions qui convertissent une page pdf en html, puis parsent le html :
from pdfminer.pdfinterp import PDFResourceManager, PDFPageInterpreter
from pdfminer.converter import HTMLConverter
from pdfminer.converter import TextConverter
from pdfminer.layout import LAParams
from pdfminer.pdfpage import PDFPage
from cStringIO import StringIO
import re
import csv
import time
from BeautifulSoup import BeautifulSoup
def convert_pdf_to_html(path,pagenos):
rsrcmgr = PDFResourceManager()
retstr = StringIO()
codec = 'utf-8'
laparams = LAParams()
device = HTMLConverter(rsrcmgr, retstr, codec=codec, laparams=laparams)
fp = file(path, 'rb')
interpreter = PDFPageInterpreter(rsrcmgr, device)
password = ""
maxpages = 0
caching = True
for page in PDFPage.get_pages(fp, pagenos, maxpages=maxpages, password=password,caching=caching, check_extractable=True):
interpreter.process_page(page)
fp.close()
device.close()
string = retstr.getvalue()
retstr.close()
return string
def clean_html(alltext):
data = []
complementposition=1000
soup = BeautifulSoup(alltext)
divs = soup.findAll('div')
splitdiv=""
splitbr=""
for div in divs:
div_string = str(div)
splitdiv+=div_string
div_string = div_string.replace('<br />', '')
div_list = div_string.split('\n')
div_list = map(lambda x: x.strip('[').strip(']'),div_list)
if checkcomplement(div_list[0]):
complementposition = int(find_between(div_list[0],"top:", "px"))
record = []
for item in div_list:
record.append(item.split(':', 1))
data.append(record)
return(data, complementposition)
headers = "STATIONS,TN,TX,TNN,D,TXX,D,H.RR,RMAX,D,INST,FXI,D".split(",")
headersdico={"STATIONS" : str, "TN" : float, "TX" : float, "TNN" : float, "D" : int,
"TXX" : float, "H.RR" : float, "RMAX" : float, "INST" : int,
"FXI" : int}
Le code de la plupart des sous-fonctions ne sera pas montré. Il varie en fonction des pdfs (entre 1999 et 2015, le format des pdfs a beaucoup évolué et ainsi l'interprétation par le module pdfminer varie énormément). L'essentiel du travail repose pourtant sur les ajustements à faire.
- Dans beaucoup de pdfs, il y a des données manquantes, marquées comme des "*".
- Parfois, ces données manquantes sont complétées dans un pdf ultérieur grâce à un tableau nommé "compléments d'information". Il faut ainsi ignorer les données situées en dessous du titre "compléments d'information".
- Afin de reconnaître des données à récupérer dans le pdf, plusieurs stratégies ont été mises en place, selon les dates. On a établi une liste de caractères acceptées afin de former des floats, on a déterminé les positions verticales acceptées afin de délimiter le tableau. Enfin, la taille de la police servait également d'indice pour différencier les headers des données elles-mêmes.
- Pour certaines données, des solutions n'ont pas été trouvées. La colonne INST est irrécupérable, car les données manquantes sont trop nombreuses et marquées par des espaces vides. La détection des espaces vides est possible. Cependant, pdfminer ne permet pas de faire la différence entre 1 et 2 valeurs manquantes consécutives.
import csv
date = 199901
while date < 201511: #pour chaque pdf, nous travaillons sur la récupération des données
path ="meteo/BCM{0}.pdf".format(date)
outpath = "meteoout/BCM{0}.csv".format(date)
with open(outpath, "a") as f:
writer = csv.writer(f)
f.write(','.join(headers) + '\n')
for pagefrompdf in date_pages[date]: #pour chaque page retenue du pdf
donnees = []
alltext = convert_pdf_to_html(path, [pagefrompdf]) #on transforme la page en html
data, complementposition=clean_html(alltext) #on nettoie une partie du code html, en repérant
#la position verticale du "complément d'informations"
data, headersposition = wordsearch(data, '>STATIONS') #on cherche la position du header.
donneestation,fontsize=getstationsdata(data, headersposition, complementposition) # on récupère les
donnees.append(donneestation) # noms des stations
donnees+=trouveur(data, fontsize, complementposition) #récupération des données
try:
writer.writerows(clean_data(donnees, date))
except IndexError: #Les erreurs sont nombreuses. Connaître les dates des erreurs permet de faire
# des modifications ciblées sur des groupes de dates (avec pour hypothèse que
# pour des dates proches, le format des tables est semblable.)
print(date)
if date%100==12:
date+=89
else:
date+=1
Dans un second temps, nous avons dû nettoyer les bases extraites. Là, le site Dataiku et le Data Science Studio nous ont été d'une aide précieuse (nous y avons un accès gratuit le temps du challenge). Grâce au DSS, nous avons pu centraliser en un seul endroit tous nos datasets et les nettoyer sans avoir recours à chaque fois à du code Python. En outre, on avait la possibilité de visualiser en un instant les changements que nous avions effectués sur les données et leur mise en forme.
Néanmoins, ce fut également une étape très chronophage pour nous. Effectivement, le nombre de bases à traiter était élevé (plus d'une cinquantaine) et nous avons dû aussi faire face à des problèmes de mise en forme des données. Par exemple, certaines bases contenaient des données sur les communes et non pas sur les départements ; nous avons donc dû réfléchir à la meilleure manière d'agréger ces données et d'en tirer des informations représentatives au niveau départemental. En outre, nous avons dû sélectionner parmi la masse d'informations dont nous disposions les variables les plus pertinentes. A la fin, nous n'avons retenu que 1 500 variables.
En outre, le DSS a parfois eu des bugs qui nous ont fait perdre du temps : pour certaines bases, quand on voulait supprimer une colonne, seul le nom était supprimé et tout le reste des colonnes se décalait. Il a donc aussi fallu porter une attention particulière à ces bugs.
Enfin, une bonne partie de notre temps a aussi été consacrée au nettoyage général des données.
Voici une image montrant une partie des contributions de Joseph Lam :
from IPython.display import Image
PATH = 'C:\\Users\\lamjo_000\\Dropbox\\ENSAE C4C\\France\\Bases to merge\\Bases_recuperees\\'
Image(filename= PATH + 'JF_contributions.png')
Voici une image montrant une partie des contributions de Benoît Choffin :
Image(filename= PATH + 'BC_contributions.png')
Et enfin, voici une image qui montre notre travail sur Dataiku :
Image(filename= PATH + 'Dataiku_Flow.png')
Là encore, nous n'allons pas vous présenter tout le code qui a servi à créer ces trois bases, mais nous allons vous en montrer un aperçu :
dico1={} #dictionnaire qui contiendra les labels à garder dans la base de départ
for label in df_femmes.columns.values:
date=label[-4:]
labelsplit = label.split('_')
if 'HOMME' in label:
df_femmes.drop(label, axis=1,inplace=True)
elif 'FEMME' in label:
dico1[label]=labelsplit[0] + "_" + labelsplit[2]
elif ('FEMME' not in label) & date.isdigit():
df_femmes[label] = df_femmes[label] * df_femmes_pop[date] #ligne servant à mutliplier les variables explicatives
#par la population féminine dans tel département cette année
df_femmes.rename(columns=dico1, inplace=True)
df_femmes.to_csv(join(directory,"Souspop\\df_femmes.csv"))
#on fait un processus symétrique pour le cas des hommes
df_femmes_total=pd.read_csv(join(directory,"Souspop\\df_femmes.csv"),index_col=0)
df_hommes_total=pd.read_csv(join(directory,"Souspop\\df_hommes.csv"),index_col=0)
for label in df_femmes_total.columns.values:
labelsplit = label.split('_')
if labelsplit[0] in ['Sein', 'Ovaires']: #pour le total, on ne garde que les cancers en commun aux deux sexes
df_femmes_total.drop(label, axis=1,inplace=True)
for label in df_hommes_total.columns.values:
labelsplit = label.split('_')
if labelsplit[0] in ['Larynx', 'Prostate']:
df_hommes_total.drop(label, axis=1,inplace=True)
df_total = df_hommes_total + df_femmes_total
df_total.to_csv(join(directory,"Souspop\\df_total.csv"))
Nous devons traiter les données manquantes. Pour cela, nous avons considéré plusieurs stratégies. Dans tous les cas, lorsqu'une colonne comportait uniquement des valeurs manquantes, elle était supprimée.
from os.path import join
directory = "C:\\Users\\lamjo_000\\Dropbox\\ENSAE C4C\\France\\Bases to merge\\Total"
base_hommes = pd.read_csv(join(directory,"Souspop\\df_hommes.csv"),index_col=0)
base_femmes = pd.read_csv(join(directory,"Souspop\\df_femmes.csv"),index_col=0)
base_tot = pd.read_csv(join(directory,"Souspop\\df_total.csv"),index_col=0)
base_H_NaN = base_hommes.dropna(axis=1,how='all') #nous retirons les colonnes vides
base_F_NaN = base_femmes.dropna(axis=1,how='all')
base_T_NaN = base_tot.dropna(axis=1,how='all')
Une première méthode naïve est de supprimer toutes les colonnes où se trouvent au moins une valeur manquante. Elle se justifie par le nombre très important de colonnes que nous avons.
base_H_NaN = base_hommes.dropna(axis=1,how='any') #nous retirons les colonnes
base_F_NaN = base_femmes.dropna(axis=1,how='any')#avec une valeur manquante
base_T_NaN = base_tot.dropna(axis=1,how='any')
Une seconde méthode non moins naïve à laquelle nous avons pensé serait de prendre la moyenne de la colonne et de remplacer les valeurs manquantes par celle-ci. Scikitlearn possède une feature pour cela, il s'agit d'Imputer. Le problème de cette méthode est que cela introduit du biais dans nos modèles.
from sklearn.preprocessing import Imputer
import numpy as np, pandas as pd
imp1 = Imputer(missing_values='NaN', strategy='mean', axis=0)
imp2 = Imputer(missing_values='NaN', strategy='mean', axis=0)
imp3 = Imputer(missing_values='NaN', strategy='mean', axis=0)
imp1.fit(base_hommes)
imp2.fit(base_femmes)
imp3.fit(base_tot)
base_hommes_v2=imp1.transform(base_hommes)
base_femmes_v2=imp2.transform(base_femmes)
base_tot_v2=imp3.transform(base_tot)
df_hommes_v2 = pd.DataFrame(base_hommes_v2, columns=base_H_NaN.columns.values, index = base_H_NaN.index)
df_femmes_v2 = pd.DataFrame(base_femmes_v2, columns=base_F_NaN.columns.values, index = base_F_NaN.index)
df_tot_v2 = pd.DataFrame(base_tot_v2, columns=base_T_NaN.columns.values, index = base_T_NaN.index)
Enfin, la méthode qui nous a paru la plus puissante est celle de la Multivariate Imputation by Chained Equations. Nous l'avons réalisée à l'aide du package R : MICE.
# Code R
library(mice)
data = read.csv("C:/Users/lamjo_000/Dropbox/ENSAE C4C/France/Bases to merge/Total/Souspop/base_reduitenopop.csv", header = TRUE)
Il ne suffit cependant pas de simplement lancer la fonction. En effet, cela nécessite une préparation auparavant des données. Ainsi, lorsqu'il y a trop de valeurs manquantes pour une colonne, le remplissage semble périlleux.
# Code R
new_DF <- data[ lapply( data, function(x) sum(is.na(x)) / length(x) ) < 0.05 ] #On supprime les colonnes avec plus de 5% de NA
new_DF$Departement <- NULL #Comme ce sont des strings, cela peut causer des problèmes.
drops <- c("TOUS.CANCERS_1985.89","TOUS.CANCERS_1990.94", "TOUS.CANCERS_1995.99")
new_DF <- new_DF[,!(names(new_DF) %in% drops)]
Puis, lors des régressions, si nous avons des régresseurs linéairement dépendants, cela peut poser problème. C'est pourquoi nous cherchons les couples de variables avec un coefficient de corrélation trop élevé et nous ne gardons que l'une des deux.
# Code R
tmp <- cor(new_DF) #On travaille sur la matrice de corrélation, sans répétition (on considère qu'un triangle)
tmp[upper.tri(tmp)] <- 0
diag(tmp) <- 0 #On ne souhaite pas supprimer toutes les variables puisqu'elles ont une corrélation = 1 avec elles-mêmes
new_DF90<- apply(tmp,2,function(x) any(x > 0.90))
new_DF90[is.na(new_DF90)] <- FALSE #Les valeurs manquantes seront gardées
new_DF <- new_DF[,!new_DF90]
new_DF[,"TOUS.CANCERS_2005.09"] = data[,"TOUS.CANCERS_2005.09"]#On remet notre variable à modéliser
Enfin, nous continuons à retirer des variables afin d'obtenir une matrice de rang le nombre de variables. Pour cela, nous procédons à une stepwise selection. Cela permettra également de ne garder que les variables essentielles pour des régressions efficaces.
# Code R
library(MASS)
fit <- lm(TOUS.CANCERS_2005.09 ~ .,data=new_DF)
step <- stepAIC(fit, direction="backward")
step$anova # display results
#Nous suivons les recommandations de la Stepwise Regression afin de garder seulement les variables nécessaires.
DF_stepped <- new_DF[, (names(new_DF) %in% c("TOUS.CANCERS_2005.09","Arsenic...Nombre.d.etablissements.concernes_2001","Arsenic...Quantites.annuelles.emises_2001","Arsenic...Quantites.annuelles.emises_2007","Arsenic...Quantites.annuelles.emises_2008","Arsenic...Quantites.annuelles.emises_2009","Arsenic..et.ses.composes....Nombre.d.etablissements.concernes_2002","Arsenic..et.ses.composes....Quantites.annuelles.emises_2005","Arsenic..et.ses.composes....Quantites.annuelles.emises_2006","Cadmium...Quantites.annuelles.emises_2002","Cadmium...Quantites.annuelles.emises_2006","Cadmium...Nombre.d.etablissements.concernes_2008","Vente.de.cigarettes.en.grammes.par.habitant"))]
Nous appliquons enfin la Multivariate Imputation by Chained Equations afin de compléter les valeurs manquantes grâce à celles présentes.
# Code R
tempData <- mice(DF_stepped,m=5,maxit=10,meth='pmm',seed=400)
summary(tempData)
completedData <- complete(tempData,1)
completedData[,"Departement"] <- data[,"Departement"]
Par application des méthodes précédentes, nous obtenons une base avec 13 régresseurs et une base avec 39 régresseurs. Elles sont au nombre de deux car pour la base avec 39 variables explicatives, nous sommes partis sans a priori avec la base que nous avions nettoyée, tandis qu'avec la base à 13 variables explicatives, nous avons effectué une présélection des variables en fonction des facteurs cancérigènes connus (cf. http://www.cancer.gov/about-cancer/causes-prevention/risk).
Ce sont ces bases que nous utiliserons par la suite dans nos analyses. Nous n'ajouterons pas la population à ces bases, car nous avons remarqué que cela réduisait le rang de notre matrice lors de la Multivariate Imputation by Chained Equations.
from os.path import join
import numpy as np, pandas as pd
import matplotlib.pyplot as plt
plt.style.use('ggplot')
import statsmodels.api as sm
%pylab inline
directory = "C:\\Users\\lamjo_000\\Dropbox\\ENSAE C4C\\France\\Bases to merge\\Total"
df_total = pd.read_csv(join(directory,"Base_totale.csv"),index_col=0) # Nous définissons les bases
df_39 = pd.read_csv(join(directory,"base_reduite_39.csv"),index_col=39) # accessibles en téléchargement
df_13 = pd.read_csv(join(directory,"base_reduite_13.csv"),index_col=13) # plus haut
Voici un aperçu des tableaux sur lesquels nous avons travaillé :
df_total.head()
Il serait inutile de faire un descriptif détaillé de chacune de nos variables, car elles sont très nombreuses. Néanmoins, nous allons vous en faire une description générale :
Toutes nos variables sont continues. Nos variables explicatives correspondent le plus souvent à des relevés/mesures dans les différents départements français et sur plusieurs années (pour chaque variable initiale, nous avons créé autant de variables qu'il y avait d'années disponibles pour cette variable), de produits polluants comme l'arsenic. Mais nous avons également des données comme la vente de cigarettes en grammes par habitant et le taux de chômage, qui ne sont pas reliées directement à la pollution.
Nos variables expliquées correspondent quant à elles au nombre de morts pour différents types de cancer (pancréas, estomac,...) pour tous les départements entre 1985 et 2009. Ce sont des données quinquennales. Pour les bases réduites, nous n'avons retenu que la variable "Tous Cancers 2005-09" comme variable expliquée.
Pour finir, nous avons supprimé toutes les variables explicatives postérieures à 2009 (nous n'avons pas de données sur le cancer à partir de ce moment-là) et toutes les variables expliquées antérieures à 2000 (trop peu de données explicatives pour ces années).
colon_sorted=df_total[['Colon et rectum_2005-09']].sort_values('Colon et rectum_2005-09',ascending=False)
hist1 = colon_sorted.plot(kind = "bar", figsize=(18,6))
hist1.set_xlabel("Departements", fontsize=16)
hist1.set_title("Nombre de morts du cancer du colon/rectum entre 2005 et 2009", fontsize=16)
hist1.legend().set_visible(False)
estomac_sorted=df_total[['Estomac_2005-09']].sort_values('Estomac_2005-09',ascending=False)
hist2 = estomac_sorted.plot(kind = "bar", figsize=(18,6))
hist2.set_xlabel("Departements", fontsize=16)
hist2.set_title("Nombre de morts du cancer de l'estomac entre 2005 et 2009", fontsize=16)
hist2.legend().set_visible(False)
Ces deux histogrammes font apparaître plusieurs choses.
Tout d'abord, il y a une forte disparité entre les départements en fonction du nombre de morts pour un type de cancer. Par exemple, le département 70 a eu un peu moins de 20 morts à cause du cancer de l'estomac entre 2005 et 2009, contre plus de 160 pour le département 59.
De même, il y a des disparités entre les types de cancer eux-mêmes : le classement des départements n'est pas identique si l'on change de type de cancer et certains cancers causent plus de morts que d'autres. Ainsi, le département 76 est en septième position pour la mortalité liée au cancer du colon alors qu'il est en cinquième position pour le cancer de l'estomac.
La disparité entre les départements peut s'expliquer en partie par la taille du département. En effet, on peut légitimement penser que plus un département a une population nombreuse, plus il a un taux d'individus atteints par le cancer élevé et plus la mortalité liée au cancer de ce département est élevée. En outre, le vieillissement d'une population est aussi un facteur influant sur la mortalité liée au cancer.
Nous utiliserons la variable de la population dans la suite de ce mémoire pour tenter d'affiner nos résultats et prendre en compte l'effet-taille qui joue dans notre étude.
from sklearn.decomposition import PCA
from sklearn.preprocessing import scale
from pandas.tools.plotting import scatter_matrix
dim_reduction = PCA() #instanciation de la PCA
Xc = dim_reduction.fit_transform(scale(df_13)) #application de la PCA à la base réduite à 13 variables
#Graphique représentant les parts de variance expliquée par les axes principaux
plt.figure(figsize=(18,6))
plt.bar(np.arange(len(dim_reduction.explained_variance_ratio_))+0.5, dim_reduction.explained_variance_ratio_)
plt.title("Variance expliquee par les axes principaux")
plt.show()
print(dim_reduction.explained_variance_ratio_)
df_pca = pd.DataFrame(Xc, columns=['comp_'+str(j+1) for j in range(13)])
#Représentation des départements sur le premier et le dernier plan factoriel pour déterminer les outliers
plt.figure(figsize=(18,6))
plt.scatter(df_pca['comp_1'], df_pca['comp_2'])
plt.title("Projection des departements sur le premier plan factoriel")
for label, x, y in zip(df_13.index, df_pca['comp_1'], df_pca['comp_2']):
plt.annotate(
label,
xy = (x, y), xytext = (-10, 10),
textcoords = 'offset points', ha = 'right', va = 'bottom',
bbox = dict(boxstyle = 'round,pad=0.5', fc = 'yellow', alpha = 0.5),
arrowprops = dict(arrowstyle = '->', connectionstyle = 'arc3,rad=0'))
plt.figure(figsize=(18,6))
plt.scatter(df_pca['comp_12'], df_pca['comp_13'])
plt.title("Projection des departements sur le dernier plan factoriel")
for label, x, y in zip(df_13.index, df_pca['comp_12'], df_pca['comp_13']):
plt.annotate(
label,
xy = (x, y), xytext = (-10, 10),
textcoords = 'offset points', ha = 'right', va = 'bottom',
bbox = dict(boxstyle = 'round,pad=0.5', fc = 'yellow', alpha = 0.5),
arrowprops = dict(arrowstyle = '->', connectionstyle = 'arc3,rad=0'))
Par le critère du coude, nous retenons les deux premiers axes.
Cette ACP nous permet non seulement de représenter les départements, mais aussi de détecter les outliers, les départements qui diffèrent de la norme.
Ces outliers forment des groupes et l'on peut d'ores et déjà en regrouper certains :
- les départements 62 (Pas-de-Calais) et 59 (Nord)
- les départements 67 (Bas-Rhin) et 68 (Haut-Rhin)
- les départements 55 (Meuse), 88 (Vosges), 54 (Meurthe-et-Moselle) et 57 (Moselle)
- les départements 52 (Haute-Marne), 10 (Aube), 51 (Marne) et 08 (Ardennes)
Le dernier plan factoriel nous permet de discerner certains outliers difficiles à détecter : ici, on peut noter les départements 27 (Eure) et 76 (Seine-Maritime).
df_13_array=df_13.values
Y_tot = df_13.ix[:,0].values
X_tot = df_13.ix[:,1:].values
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_samples, silhouette_score
import matplotlib.pyplot as plt
import matplotlib.cm as cm
range_n_clusters = [5,6,7,8,9] #La liste des nombres de clusters testés
for n_clusters in range_n_clusters:
#Nous allons créer 2 subplots côte à côte
fig, (ax1, ax2) = plt.subplots(1, 2)
fig.set_size_inches(18, 7)
# Le premier est le tracé de la silhouette
ax1.set_xlim([-0.1, 1])
ax1.set_ylim([0, df_13_array.shape[0] + (n_clusters + 1) * 10])
# On initialise l'algorithme sur les composantes principales de l'ACP précédente
clusterer = KMeans(init='k-means++', n_clusters=n_clusters, n_init=10)
cluster_labels = clusterer.fit_predict(df_pca)
# On affiche le silhouette_score.
# Cela donne une idée du niveau de séparation des clusters formés.
silhouette_avg = silhouette_score(df_pca, cluster_labels)
print("For n_clusters =", n_clusters,
"The average silhouette_score is :", silhouette_avg)
sample_silhouette_values = silhouette_samples(df_pca, cluster_labels)
y_lower = 10
for i in range(n_clusters):
# On rassemble les silhouette_scores pour le ith cluster, puis on trie.
ith_cluster_silhouette_values = sample_silhouette_values[cluster_labels == i]
ith_cluster_silhouette_values.sort()
size_cluster_i = ith_cluster_silhouette_values.shape[0]
y_upper = y_lower + size_cluster_i
color = cm.spectral(float(i) / n_clusters)
ax1.fill_betweenx(np.arange(y_lower, y_upper),
0, ith_cluster_silhouette_values,
facecolor=color, edgecolor=color, alpha=0.7)
# On indique le numéro du cluster pour chaque silhouette
ax1.text(-0.05, y_lower + 0.5 * size_cluster_i, str(i))
# On déplace le curseur vers le haut
y_lower = y_upper + 10 # 10 for the 0 samples
ax1.set_title("The silhouette plot for the various clusters.")
ax1.set_xlabel("The silhouette coefficient values")
ax1.set_ylabel("Cluster label")
# On trace le silhouette_score moyen.
ax1.axvline(x=silhouette_avg, color="red", linestyle="--")
#ax1.set_yticks([])
ax1.set_xticks([-0.1, 0, 0.2, 0.4, 0.6, 0.8, 1])
# Le 2e plot donne la projection des clusters
colors = cm.spectral(cluster_labels.astype(float) / n_clusters)
ax2.scatter(df_13_array[:, 0], df_13_array[:, 1], marker='.', s=30, lw=0, alpha=0.7,
c=colors)
centers = clusterer.cluster_centers_
# On marque le centre du cluster
ax2.scatter(centers[:, 0], centers[:, 1],
marker='o', c="white", alpha=1, s=200)
for i, c in enumerate(centers):
ax2.scatter(c[0], c[1], marker='$%d$' % i, alpha=1, s=50)
ax2.set_title("The visualization of the clustered data.")
ax2.set_xlabel("Feature space for the 1st feature")
ax2.set_ylabel("Feature space for the 2nd feature")
plt.suptitle(("Silhouette analysis for KMeans clustering on sample data "
"with n_clusters = %d" % n_clusters),
fontsize=14, fontweight='bold')
plt.show()
Les graphes ci-dessus présentent les silhouettes pour différents choix de nombres de clusters. En parallèle, on peut voir les départements projetés sur les deux premiers axes (les deux premières variables) et les clusters formés.
La silhouette est une mesure de la distance qui sépare chaque point d'un cluster des points des clusters environnants. Plus la silhouette est proche de 1 et plus les clusters sont bien délimités entre eux.
En appliquant le K-Means sur l'ACP, nous appliquons les différents critères se basant sur la silhouette_score, l'équilibrage des classes et leur bonne séparation. Nous retenons ainsi un découpage en 7 classes. La répartition des départements dans ces classes est explicitée dans le tableau suivant.
clusterer = KMeans(init='k-means++', n_clusters=7, n_init=10)
clusterer.fit(df_pca)
pd.DataFrame(data=clusterer.labels_, index=df_13.index.values, columns=["Classe"])
Nous examinons les coordonnées décrivant chaque classe.
clusterer.cluster_centers_
Ainsi, pour les outliers par exemple, nous pouvons observer comment ils se situent par rapport aux classes précédemment décrites et par rapport à la moyenne de la base en général.
outliers_base=df_13.loc[['62','59','67','68','55','88','54','57','52','10','51','08','27','76']]
study_mean = df_13.mean(axis=0)
outliers_deviation=outliers_base.copy()
for dept in range(outliers_base.shape[0]):
outliers_deviation.ix[dept,:] = study_mean - outliers_base.ix[dept,:]
outliers_deviation.ix["std",:] = df_13.std(axis=0) #std of whole study base
outliers_deviation
Ce tableau ci-dessus nous permet, pour les départements repérés grâce à l'ACP réalisée plus haut, de savoir pourquoi ce sont des outliers et plus précisément, ce qui fait leur spécificité au sein des autres départements.
Nous n'allons pas analyser chacun des outliers dans le détail, mais nous donnons ci-après les spécificités de certains d'entre eux.
On peut par exemple étudier les deux départements du Nord-Pas-de-Calais, le Nord (59) et le Pas-de-Calais (62). Ces deux départements se distinguent par une mortalité liée au cancer et des émissions de produits polluants (arsenic, cadmium,...) bien plus importantes que la moyenne (un signe négatif signifie que la valeur prise par le département pour cette variable est supérieure à la moyenne). Le Nord-Pas-de-Calais est en effet une ancienne région industrielle, aujourd'hui sinistrée par la désindustrialisation. On peut supposer que dans cette région, de nombreux produits polluants ont été par le passé dégagés ; nous avons de plus fait l'hypothèse que la survenue du cancer dépendait de l'exposition prolongée à des facteurs de risque. On peut donc proposer comme interprétation que dans cette région la pollution joue un rôle plus important que dans d'autres (du fait de la présence accrue de la pollution) dans la mortalité liée au cancer.
import seaborn as sns
sns.set(style="white")
corr1=df_13.corr() #création de la matrice de corrélation sur notre base réduite à 13 variables
f, ax = plt.subplots(figsize=(11, 9))
cmap = sns.diverging_palette(220, 10, as_cmap=True)
sns.heatmap(corr1, vmax=.8, square=True)
Voici de plus la matrice de corrélation de la base réduite à 13 variables :
corr1
sns.set(style="white")
corr2=df_39.corr() #création de la matrice de corrélation sur notre base réduite à 13 variables
f, ax = plt.subplots(figsize=(11, 9))
cmap = sns.diverging_palette(220, 10, as_cmap=True)
sns.heatmap(corr2, vmax=.8, square=True)
Et voici la matrice de corrélation pour notre base réduite à 39 variables :
corr2.head()
Nous avons ôté précédemment toutes les variables explicatives qui semblaient trop corrélées entre elles (par exemple une même variable sur plusieurs années, potentiellement sujette à de l'autocorrélation) car cela pouvait fausser nos régressions. Néanmoins, il faut tout de même que les variables explicatives soient corrélées avec la variable expliquée afin que nos modèles soient pertinents.
C'est le cas ici pour la plupart de nos variables : par exemple, entre le nombre total de morts du cancer entre 2005 et 2009 et le tonnage de collectes de déchets dangereux en 2009, il y a un coefficient de corrélation de 0,48. Cela peut sembler peu, mais il ne faut pas oublier que nos variables sont toutes agrégées et que la pollution n'est pas a priori un des facteurs principaux de la mortalité dûe au cancer.
fig, axs = plt.subplots(1, 3, sharey=True)
df_39.plot(kind='scatter', x='Vente.de.cigarettes.en.grammes.par.habitant', y='TOUS.CANCERS_2005.09', ax=axs[0], figsize=(16, 8))
df_39.plot(kind='scatter', x='Consommation.de.combustibles.fossiles.dans.les.reseaux.de.chaleur_2007', y='TOUS.CANCERS_2005.09', ax=axs[1])
df_39.plot(kind='scatter', x='Mercure..et.ses.composants....Quantites.annuelles.emises_2000', y='TOUS.CANCERS_2005.09', ax=axs[2])
Sur le premier graphique, qui donne la mortalité dûe au cancer entre 2005 et 2009 dans les départements français en fonction de la vente de cigarettes en grammes par habitant, il semble y avoir une relation linéaire entre les deux variables. Dans le deuxième graphique (Mortalité en fonction de la consommation de combustibles fossiles dans les réseaux de chaleur en 2007), la relation est moins nette, mais il est possible qu'elle existe. En revanche, dans le troisième graphique (Quantités annuelles émises de mercure et de ses composants en 2000), la relation semble totalement inexistante.
Notre approche a bien entendu des limites, parmi lesquelles figurent l'agrégation forte de nos données et notre choix de se restreindre quasi-exclusivement à des variables décrivant la pollution. De même, la présence de valeurs manquantes n'améliore pas nos prédictions.
Pour autant, nous allons par la suite tirer le maximum d'informations possibles de nos données et essayer de prédire au mieux la mortalité liée au cancer à partir des variables en notre possession.
from sklearn import cross_validation
Ici, nous séparons notre base entre base d'apprentissage et base de test. La première méthode utilisée est la plus simple, elle consiste à séparer la base de test de la base d'apprentissage de façon aléatoire. On choisit ici une taille de l'échantillon de test égale à 40% de la base totale.
#on utilise une première méthode simple de splitting
X_train1, X_test1, y_train1, y_test1 = cross_validation.train_test_split(df_13.ix[:,1:], df_13.ix[:,0], test_size=0.4, random_state=0)
test_index1 = X_test1.index.values.tolist()
train_index1 = X_train1.index.values.tolist()
La deuxième méthode est celle du KFold. En effet, la méthode du KFold consiste à séparer la base en k échantillons de même taille et effectuer k-1 fois le processus d'apprentissage et de test, mais à chaque fois en changeant l'échantillon de test. Cette technique doit permettre d'améliorer la robustesse de nos prédictions.
#et une deuxième, celle du KFold
kf = cross_validation.KFold(df_13.shape[0],n_folds=3,shuffle=True,random_state=True)
Comme nous étudions uniquement des variables continues, nous nous tournons tout naturellement vers la méthode supervisée pour les variables continues : la régression. Notre premier modèle sera donc la régression multilinéaire. Nous souhaitons en effet tester des modèles à la complexité croissante pour voir si les prédictions s'améliorent à raison de la complexité du modèle.
#Nous appliquons notre modèle à la base réduite à 13 variables dans un premier temps
from sklearn import linear_model
regr = linear_model.LinearRegression() #création de l'objet régression linéaire
regr.fit(X_train1, y_train1) #application à notre base splittée de la manière la plus simple
print('Coefficients: %s' %regr.coef_)
print("Residual sum of squares: %.2f"
% np.mean((regr.predict(X_test1) - y_test1) ** 2))
print('Variance score: %.2f' % regr.score(X_test1, y_test1))
from matplotlib.legend_handler import HandlerLine2D
y_pred1 = regr.predict(X_test1)
xplot = [i for i in range(y_pred1.shape[0])]
plt.figure(figsize=(18,6))
plt.title('Comparaison entre la valeur attendue et la valeur predite avec le modele de regression multilineaire (13 variables)')
plt.xticks(xplot, df_13.index[test_index1].values.tolist())
plt.scatter(xplot, y_test1.tolist(), color='red', label='Valeur attendue')
plt.plot(xplot, y_pred1.tolist(), color='blue', label='Valeur predite')
plt.legend()
plt.show()
Le score de la variance (ou coefficient de détermination) est égal à -0.27, ce qui est très mauvais. Peut-être avons nous oublié des régresseurs importants : nous allons essayer la régression multilinéaire sur la base réduite à 39 variables.
X_train2, X_test2, y_train2, y_test2 = cross_validation.train_test_split(df_39.ix[:,1:], df_39.ix[:,0], test_size=0.4, random_state=0)
test_index2 = X_test2.index.values.tolist()
train_index2 = X_train2.index.values.tolist()
#séparation entre base d'apprentissage et base de test pour la base réduite à 39 variables
regr1 = linear_model.LinearRegression() #création de l'objet régression linéaire
regr1.fit(X_train2, y_train2) #application à notre base splittée de la manière la plus simple
print('Coefficients: %s' %str(regr1.coef_))
print("Residual sum of squares: %.2f"
% np.mean((regr1.predict(X_test2) - y_test2) ** 2))
print('Variance score: %.2f' % regr1.score(X_test2, y_test2))
y_pred2 = regr1.predict(X_test2)
xplot = [i for i in range(y_pred2.shape[0])]
plt.figure(figsize=(18,6))
plt.title('Comparaison entre la valeur attendue et la valeur predite avec le modele de regression multilineaire (39 variables)')
plt.xticks(xplot, df_39.index[test_index2].values.tolist())
plt.scatter(xplot, y_test2.tolist(), color='red', label='Valeur attendue')
plt.plot(xplot, y_pred2.tolist(), color='blue', label='Valeur predite')
plt.legend()
plt.show()
La situation a l'air de s'améliorer. Le coefficient de détermination augmente (mécaniquement) du fait de l'ajout de régresseurs. Notre somme des carrés résiduels, elle aussi, diminue à environ 740 000. Enfin, la mortalité de certains départements est mieux prédite avec les 39 variables qu'avec seulement 13 régresseurs. Par exemple, les départements 05, 84 et 55, dont les prédictions étaient au départ très éloignées de leurs vraies valeurs, sont maintenant bien prédits. D'autres départements (comme le 78), moins nombreux, sont moins bien prédits.
Mais ce qui est ici remaraquable, c'est qu'en voyant nos deux graphes, on se rend compte que notre modèle prédit bien en général, mais qu'il a beaucoup de mal à prédire les outliers. En effet, les départements qui lui posent le plus de problèmes (55, 68, 62, 76) sont en règle générale les outliers que nous avons précédemment déterminés. Cela justifie donc l'utilisation de modèles de machine learning pour les départements "dans la norme" mais une étude spécifique pour chacun des outliers.
Essayons maintenant d'utiliser la méthode du K-Fold pour améliorer notre score de prédiction et la robustesse de nos résultats.
kf2 = cross_validation.KFold(df_39.shape[0],n_folds=3,shuffle=True,random_state=True)
df_39_array=df_39.values
Y_tot1 = df_39.ix[:,0].values
X_tot1 = df_39.ix[:,1:].values
regr3 = linear_model.LinearRegression()
listy_pred1=pd.DataFrame(data = np.zeros((df_39.shape[0],1)), index = df_13.index)
for train_index, test_index in kf2:
X_train3, X_test3 = X_tot1[train_index], X_tot1[test_index]
y_train3, y_test3 = Y_tot1[train_index], Y_tot1[test_index]
regr3.fit(X_train3,y_train3)
temp_df = pd.DataFrame(data = np.zeros((df_39.shape[0],1)), index = df_39.index)
temp_df.ix[test_index] = regr3.predict(X_test3).reshape(32,1)
listy_pred1.ix[test_index] += temp_df.ix[test_index]
y_pred1 = listy_pred1.values
xplot = [i for i in range(listy_pred1.shape[0])]
plt.figure(figsize=(18,6))
plt.xticks(xplot, df_39.index.values.tolist())
plt.scatter(xplot, Y_tot1.tolist(), color='red')
plt.plot(xplot, y_pred1.tolist(), color='blue')
plt.show()
print('Coefficients: %s' %regr.coef_)
print("Residual sum of squares: %.2f"
% np.mean((y_pred1 - Y_tot1) ** 2))
print(metrics.mean_squared_error(df_13['TOUS.CANCERS_2005.09'].tolist(), np.array(y_pred1)))
Par le K-Fold, nous obtenons bien des courbes plus proches de la base initiale.
L'amélioration des résultats passe par l'ajustement des hyperparamètres. Leur choix se fait par intuition, et surtout au travers de différents essais comparés. Par exemple, maximiser le nombre d'arbres n_estimators est positif pour la régression avec un random forest. Cependant, nous souhaitons optimiser avec un arbitrage en temps notamment. Le nombre d'arbres devrait ainsi se situer entre 64 et 128 d'après l'article suivant : https://www.researchgate.net/publication/230766603_How_Many_Trees_in_a_Random_Forest. Des recherches nous poussent à prendre : max_features inclus [30,50%] du nombre de variables, max_depth aux alentours de 5-10 avant d'essayer plus élevé, et min_samples_leaf > 1.
En ce qui concerne les essais, il est possible de lancer une comparaison exhaustive entre des choix de paramètres définis à l'avance grâce à GridSearch.
from sklearn import datasets
from sklearn.grid_search import GridSearchCV
from sklearn.ensemble import RandomForestRegressor
tuned_parameters = [{'max_depth': [5,6,7,8,9,10], 'max_features':[0.3,0.4,0.5]}]
for train_index, test_index in kf2:
clf = GridSearchCV(RandomForestRegressor(n_estimators = 128), tuned_parameters, cv=5,
scoring='mean_squared_error')
X_train, X_test = X_tot1[train_index], X_tot1[test_index]
y_train, y_test = Y_tot1[train_index], Y_tot1[test_index]
clf.fit(X_train, y_train)
print("Best parameters set found on development set:")
print()
print(clf.best_params_)
print()
print("Grid scores on development set:")
print()
for params, mean_score, scores in clf.grid_scores_:
print("%0.3f (+/-%0.03f) for %r"
% (mean_score, scores.std() * 2, params))
print()
print("Detailed classification report:")
print()
print("The model is trained on the full development set.")
print("The scores are computed on the full evaluation set.")
print()
y_true, y_pred = y_test, clf.predict(X_test)
Nous prenons : {'max_features': 0.4, 'max_depth': 7}.
from sklearn import metrics
from sklearn.ensemble import RandomForestRegressor
from sklearn.grid_search import GridSearchCV
clf3 = RandomForestRegressor(n_estimators=128,max_features=0.4,max_depth=7)
listy_pred2=pd.DataFrame(data = numpy.zeros((df_39.shape[0],1)), index = df_39.index)
for train_index, test_index in kf2:
X_train, X_test = X_tot1[train_index], X_tot1[test_index]
y_train, y_test = Y_tot1[train_index], Y_tot1[test_index]
clf3.fit(X_train,y_train)
temp_df = pd.DataFrame(data = numpy.zeros((df_39.shape[0],1)), index = df_39.index)
temp_df.ix[test_index] = clf3.predict(X_test).reshape(32,1)
listy_pred2.ix[test_index] += temp_df.ix[test_index]
y_pred2 = listy_pred2.values
xplot = [i for i in range(y_pred2.shape[0])]
plt.figure(figsize=(18,6))
plt.title('Comparaison entre la valeur attendue et la valeur predite avec le modele du Random Forest')
plt.xticks(xplot, df_39.index.values.tolist())
plt.scatter(xplot, df_39['TOUS.CANCERS_2005.09'].tolist(), color='red', label='Valeur attendue')
plt.plot(xplot, y_pred2.T[0].tolist(), color='blue', label='Valeur predite')
plt.legend()
plt.show()
print(metrics.mean_squared_error(df_13['TOUS.CANCERS_2005.09'].tolist(), np.array(y_pred2)))
Par comparaison des mean squared errors, nous trouvons un meilleur résultat pour la régression multilinéaire.
De même, nous recherchons les hyperparamètres par GridSearch.
from sklearn.ensemble import GradientBoostingRegressor
param_grid = {'learning_rate': [0.1, 0.05, 0.02, 0.01],
'max_depth': [4, 6],
'min_samples_leaf': [3, 5, 9, 17],
'max_features': [1.0, 0.3, 0.1]
}
for train_index, test_index in kf2:
clf4 = GridSearchCV(GradientBoostingRegressor(n_estimators=3000), param_grid, cv=5,
scoring='mean_squared_error')
X_train, X_test = X_tot1[train_index], X_tot1[test_index]
y_train, y_test = Y_tot1[train_index], Y_tot1[test_index]
clf4.fit(X_train, y_train)
print("Best parameters set found on development set:")
print()
print(clf.best_params_)
print()
print("Grid scores on development set:")
print()
for params, mean_score, scores in clf4.grid_scores_:
print("%0.3f (+/-%0.03f) for %r"
% (mean_score, scores.std() * 2, params))
print()
print("Detailed classification report:")
print()
print("The model is trained on the full development set.")
print("The scores are computed on the full evaluation set.")
print()
y_true, y_pred = y_test, clf4.predict(X_test)
Nous prenons {'max_features': 0.3, 'learning_rate': 0.02, 'max_depth': 6, 'min_samples_leaf': 17}
En plus des prédictions, nous ajoutons les intervalles de prédiction à 90%.
tuned_params = {'n_estimators' : 3000,'max_features': 0.3, 'learning_rate': 0.02, 'max_depth': 6, 'min_samples_leaf': 17}
listy_pred3=pd.DataFrame(data = np.zeros((df_39.shape[0],1)), index = df_39.index)
listy_predupper=pd.DataFrame(data = np.zeros((df_39.shape[0],1)), index = df_39.index)
listy_predlower=pd.DataFrame(data = np.zeros((df_39.shape[0],1)), index = df_39.index)
# Cadrille l'espace afin de tracer la courbe
xx = np.atleast_2d(np.linspace(1, 96, 96)).T
xx = xx.astype(np.float32)
clf = GradientBoostingRegressor(**tuned_params)
for train_index, test_index in kf2:
#Découpe la base
X_train, X_test = X_tot1[train_index], X_tot1[test_index]
y_train = Y_tot1[train_index]
#On initialise le régresseur
alpha = 0.95
clf.set_params(loss='quantile', alpha=alpha)
clf.fit(X_train, y_train)
#On prédit la borne supérieure de l'intervalle
temp_df = pd.DataFrame(data = np.zeros((df_39.shape[0],1)), index = df_39.index)
temp_df.ix[test_index] = clf.predict(X_test).reshape(32,1)
listy_predupper.ix[test_index] += temp_df.ix[test_index]
#On réinitialise
clf.set_params(alpha=1.0 - alpha)
clf.fit(X_train, y_train)
#On prédit la borne inférieure de l'intervalle
temp_df = pd.DataFrame(data = np.zeros((df_39.shape[0],1)), index = df_39.index)
temp_df.ix[test_index] = clf.predict(X_test).reshape(32,1)
listy_predlower.ix[test_index] += temp_df.ix[test_index]
#On réinitialise
clf.set_params(loss='ls')
clf.fit(X_train, y_train)
# On prédit la valeur.
temp_df = pd.DataFrame(data = np.zeros((df_39.shape[0],1)), index = df_39.index)
temp_df.ix[test_index] = clf.predict(X_test).reshape(32,1)
listy_pred3.ix[test_index] += temp_df.ix[test_index]
y_upper = listy_predupper.values
y_lower = listy_predlower.values
y_pred3 = listy_pred3.values
xplot = [i for i in range(y_pred3.shape[0])]
plt.figure(figsize=(18,6))
plt.title('Comparaison entre la valeur attendue et la valeur predite avec le modele du Gradient Boosting')
plt.xticks(xplot, df_39.index.values.tolist())
plt.scatter(xplot, df_39['TOUS.CANCERS_2005.09'].tolist(), color='red', label='Valeur attendue')
plt.plot(xplot, y_pred3.T[0].tolist(), color='blue', label='Valeur predite')
plt.legend()
plt.show()
print(metrics.mean_squared_error(Y_tot1, y_pred3))
# Nous traçons les observations et leur prédiction avec un intervalle de prédiction de 90%
fig = plt.figure(figsize=(18,6))
plt.plot(xx, Y_tot1, 'b.', markersize=10, label=u'Observations')
plt.plot(xx, y_pred3, 'r-', label=u'Prediction')
plt.plot(xx, y_upper, 'k-')
plt.plot(xx, y_lower, 'k-')
plt.fill(np.concatenate([xx, xx[::-1]]),
np.concatenate([y_upper, y_lower[::-1]]),
alpha=.5, fc='b', ec='None', label='intervalle de prediction a 90%')
plt.xticks(xx,df_39.index.values.tolist())
plt.xlabel('Departements')
plt.ylabel('Mortalite tous cancers 2005-2009')
plt.legend(loc='upper left')
plt.show()
xplot = [i for i in range(df_13.shape[0])]
plt.figure(figsize=(18,6))
plt.title('Synthese generale des performances de nos modeles')
plt.xticks(xplot, df_13.index.values.tolist())
plt.scatter(xplot, df_13['TOUS.CANCERS_2005.09'].tolist(), color='red', label='Valeur attendue')
plt.plot(xplot, y_pred1, color='blue', label='Regression multilineaire')
plt.plot(xplot, y_pred2, color='green', label='Random Forest')
plt.plot(xplot, y_pred3, color='magenta', label='Gradient Boosting')
plt.legend()
plt.show()
Il a été assez surprenant de constater que les deux modèles les plus simples utilisés (le Random Forest et la Régression multilinéaire) ont été les modèles les plus efficaces dans notre projet.
Ainsi ces deux modèles présentaient une erreur quadratique moyenne significativement inférieure au Gradient Boosting. Nous les préférons donc à lui.
La question se pose maintenant du modèle à retenir entre le Random Forest et la régression multilinéaire. Les erreurs quadratiques moyennes étant proches, nous étudions les graphes des modèles. De cette étude, il en ressort que le Random Forest est préféré, car nous n'avons pas dans notre projet épuisé toutes les possibilités offertes par l'hyperparamétrisation du Random Forest. En outre, le Random Forest prédit généralement mieux les départements normaux, mais a plus de difficulté à traiter les départements atypiques. Sachant cela, nous pouvons procéder à une étude spécifique à chaque outlier détecté et se fier aux prédictions du Random Forest pour les autres départements.
Les erreurs restent tout de même élevées ; de ce fait, on peut penser qu'un travail sur les hyperparamètres ne suffise pas. D'autres modèles et d'autres données (non agrégées par exemple) sont probablement à considérer. Toutefois, il semble assez difficile, encore aujourd'hui, d'obtenir des données individuelles sur la santé ; la question éthique est encore cruciale dans ce domaine et peut constituer un frein (dont nous ne discutons pas ici le caractère salutaire) pour certaines recherches.
Néanmoins, notre but initial dans le cadre du challenge n'est pas totalement atteint : nous souhaitions en effet également exhiber des mécanismes de pollution qui seraient (selon notre étude) cancérigènes. Il faudrait alors dans un second temps porter une plus grande attention à ces mécanismes, notamment par le biais des coefficients dans les régressions et de leur significativité. De plus, les différents niveaux de détails restent à être davantage exploités. Ainsi, il est possible de nuancer les résultats, notamment par sexe ou par type de cancers.
Nous concluons donc ici le rapport de notre projet de Machine Learning, mais le challenge, lui, continue. Le début de l'année 2016 sera l'occasion d'approfondir nos recherches.
Benoît Choffin - Joseph Lam