Key Word(s): Decision Tree, Bootstrap-Aggregating (Bagging), Random Forest, Boosting, Adaboost
CS-109A Introduction to Data Science
Lab 9: Decision Trees, Bagged Trees, Random Forests and Boosting - Student Version¶
Harvard University
Fall 2018
Instructors: Pavlos Protopapas and Kevin Rader
Lab Instructor: Kevin Rader (today at least)
Authors: Kevin Rader, Rahul Dave
## RUN THIS CELL TO PROPERLY HIGHLIGHT THE EXERCISES
import requests
from IPython.core.display import HTML
styles = requests.get("https://raw.githubusercontent.com/Harvard-IACS/2018-CS109A/master/content/styles/cs109.css").text
HTML(styles)
We will look here into the practicalities of fitting regression trees, random forests, and boosted trees. These involve out-of-bound estmates and cross-validation, and how you might want to deal with hyperparameters in these models. Along the way we will play a little bit with different loss functions, so that you start thinking about what goes in general into cooking up a machine learning model.
%matplotlib inline
import numpy as np
import scipy as sp
import matplotlib as mpl
import matplotlib.cm as cm
import matplotlib.pyplot as plt
import pandas as pd
pd.set_option('display.width', 500)
pd.set_option('display.max_columns', 100)
pd.set_option('display.notebook_repr_html', True)
import seaborn as sns
sns.set_style("whitegrid")
sns.set_context("poster")
Dataset¶
First, the data. We will be attempting to predict the presidential election results (at the county level) from 2016, measured as 'votergap' = (trump - clinton) in percentage points, based mostly on demographic features of those counties. Let's quick take a peak at the data:
elect_df = pd.read_csv("county_level_election.csv")
elect_df.head()
from sklearn.model_selection import train_test_split
# split 80/20 train-test
X = elect_df[['population','hispanic','minority','female','unemployed','income','nodegree','bachelor','inactivity','obesity','density','cancer']]
response = elect_df['votergap']
Xtrain, Xtest, ytrain, ytest = train_test_split(X,response,test_size=0.2)
plt.hist(ytrain)
Xtrain.hist(column=['minority', 'population','hispanic','female']);
How would you describe these variables?
print(elect_df.shape)
print(Xtrain.shape)
print(Xtest.shape)
Regression Trees¶
We will start by using a simple Decision Tree Regressor to predict votergap. Thats not the aim of this lab, so we'll run a few of these models without any cross-validation or 'regularization' just to illustrate what is going on.
This is what you ought to keep in mind about decision trees.
from the docs:
max_depth : int or None, optional (default=None)
The maximum depth of the tree. If None, then nodes are expanded until all leaves are pure or until all leaves contain less than min_samples_split samples.
min_samples_split : int, float, optional (default=2)- The deeper the tree, the more prone you are to overfitting.
- The smaller
min_samples_split, the more the overfitting. One may usemin_samples_leafinstead. More samples per leaf, the higher the bias.
from sklearn.tree import DecisionTreeRegressor
#x = np.arange(0, 2*np.pi, 0.1)
#y = np.sin(x) + 0.1*np.random.normal(size=x.shape[0])
x = Xtrain['minority'].values
o = np.argsort(x)
x = x[o]
y = ytrain.values
y = y[o]
plt.plot(x,y, '.');
plt.plot(np.log(x),y, '.');
plt.plot(np.log(x),y,'.')
xx = np.log(x).reshape(-1,1)
for i in [1,2]:
dtree = DecisionTreeRegressor(max_depth=i)
dtree.fit(xx, y)
plt.plot(np.log(x), dtree.predict(xx), label=str(i), alpha=1-i/10, lw=4)
plt.legend();
plt.plot(np.log(x),y,'.')
xx = np.log(x).reshape(-1,1)
for i in [500,200,100,20]:
dtree = DecisionTreeRegressor(min_samples_split=i)
dtree.fit(xx, y)
plt.plot(np.log(x), dtree.predict(xx), label=str(i), alpha=0.8, lw=4)
plt.legend();
plt.plot(np.log(x),y,'.')
xx = np.log(x).reshape(-1,1)
for i in [500,200,100,20]:
dtree = DecisionTreeRegressor(max_depth=6, min_samples_split=i)
dtree.fit(xx, y)
plt.plot(np.log(x), dtree.predict(xx), label=str(i), alpha=0.8, lw=4)
plt.legend();
#let's also include logminority as a predictor going forward
xtemp = np.log(Xtrain['minority'].values)
Xtrain = Xtrain.assign(logminority = xtemp)
Xtest = Xtest.assign(logminority = np.log(Xtest['minority'].values))
Xtrain.head()
- Perform 5-fold cross-validation to determine what the best
max_depthwould be for a single regression tree using the entire 'Xtrain' feature set. - Visualize the results with mean +/- 2 sd's across the validation sets.
from sklearn.model_selection import cross_val_score
# your code here
depths = list(range(1, 21))
train_scores = []
cvmeans = []
cvstds = []
cv_scores = []
for depth in depths:
dtree = DecisionTreeRegressor(max_depth=depth)
# Perform 5-fold cross validation and store results
cv_scores = cross_val_score(dtree, Xtrain,ytrain,cv=5)
cvmeans = append()
cvmeans = np.array(cvmeans)
cvstds = np.array(cvstds)
# your code here
# plot means and shade the 2 SD interval
#plt.fill_between(depths, cvmeans - 2*cvstds, cvmeans + 2*cvstds, alpha=0.3);
Ok with this discussion in mind, lets improve this model by Bagging.
Bootstrap-Aggregating (called Bagging)¶
Whats the basic idea?
- A Single Decision tree is likely to overfit.
- So lets introduce replication through Bootstrap sampling.
- Bagging uses bootstrap resampling to create different training datasets. This way each training will give us a different tree.
- Added bonus: the left off points can be used to as a natural "validation" set, so no need to
- Since we have many trees that we will average over for prediction, we can choose a large
max_depthand we are ok as we will rely on the law of large numbers to shrink this large variance, low bias approach for each individual tree.
from sklearn.utils import resample
ntrees = 500
estimators = []
R2s = []
yhats_test = np.zeros((Xtest.shape[0], ntrees))
plt.plot(np.log(x),y,'.')
for i in range(ntrees):
simpletree = DecisionTreeRegressor(max_depth=3)
boot_xx, boot_y = resample(Xtrain[['logminority']], ytrain)
estimators = np.append(estimators,simpletree.fit(boot_xx, boot_y))
R2s = np.append(R2s,simpletree.score(Xtest[['logminority']], ytest))
yhats_test[:,i] = simpletree.predict(Xtest[['logminority']])
plt.plot(np.log(x), simpletree.predict(np.log(x).reshape(-1,1)), 'red', alpha=0.05)
yhats_test.shape
- Edit the code below (which is just copied from above) to refit many bagged trees on the entire xtrain feature set (without the plot...lots of predictors now so difficult to plot).
- Summarize how each of the separate trees performed (both numerically and visually) using $R^2$ as the metric. How do they perform on average?
- Combine the trees into one prediction and evaluate it using $R^2$.
- Briefly discuss the results. How will the results above change if 'max_depth=4' is increased? What if it is decreased?
from sklearn.metrics import r2_score
ntrees = 500
estimators = []
R2s = []
yhats_test = np.zeros((Xtest.shape[0], ntrees))
for i in range(ntrees):
dtree = DecisionTreeRegressor(max_depth=3)
boot_xx, boot_y = resample(Xtrain[['logminority']], ytrain)
estimators = np.append(estimators,dtree.fit(boot_xx, boot_y))
R2s = np.append(R2s,dtree.score(Xtest[['logminority']], ytest))
yhats_test[:,i] = dtree.predict(Xtest[['logminority']])
# your code here
Your answer here¶
Random Forests¶
What's the basic idea?
Bagging alone is not enough randomization, because even after bootstrapping, we are mainly training on the same data points using the same variablesn, and will retain much of the overfitting.
So we will build each tree by splitting on "random" subset of predictors at each split (hence, each is a 'random tree'). This can't be done in with just one predcitor, but with more predictors we can choose what predictors to split on randomly and how many to do this on. Then we combine many 'random trees' together by averaging their predictions, and this gets us a forest of random trees: a random forest.
Below we create a hyper-param Grid. We are preparing to use the bootstrap points not used in training for validation.
max_features : int, float, string or None, optional (default=”auto”)
- The number of features to consider when looking for the best split.max_features: Default splits on all the features and is probably prone to overfitting. You'll want to validate on this.- You can "validate" on the trees
n_estimatorsas well but many a times you will just look for the plateau in the trees as seen below. - From decision trees you get the
max_depth,min_samples_split, andmin_samples_leafas well but you might as well leave those at defaults to get a maximally expanded tree.
from sklearn.ensemble import RandomForestRegressor
# code from
# Adventures in scikit-learn's Random Forest by Gregory Saunders
from itertools import product
from collections import OrderedDict
param_dict = OrderedDict(
n_estimators = [400, 600, 800],
max_features = [0.2, 0.4, 0.6, 0.8]
)
param_dict.values()
Using the OOB score.¶
We have been putting "validate" in quotes. This is because the bootstrap gives us left-over points! So we'll now engage in our very own version of a grid-search, done over the out-of-bag scores that sklearn gives us for free
from itertools import product
#make sure ytrain is the correct data type...in case you have warnings
#print(yytrain.shape,ytrain.shape,Xtrain.shape)
#ytrain = np.ravel(ytrain)
#Let's Cross-val. on the two 'hyperparameters' we based our grid on earlier
results = {}
estimators= {}
for ntrees, maxf in product(*param_dict.values()):
params = (ntrees, maxf)
est = RandomForestRegressor(oob_score=True,
n_estimators=ntrees, max_features=maxf, max_depth=50, n_jobs=-1)
est.fit(Xtrain, ytrain)
results[params] = est.oob_score_
estimators[params] = est
outparams = max(results, key = results.get)
outparams
rf1 = estimators[outparams]
results
rf1.score(Xtest, ytest)
Finally you can find the feature importance of each predictor in this random forest model. Whenever a feature is used in a tree in the forest, the algorithm will log the decrease in the splitting criterion (such as gini). This is accumulated over all trees and reported in est.feature_importances_
pd.Series(rf1.feature_importances_,index=list(Xtrain)).sort_values().plot(kind="barh")
Since our response isn't very symmetric, we may want to suppress outliers by using the mean_absolute_error instead.
from sklearn.metrics import mean_absolute_error
mean_absolute_error(ytest, rf1.predict(Xtest))
- Tune the 'RandomForestRegressor' above to minimize 'mean_absolute_error' instead of the default __
Note: sklearn supports this (criterion='mae') since v0.18, but does not have completely arbitrary loss functions for Random Forests.
# your code here
Note: instead of using oob scoring, we could do cross-validation, and a cv of 3 will roughly be comparable (same approximate train-vs.-validation set sizes). But this will take much more time as you are doing the fit 3 times plus the bootstraps. So at least three times as long!
param_dict2 = OrderedDict(
n_estimators = [600],
max_features = [0.2, 0.4, 0.6, 0.8]
)
from sklearn.model_selection import GridSearchCV
est2 = RandomForestRegressor(oob_score=False)
gs = GridSearchCV(est2, param_grid = param_dict2, cv=3, n_jobs=-1)
gs.fit(Xtrain, ytrain)
rf2 = gs.best_estimator_
rf2
gs.best_score_
- What are the 3 hyperparameters for a random forest (one of the hyperparameters come in many flavors)?
- Which hyperparameters need to be tuned?
- How would you go about tuning these hyperparemeters?
Your answer here¶
Seeing error as a function of the proportion of predictors at each split¶
Let's look to see how max_features affects performance on the training set.
# from http://scikit-learn.org/stable/auto_examples/ensemble/plot_ensemble_oob.html
feats = param_dict['max_features']
#
error_rate = OrderedDict((label, []) for label in feats)
# Range of `n_estimators` values to explore.
min_estimators = 200
step_estimators = 100
num_steps = 3
max_estimators = min_estimators + step_estimators*num_steps
for label in feats:
for i in range(min_estimators, max_estimators+1, step_estimators):
clf = RandomForestRegressor(oob_score=True, max_features=label)
clf.set_params(n_estimators=i)
clf.fit(Xtrain, ytrain)
# Record the OOB error for each `n_estimators=i` setting.
oob_error = 1 - clf.oob_score_
error_rate[label].append((i, oob_error))
print(feats)
# Generate the "OOB error rate" vs. "n_estimators" plot.
for label, clf_err in error_rate.items():
xs, ys = zip(*clf_err)
plt.plot(xs, ys, label=label)
plt.xlim(min_estimators, max_estimators)
plt.xlabel("n_estimators")
plt.ylabel("OOB error rate")
plt.legend(loc="upper right")
plt.show()
Boosting Regression Trees¶
Adaboost Classification, which you will be doing in your homework, is a special case of a gradient-boosted algorithm. Gradient Bossting is very state of the art, and has major connections to logistic regression, gradient descent in a functional space, and search in information space. See Schapire and Freund's MIT Press book for details (Google is a wonderful thing).
But briefly, let us cover the idea here. The idea is that we will use a bunch of weak 'learners' (aka, models) which are fit sequentially. The first one fits the signal, the second one the first model's residual, the third the second model's residual, and so on. At each stage we upweight the places that our previous model did badly on. First let us illustrate.
from sklearn.ensemble import AdaBoostRegressor
estab = AdaBoostRegressor(base_estimator=DecisionTreeRegressor(max_depth=1), n_estimators=200, learning_rate=1.0)
estab.fit(xx, y)
staged_predict_generator = estab.staged_predict(xx)
# code from http://nbviewer.jupyter.org/github/pprett/pydata-gbrt-tutorial/blob/master/gbrt-tutorial.ipynb
import time
from IPython import display
plt.plot(xx, y, '.');
counter = 0
for stagepred in staged_predict_generator:
counter = counter + 1
if counter in [1, 2, 4, 8, 10, 50, 100, 200]:
plt.plot(xx, stagepred, alpha=0.7, label=str(counter), lw=4)
plt.legend();
display.display(plt.gcf())
display.clear_output(wait=True)
time.sleep(1)
print(counter, i)
Ok, so this demonstration helps us understand some things about boosting.
n_estimatorsis the number of trees, and thus the stage in the fitting. It also controls the complexity for us. The more trees we have the more we fit to the tiny details.staged_predictgives us the prediction at each step- once again
max_depthfrom the underlying decision tree tells us the depth of the tree. But here it tells us the amount of features interactions we have, not just the scale of our fit. But clearly it increases the variance again.
Ideas from decision trees remain. For example, increase min_samples_leaf or decrease max_depth to reduce variance and increase the bias.
- What do you expect to happen if you increase
max_depthto 5? Edit the code above to explore the result. - What do you expect to happen if you put
max_depthback to 1 and decrease thelearning_rateto 0.1? Edit the code above to explore the result. - Do a little work to find some sort of 'best' values of
max_depthandlearning_rate. Does this result make sense?
# your code here
Note: what's the relationship between residuals and the gradient?¶
Pavlos showed in class that for the squared loss, taking the gradient in the "data point functional space", ie a N-d space for N data points with each variable being $f(x_i)$ just gives us the residuals. It turns out that the gradient descent is a more general idea, and one can use this for different losses. And the upweighting of poorly fit points in AdaBoost is simply a weighing by gradient. If the gradient (or residual) is high it means you are far away from optimum in this functional space, and if you are at 0, you have a flat gradient!
The ideas from the general theory of gradient descent tell us this: we can slow the learning by shrinking the predictions of each tree by some small number, which is called the learning_rate (learning_rate). This "shrinkage" helps us not overshoot, but for a finite number of iterations also simultaneously ensures we dont overfit by being in the neighboorhood of the minimum rather than just at it! But we might need to increase the iterations some to get into the minimum area.