#!/usr/bin/env python # coding: utf-8 # # 07 - Ensemble Methods # # by [Alejandro Correa Bahnsen](http://www.albahnsen.com/) & [Iván Torroledo](http://www.ivantorroledo.com/) # # version 1.3, June 2018 # # ## Part of the class [Applied Deep Learning](https://github.com/albahnsen/AppliedDeepLearningClass) # # # This notebook is licensed under a [Creative Commons Attribution-ShareAlike 3.0 Unported License](http://creativecommons.org/licenses/by-sa/3.0/deed.en_US). Special thanks goes to [Kevin Markham](https://github.com/justmarkham) # Why are we learning about ensembling? # # - Very popular method for improving the predictive performance of machine learning models # - Provides a foundation for understanding more sophisticated models # ## Lesson objectives # # Students will be able to: # # - Define ensembling and its requirements # - Identify the two basic methods of ensembling # - Decide whether manual ensembling is a useful approach for a given problem # - Explain bagging and how it can be applied to decision trees # - Explain how out-of-bag error and feature importances are calculated from bagged trees # - Explain the difference between bagged trees and Random Forests # - Build and tune a Random Forest model in scikit-learn # - Decide whether a decision tree or a Random Forest is a better model for a given problem # # Part 1: Introduction # Ensemble learning is a widely studied topic in the machine learning community. The main idea behind # the ensemble methodology is to combine several individual base classifiers in order to have a # classifier that outperforms each of them. # # Nowadays, ensemble methods are one # of the most popular and well studied machine learning techniques, and it can be # noted that since 2009 all the first-place and second-place winners of the KDD-Cup https://www.sigkdd.org/kddcup/ used ensemble methods. The core # principle in ensemble learning, is to induce random perturbations into the learning procedure in # order to produce several different base classifiers from a single training set, then combining the # base classifiers in order to make the final prediction. In order to induce the random permutations # and therefore create the different base classifiers, several methods have been proposed, in # particular: # * bagging # * pasting # * random forests # * random patches # # Finally, after the base classifiers # are trained, they are typically combined using either: # * majority voting # * weighted voting # * stacking # # There are three main reasons regarding why ensemble # methods perform better than single models: statistical, computational and representational . First, from a statistical point of view, when the learning set is too # small, an algorithm can find several good models within the search space, that arise to the same # performance on the training set $\mathcal{S}$. Nevertheless, without a validation set, there is # a risk of choosing the wrong model. The second reason is computational; in general, algorithms # rely on some local search optimization and may get stuck in a local optima. Then, an ensemble may # solve this by focusing different algorithms to different spaces across the training set. The last # reason is representational. In most cases, for a learning set of finite size, the true function # $f$ cannot be represented by any of the candidate models. By combining several models in an # ensemble, it may be possible to obtain a model with a larger coverage across the space of # representable functions. # ![s](images/ch9_fig1.png) # ## Example # # Let's pretend that instead of building a single model to solve a binary classification problem, you created **five independent models**, and each model was correct about 70% of the time. If you combined these models into an "ensemble" and used their majority vote as a prediction, how often would the ensemble be correct? # In[1]: import numpy as np # set a seed for reproducibility np.random.seed(1234) # generate 1000 random numbers (between 0 and 1) for each model, representing 1000 observations mod1 = np.random.rand(1000) mod2 = np.random.rand(1000) mod3 = np.random.rand(1000) mod4 = np.random.rand(1000) mod5 = np.random.rand(1000) # each model independently predicts 1 (the "correct response") if random number was at least 0.3 preds1 = np.where(mod1 > 0.3, 1, 0) preds2 = np.where(mod2 > 0.3, 1, 0) preds3 = np.where(mod3 > 0.3, 1, 0) preds4 = np.where(mod4 > 0.3, 1, 0) preds5 = np.where(mod5 > 0.3, 1, 0) # print the first 20 predictions from each model print(preds1[:20]) print(preds2[:20]) print(preds3[:20]) print(preds4[:20]) print(preds5[:20]) # In[2]: # average the predictions and then round to 0 or 1 ensemble_preds = np.round((preds1 + preds2 + preds3 + preds4 + preds5)/5.0).astype(int) # print the ensemble's first 20 predictions print(ensemble_preds[:20]) # In[3]: # how accurate was each individual model? print(preds1.mean()) print(preds2.mean()) print(preds3.mean()) print(preds4.mean()) print(preds5.mean()) # In[4]: # how accurate was the ensemble? print(ensemble_preds.mean()) # **Note:** As you add more models to the voting process, the probability of error decreases, which is known as [Condorcet's Jury Theorem](http://en.wikipedia.org/wiki/Condorcet%27s_jury_theorem). # ## What is ensembling? # # **Ensemble learning (or "ensembling")** is the process of combining several predictive models in order to produce a combined model that is more accurate than any individual model. # # - **Regression:** take the average of the predictions # - **Classification:** take a vote and use the most common prediction, or take the average of the predicted probabilities # # For ensembling to work well, the models must have the following characteristics: # # - **Accurate:** they outperform the null model # - **Independent:** their predictions are generated using different processes # # **The big idea:** If you have a collection of individually imperfect (and independent) models, the "one-off" mistakes made by each model are probably not going to be made by the rest of the models, and thus the mistakes will be discarded when averaging the models. # # There are two basic **methods for ensembling:** # # - Manually ensemble your individual models # - Use a model that ensembles for you # ### Theoretical performance of an ensemble # If we assume that each one of the $T$ base classifiers has a probability $\rho$ of # being correct, the probability of an ensemble making the correct decision, assuming independence, # denoted by $P_c$, can be calculated using the binomial distribution # # $$P_c = \sum_{j>T/2}^{T} {{T}\choose{j}} \rho^j(1-\rho)^{T-j}.$$ # # Furthermore, as shown, if $T\ge3$ then: # # $$# \lim_{T \to \infty} P_c= \begin{cases} # 1 &\mbox{if } \rho>0.5 \\ # 0 &\mbox{if } \rho<0.5 \\ # 0.5 &\mbox{if } \rho=0.5 , # \end{cases} #$$ # leading to the conclusion that # $$# \rho \ge 0.5 \quad \text{and} \quad T\ge3 \quad \Rightarrow \quad P_c\ge \rho. #$$ # # Part 2: Manual ensembling # # What makes a good manual ensemble? # # - Different types of **models** # - Different combinations of **features** # - Different **tuning parameters** # ![Machine learning flowchart](https://raw.githubusercontent.com/justmarkham/DAT8/master/notebooks/images/crowdflower_ensembling.jpg) # # *Machine learning flowchart created by the [winner](https://github.com/ChenglongChen/Kaggle_CrowdFlower) of Kaggle's [CrowdFlower competition](https://www.kaggle.com/c/crowdflower-search-relevance)* # In[5]: # read in and prepare the vehicle training data import zipfile import pandas as pd with zipfile.ZipFile('../datasets/vehicles_train.csv.zip', 'r') as z: f = z.open('vehicles_train.csv') train = pd.io.parsers.read_table(f, index_col=False, sep=',') with zipfile.ZipFile('../datasets/vehicles_test.csv.zip', 'r') as z: f = z.open('vehicles_test.csv') test = pd.io.parsers.read_table(f, index_col=False, sep=',') train['vtype'] = train.vtype.map({'car':0, 'truck':1}) # read in and prepare the vehicle testing data test['vtype'] = test.vtype.map({'car':0, 'truck':1}) # In[6]: train.head() # ### Train different models # In[7]: from sklearn.linear_model import LinearRegression from sklearn.tree import DecisionTreeRegressor from sklearn.naive_bayes import GaussianNB from sklearn.neighbors import KNeighborsRegressor models = {'lr': LinearRegression(), 'dt': DecisionTreeRegressor(), 'nb': GaussianNB(), 'nn': KNeighborsRegressor()} # In[8]: # Train all the models X_train = train.iloc[:, 1:] X_test = test.iloc[:, 1:] y_train = train.price y_test = test.price for model in models.keys(): models[model].fit(X_train, y_train) # In[9]: # predict test for each model y_pred = pd.DataFrame(index=test.index, columns=models.keys()) for model in models.keys(): y_pred[model] = models[model].predict(X_test) # In[10]: # Evaluate each model from sklearn.metrics import mean_squared_error for model in models.keys(): print(model,np.sqrt(mean_squared_error(y_pred[model], y_test))) # ### Evaluate the error of the mean of the predictions # In[11]: np.sqrt(mean_squared_error(y_pred.mean(axis=1), y_test)) # ## Comparing manual ensembling with a single model approach # # **Advantages of manual ensembling:** # # - Increases predictive accuracy # - Easy to get started # # **Disadvantages of manual ensembling:** # # - Decreases interpretability # - Takes longer to train # - Takes longer to predict # - More complex to automate and maintain # - Small gains in accuracy may not be worth the added complexity # # Part 3: Bagging # # The primary weakness of **decision trees** is that they don't tend to have the best predictive accuracy. This is partially due to **high variance**, meaning that different splits in the training data can lead to very different trees. # # **Bagging** is a general purpose procedure for reducing the variance of a machine learning method, but is particularly useful for decision trees. Bagging is short for **bootstrap aggregation**, meaning the aggregation of bootstrap samples. # # What is a **bootstrap sample**? A random sample with replacement: # In[12]: # set a seed for reproducibility np.random.seed(1) # create an array of 1 through 20 nums = np.arange(1, 21) print(nums) # sample that array 20 times with replacement print(np.random.choice(a=nums, size=20, replace=True)) # **How does bagging work (for decision trees)?** # # 1. Grow B trees using B bootstrap samples from the training data. # 2. Train each tree on its bootstrap sample and make predictions. # 3. Combine the predictions: # - Average the predictions for **regression trees** # - Take a vote for **classification trees** # # Notes: # # - **Each bootstrap sample** should be the same size as the original training set. # - **B** should be a large enough value that the error seems to have "stabilized". # - The trees are **grown deep** so that they have low bias/high variance. # # Bagging increases predictive accuracy by **reducing the variance**, similar to how cross-validation reduces the variance associated with train/test split (for estimating out-of-sample error) by splitting many times an averaging the results. # # In[13]: # set a seed for reproducibility np.random.seed(123) n_samples = train.shape[0] n_B = 10 # create ten bootstrap samples (will be used to select rows from the DataFrame) samples = [np.random.choice(a=n_samples, size=n_samples, replace=True) for _ in range(1, n_B +1 )] samples # In[14]: # show the rows for the first decision tree train.iloc[samples[0], :] # Build one tree for each sample # In[15]: from sklearn.tree import DecisionTreeRegressor # grow each tree deep treereg = DecisionTreeRegressor(max_depth=None, random_state=123) # DataFrame for storing predicted price from each tree y_pred = pd.DataFrame(index=test.index, columns=[list(range(n_B))]) # grow one tree for each bootstrap sample and make predictions on testing data for i, sample in enumerate(samples): X_train = train.iloc[sample, 1:] y_train = train.iloc[sample, 0] treereg.fit(X_train, y_train) y_pred[i] = treereg.predict(X_test) # In[16]: y_pred # Results of each tree # In[17]: for i in range(n_B): print(i, np.sqrt(mean_squared_error(y_pred[i], y_test))) # Results of the ensemble # In[18]: y_pred.mean(axis=1) # In[19]: np.sqrt(mean_squared_error(y_test, y_pred.mean(axis=1))) # ## Bagged decision trees in scikit-learn (with B=500) # In[20]: # define the training and testing sets X_train = train.iloc[:, 1:] y_train = train.iloc[:, 0] X_test = test.iloc[:, 1:] y_test = test.iloc[:, 0] # In[21]: # instruct BaggingRegressor to use DecisionTreeRegressor as the "base estimator" from sklearn.ensemble import BaggingRegressor bagreg = BaggingRegressor(DecisionTreeRegressor(), n_estimators=500, bootstrap=True, oob_score=True, random_state=1) # In[22]: # fit and predict bagreg.fit(X_train, y_train) y_pred = bagreg.predict(X_test) y_pred # In[23]: # calculate RMSE np.sqrt(mean_squared_error(y_test, y_pred)) # ## Estimating out-of-sample error # # For bagged models, out-of-sample error can be estimated without using **train/test split** or **cross-validation**! # # On average, each bagged tree uses about **two-thirds** of the observations. For each tree, the **remaining observations** are called "out-of-bag" observations. # In[24]: # show the first bootstrap sample samples[0] # In[25]: # show the "in-bag" observations for each sample for sample in samples: print(set(sample)) # In[26]: # show the "out-of-bag" observations for each sample for sample in samples: print(sorted(set(range(n_samples)) - set(sample))) # How to calculate **"out-of-bag error":** # # 1. For every observation in the training data, predict its response value using **only** the trees in which that observation was out-of-bag. Average those predictions (for regression) or take a vote (for classification). # 2. Compare all predictions to the actual response values in order to compute the out-of-bag error. # # When B is sufficiently large, the **out-of-bag error** is an accurate estimate of **out-of-sample error**. # In[27]: # compute the out-of-bag R-squared score (not MSE, unfortunately!) for B=500 bagreg.oob_score_ # ## Estimating feature importance # # Bagging increases **predictive accuracy**, but decreases **model interpretability** because it's no longer possible to visualize the tree to understand the importance of each feature. # # However, we can still obtain an overall summary of **feature importance** from bagged models: # # - **Bagged regression trees:** calculate the total amount that **MSE** is decreased due to splits over a given feature, averaged over all trees # - **Bagged classification trees:** calculate the total amount that **Gini index** is decreased due to splits over a given feature, averaged over all trees # # Part 4: Random Forests # # Random Forests is a **slight variation of bagged trees** that has even better performance: # # - Exactly like bagging, we create an ensemble of decision trees using bootstrapped samples of the training set. # - However, when building each tree, each time a split is considered, a **random sample of m features** is chosen as split candidates from the **full set of p features**. The split is only allowed to use **one of those m features**. # - A new random sample of features is chosen for **every single tree at every single split**. # - For **classification**, m is typically chosen to be the square root of p. # - For **regression**, m is typically chosen to be somewhere between p/3 and p. # # What's the point? # # - Suppose there is **one very strong feature** in the data set. When using bagged trees, most of the trees will use that feature as the top split, resulting in an ensemble of similar trees that are **highly correlated**. # - Averaging highly correlated quantities does not significantly reduce variance (which is the entire goal of bagging). # - By randomly leaving out candidate features from each split, **Random Forests "decorrelates" the trees**, such that the averaging process can reduce the variance of the resulting model. # # Part 5: Building and tuning decision trees and Random Forests # # - Major League Baseball player data from 1986-87: [data](https://github.com/justmarkham/DAT8/blob/master/data/hitters.csv), [data dictionary](https://cran.r-project.org/web/packages/ISLR/ISLR.pdf) (page 7) # - Each observation represents a player # - **Goal:** Predict player salary # In[28]: # read in the data with zipfile.ZipFile('../datasets/hitters.csv.zip', 'r') as z: f = z.open('hitters.csv') hitters = pd.read_csv(f, sep=',', index_col=False) # remove rows with missing values hitters.dropna(inplace=True) hitters.head() # In[29]: # encode categorical variables as integers hitters['League'] = pd.factorize(hitters.League)[0] hitters['Division'] = pd.factorize(hitters.Division)[0] hitters['NewLeague'] = pd.factorize(hitters.NewLeague)[0] hitters.head() # In[30]: # allow plots to appear in the notebook get_ipython().run_line_magic('matplotlib', 'inline') import matplotlib.pyplot as plt plt.style.use('fivethirtyeight') # In[31]: # scatter plot of Years versus Hits colored by Salary hitters.plot(kind='scatter', x='Years', y='Hits', c='Salary', colormap='jet', xlim=(0, 25), ylim=(0, 250)) # In[32]: # define features: exclude career statistics (which start with "C") and the response (Salary) feature_cols = hitters.columns[hitters.columns.str.startswith('C') == False].drop('Salary') feature_cols # In[33]: # define X and y X = hitters[feature_cols] y = hitters.Salary # ## Predicting salary with a decision tree # # Find the best **max_depth** for a decision tree using cross-validation: # In[34]: # list of values to try for max_depth max_depth_range = range(1, 21) # list to store the average RMSE for each value of max_depth RMSE_scores = [] # use 10-fold cross-validation with each value of max_depth from sklearn.model_selection import cross_val_score for depth in max_depth_range: treereg = DecisionTreeRegressor(max_depth=depth, random_state=1) MSE_scores = cross_val_score(treereg, X, y, cv=10, scoring='neg_mean_squared_error') RMSE_scores.append(np.mean(np.sqrt(-MSE_scores))) # In[35]: # plot max_depth (x-axis) versus RMSE (y-axis) plt.plot(max_depth_range, RMSE_scores) plt.xlabel('max_depth') plt.ylabel('RMSE (lower is better)') # In[36]: # show the best RMSE and the corresponding max_depth sorted(zip(RMSE_scores, max_depth_range))[0] # In[37]: # max_depth=2 was best, so fit a tree using that parameter treereg = DecisionTreeRegressor(max_depth=2, random_state=1) treereg.fit(X, y) # In[38]: # compute feature importances pd.DataFrame({'feature':feature_cols, 'importance':treereg.feature_importances_}).sort_values('importance') # ## Predicting salary with a Random Forest # In[39]: from sklearn.ensemble import RandomForestRegressor rfreg = RandomForestRegressor() rfreg # ### Tuning n_estimators # # One important tuning parameter is **n_estimators**, which is the number of trees that should be grown. It should be a large enough value that the error seems to have "stabilized". # In[40]: # list of values to try for n_estimators estimator_range = range(10, 310, 10) # list to store the average RMSE for each value of n_estimators RMSE_scores = [] # use 5-fold cross-validation with each value of n_estimators (WARNING: SLOW!) for estimator in estimator_range: rfreg = RandomForestRegressor(n_estimators=estimator, random_state=1, n_jobs=-1) MSE_scores = cross_val_score(rfreg, X, y, cv=5, scoring='neg_mean_squared_error') RMSE_scores.append(np.mean(np.sqrt(-MSE_scores))) # In[41]: # plot n_estimators (x-axis) versus RMSE (y-axis) plt.plot(estimator_range, RMSE_scores) plt.xlabel('n_estimators') plt.ylabel('RMSE (lower is better)') # ### Tuning max_features # # The other important tuning parameter is **max_features**, which is the number of features that should be considered at each split. # In[42]: # list of values to try for max_features feature_range = range(1, len(feature_cols)+1) # list to store the average RMSE for each value of max_features RMSE_scores = [] # use 10-fold cross-validation with each value of max_features (WARNING: SLOW!) for feature in feature_range: rfreg = RandomForestRegressor(n_estimators=150, max_features=feature, random_state=1, n_jobs=-1) MSE_scores = cross_val_score(rfreg, X, y, cv=10, scoring='neg_mean_squared_error') RMSE_scores.append(np.mean(np.sqrt(-MSE_scores))) # In[43]: # plot max_features (x-axis) versus RMSE (y-axis) plt.plot(feature_range, RMSE_scores) plt.xlabel('max_features') plt.ylabel('RMSE (lower is better)') # In[44]: # show the best RMSE and the corresponding max_features sorted(zip(RMSE_scores, feature_range))[0] # ### Fitting a Random Forest with the best parameters # In[ ]: # max_features=8 is best and n_estimators=150 is sufficiently large rfreg = RandomForestRegressor(n_estimators=150, max_features=8, max_depth=3, oob_score=True, random_state=1) rfreg.fit(X, y) # In[46]: # compute feature importances pd.DataFrame({'feature':feature_cols, 'importance':rfreg.feature_importances_}).sort_values('importance') # In[47]: # compute the out-of-bag R-squared score rfreg.oob_score_ # ### Reducing X to its most important features # # In[48]: # check the shape of X X.shape # In[49]: rfreg # In[50]: # set a threshold for which features to include from sklearn.feature_selection import SelectFromModel print(SelectFromModel(rfreg, threshold=0.1, prefit=True).transform(X).shape) print(SelectFromModel(rfreg, threshold='mean', prefit=True).transform(X).shape) print(SelectFromModel(rfreg, threshold='median', prefit=True).transform(X).shape) # In[51]: # create a new feature matrix that only includes important features X_important = SelectFromModel(rfreg, threshold='mean', prefit=True).transform(X) # In[53]: # check the RMSE for a Random Forest that only includes important features rfreg = RandomForestRegressor(n_estimators=150, max_features=3, random_state=1) scores = cross_val_score(rfreg, X_important, y, cv=10, scoring='neg_mean_squared_error') np.mean(np.sqrt(-scores)) # ## Comparing Random Forests with decision trees # # **Advantages of Random Forests:** # # - Performance is competitive with the best supervised learning methods # - Provides a more reliable estimate of feature importance # - Allows you to estimate out-of-sample error without using train/test split or cross-validation # # **Disadvantages of Random Forests:** # # - Less interpretable # - Slower to train # - Slower to predict