import pandas as pd # pour importer le jeu de données
import numpy as np  # pour les simple manips de base
from sklearn.ensemble import RandomForestClassifier, RandomForestRegressor # pour faire des RF
from sklearn.model_selection import GridSearchCV # pour faire de la validation croisée
import matplotlib.pyplot as plt # faire une figure ou deux 
%matplotlib inline

prostate = pd.read_csv('../data/Prostate_Cancer.txt', delimiter=',')
prostate.head()

prostate = prostate.drop(['id'], axis=1)
prostate.shape

(97, 10)

Jeux de données apprentissage et test¶

# binariser gleason 
prostate['gleason'][prostate['gleason']==6] = 0
prostate['gleason'][prostate['gleason']>6] = 1

# Train-test split
y_train = np.array(prostate[prostate.train == "T"]['lpsa'])
y_test = np.array(prostate[prostate.train == "F"]['lpsa'])
X_train = np.array(prostate[prostate.train == "T"].drop(['lpsa', 'train'], axis=1))
X_test = np.array(prostate[prostate.train == "F"].drop(['lpsa', 'train'], axis=1))

Note favorable à python (d'après Philippe Besse)¶

La librairie randomForest de R utilise le programme historique développé par Breiman et Cutler (2001) et interfacé par Liaw et Wiener. Cette interface est toujours mise à jour mais il n'est pas sûr que le programme original continue d'évoluer depuis 2004. Pour des tailles importantes d'échantillons, quelques milliers, cette implémentation atteint des temps d'exécution rédhibitoires (cf. cet exemple) au contraire de celle en Python dont gestion mémoire et capacité de parallélisation ont été finement optimisées par Louppe et al. (2014).

De même que le boosting, deux fonctions de forêt sont proposés dans scikit-learn; une pour la régression et une pour la classification ainsi qu'une version "plus aléatoire". Par rapport à la version originale de R, moins d'options sont proposées mais l'utilisation de base est très similaire avec le même jeu de paramètres.

forest = RandomForestRegressor(n_estimators=500, 
   criterion='mse', max_depth=None,
   min_samples_split=2, min_samples_leaf=1, 
   max_features='auto', max_leaf_nodes=None,
   bootstrap=True, oob_score=True)
# apprentissage
rfFit = forest.fit(X_train, y_train)
## affichage du RMSE calculé par OOB
print(np.sqrt(rfFit.oob_score_))

0.6739066604363924

Comparer l'erreur out-of-bag ci-dessus avec celle sur l'échantillon test.¶

# erreur de prédiction du jeu de données test
print(np.sqrt(rfFit.score(X_test,y_test)))

0.7869590276207928

Optimisation par validation croisée du nombre de variables tirés aléatoirement lors de la construction de chaque noeud.¶

Le résultat optimal est-il stable ?¶

Plusieurs exécutions, rendues aléatoires par la validation croisée, peuvent conduire à des valeurs "optimales" différentes de ce paramètre sans pour autant nuire à la qualité de prévision sur l'échantillon test.

param=[{"max_features":list(range(2,8,1))}]
rf= GridSearchCV(RandomForestRegressor(n_estimators=500),
        param,cv=3,n_jobs=-1)
rfOpt=rf.fit(X_train, y_train)
# paramètre optimal
rfOpt.best_params_['max_features']

2

nb_var=rfOpt.best_params_['max_features']
rf = RandomForestRegressor(n_estimators=500,
   criterion='mse', max_depth=None,
   min_samples_split=2, min_samples_leaf=1, 
   max_features=nb_var,max_leaf_nodes=None,
   bootstrap=True, oob_score=True)
# apprentissage
rfFit = rf.fit(X_train,y_train)
# erreur de prévision sur le test
np.sqrt(rfFit.score(X_test,y_test))

0.7017630764560215

Comme avec R, il est possible de calculer un indicateur d'importance des variables pour aider à une forme d'interprétation. Celui-ci dépend de la position de la variable dans l'arbre et correspond donc au mean decrease in mse.

# Importance décroissante des variables
importances = rf.feature_importances_
indices = np.argsort(importances)[::-1]
for f in range(X_train.shape[1]):
    print(prostate.columns[indices[f]], importances[indices[f]])

lcavol 0.28732125195976776
lweight 0.2055794048229232
pgg45 0.10193233790173001
age 0.09495552160844897
lcp 0.0939625261274299
svi 0.0832720708938718
lbph 0.08288603835841744
gleason 0.05009084832741083

# Graphe des importances
plt.figure()
plt.title("Importances des variables")
plt.bar(range(X_train.shape[1]), importances[indices])
plt.xticks(range(X_train.shape[1]), indices)
plt.xlim([-1, X_train.shape[1]])
plt.show()

	id	lcavol	lweight	age	lbph	lcp	gleason	pgg45	lpsa	train
0	1	-0.579818	2.769459	50	-1.386294	-1.386294	6	0	-0.430783	T
1	2	-0.994252	3.319626	58	-1.386294	-1.386294	6	0	-0.162519	T
2	3	-0.510826	2.691243	74	-1.386294	-1.386294	7	20	-0.162519	T
3	4	-1.203973	3.282789	58	-1.386294	-1.386294	6	0	-0.162519	T
4	5	0.751416	3.432373	62	-1.386294	-1.386294	6	0	0.371564	T