Ensemble Methods

All statistical learning and machine learning models face the challenge of balancing bias and variance.

The term "bias" refers to how much, on average, a model's predictions deviate from the true values. It reflects the model's ability to capture the true relationship between the predictors and the response variable. For example, if the relationship follows a nonlinear pattern, a linear regression model, regardless of how much data is available, will not be able to model the relationship properly and will have high bias.

On the other hand, the term "variance" refers to how much the model changes depending on the data used for training. Ideally, a model should not change much with small variations in the training data. If it does, it indicates that the model is memorizing the data instead of learning the true relationship between the predictors and the response variable. For example, a tree model with many nodes tends to change its structure even with small variations in the training data, indicating high variance.

As the complexity of a model increases, it becomes more flexible and can better fit the observations, which reduces bias and improves predictive ability. However, after a certain level of complexity, the problem of overfitting (overfitting) arises. This phenomenon occurs when the model fits the training data so closely that it cannot correctly predict new observations. The optimal model is the one that finds a suitable balance between bias and variance.

How are bias and variance controlled in tree-based models? Generally, small trees with few branches tend to have low variance but may not accurately capture the relationship between variables, resulting in high bias. On the other hand, large trees fit the training data very closely, reducing bias but increasing variance. One way to address this problem is with ensemble methods.

Ensemble methods combine multiple models into a new one with the goal of achieving a balance between bias and variance, thus obtaining better predictions than any of the individual original models. Two of the most widely used types of ensemble are:

• Bagging: Multiple models are fitted, each with a different subset of the training data. For prediction, all models in the ensemble contribute their prediction. The final value is the mean of all predictions (for continuous variables) or the most frequent class (for categorical variables). Random Forest models fall into this category.

• Boosting: Multiple simple models, called weak learners, are fitted sequentially so that each model learns from the errors of the previous one. As with bagging, the final value is the mean of all predictions (for continuous variables) or the most frequent class (for categorical variables). Three of the most commonly used boosting methods are AdaBoost, Gradient Boosting, and Stochastic Gradient Boosting.

Although the ultimate goal is the same—to achieve an optimal balance between bias and variance—there are two important differences:

• How they reduce total error. The total error of a model can be decomposed as $bias+variance+ϵ$. In bagging, models with very low bias but high variance are used; by aggregating them, variance is reduced without significantly increasing bias. In boosting, models with very low variance but high bias are used; by fitting models sequentially, bias is reduced. Therefore, each strategy reduces a different part of the total error.

• How variations are introduced in the models that make up the ensemble. In bagging, each model is different from the rest because each is trained with a different sample obtained by bootstrapping). In boosting, models are fitted sequentially and the importance (weight) of the observations changes in each iteration, resulting in different fits.

The key for ensemble methods to achieve better results than any of their individual models is that the models that make them up are as diverse as possible (their errors are not correlated). An analogy that illustrates this concept is the following: imagine a trivia game where teams have to answer questions on various topics. A team made up of many players, each an expert in a different topic, will have a better chance of winning than a team made up of players who are all experts in a single topic or a single player who knows a little about everything.

Next, the bagging strategy, on which the Random Forest model is based, is described in more detail.

Bagging

The term bagging is short for bootstrap aggregation, and refers to the use of repeated sampling with replacement bootstrapping) to reduce the variance of some statistical learning models, including tree-based models.

Given $n$ samples of independent observations $Z_1$, ..., $Z_n$, each with variance ${\sigma }^2$, the variance of the mean of the observations $\overline{Z}$ is ${\sigma }^2/n$. In other words, averaging a set of observations reduces variance. Based on this idea, one way to reduce variance and increase the accuracy of a predictive method is to obtain multiple samples from the population, fit a different model to each, and average the resulting predictions (or take the mode for qualitative variables). Since in practice we usually do not have access to multiple samples, the process can be simulated using bootstrapping), thus generating pseudo-samples to fit different models and then aggregate them. This process is known as bagging and is applicable to a wide variety of regression methods.

In the particular case of decision trees, given their low bias and high variance nature, bagging has shown very good results. The way to apply it is:

1. Generate $B$ pseudo-training sets using bootstrapping from the original training sample.

2. Train a tree with each of the $B$ samples from step 1. Each tree is created with almost no restrictions and is not pruned, so it has high variance but low bias. In most cases, the only stopping rule is the minimum number of observations that terminal nodes must have. The optimal value of this hyperparameter can be obtained by comparing the out of bag error or by cross-validation.

3. For each new observation, obtain the prediction from each of the $B$ trees. The final prediction is the mean of the $B$ predictions for quantitative variables and the most frequent predicted class (mode) for qualitative variables.

In the bagging process, the number of trees created is not a critical hyperparameter, since increasing the number does not increase the risk of overfitting. After a certain number of trees, the reduction in test error stabilizes. Nevertheless, each tree takes up memory, so it is not advisable to store more than necessary.

Training Random Forest

The Random Forest algorithm is a modification of the bagging process that improves results by further decorrelating the trees generated in the process.

As mentioned earlier, the benefits of bagging are based on the fact that averaging a set of models reduces variance. This is true as long as the aggregated models are not correlated. If the correlation is high, the reduction in variance that can be achieved is small.

Suppose a dataset has one very influential predictor, along with others that are moderately influential. In this scenario, all or almost all trees created in the bagging process will be dominated by the same predictor and will be very similar to each other. As a result of the high correlation between the trees, bagging will barely reduce variance and thus will not improve the model. Random forest avoids this problem by randomly selecting $m$ predictors before evaluating each split. In this way, on average, $(p-m)/p$ splits will not consider the influential predictor, allowing other predictors to be selected. By adding this extra step, the trees are further decorrelated, so their aggregation achieves a greater reduction in variance.

Random forest and bagging methods follow the same algorithm with the only difference being that, in random forest, before each split, $m$ predictors are randomly selected. The difference in the result will depend on the chosen value of $m$. If $m=p$, the results of random forest and bagging are equivalent. Some recommendations are:

• The square root of the total number of predictors for classification problems. $m\approx \sqrt{p}$
• One third of the number of predictors for regression problems. $m\approx \frac{p}{3}$
• If the predictors are highly correlated, small values of $m$ achieve better results.

However, the best way to find the optimal value of $m$ is to evaluate the out-of-bag-error or use cross-validation.

As with bagging, random forest does not suffer from overfitting by increasing the number of trees created in the process. After a certain number, the reduction in test error stabilizes.

Random Forest Prediction

The prediction of a Random Forest model is the mean of the predictions of all the trees that compose it.

Suppose there are 10 observations, each with a response variable $Y$ and some predictors $X$.

The following image shows how the model would predict for a new observation. In each tree, the path to the final node is highlighted. In each terminal node, the indices of the training observations that are part of it are detailed.

The value predicted by each tree is the mean of the response variable $Y$ in the terminal node. According to the image, the predictions of each of the three trees (from left to right) are:

The final prediction of the model is the mean of all individual predictions:

Although the above is the most common way to obtain predictions from a Random Forest model, there is another approach. The prediction of a regression tree can be seen as a variant of nearest neighbors in which only the observations that are part of the same terminal node as the predicted observation have influence. Following this approach, the tree prediction is defined as the weighted mean of all training observations, where the weight of each observation depends only on whether it is part of the same terminal node. For Random Forest, this is equivalent to the weighted mean of all observations, using as weights $\text{w}$ the mean of the weight vectors of all trees.

According to the previous image, the weight vector for each of the three trees (from left to right) is:

The mean of all weight vectors is:

Once the average weight vector is obtained, the prediction can be calculated as the weighted mean of all training observations:

Out-of-Bag Error

Due to the nature of the bagging process, it is possible to estimate the test error without the need to use cross-validation methods. The fact that trees are fitted using samples generated by bootstrapping means that, on average, each fit uses only about two-thirds of the original observations. The remaining third is called out-of-bag (OOB).

If for each tree fitted in the bagging process the observations used are recorded, the response for observation $i$ can be predicted using only those trees in which that observation was excluded and averaging them (the mode in the case of classification trees). Following this process, predictions can be obtained for all $n$ observations and the OOB-mean square error (for regression) or the OOB-classification error (for classification trees) can be calculated. Since the response variable for each observation is predicted using only the trees in which it did not participate, the OOB-error serves as an estimate of the test error. In fact, if the number of trees is high enough, the OOB-error is practically equivalent to the leave-one-out cross-validation error.

This is an added advantage of bagging methods, and therefore of Random Forest, as it avoids the need to use the (computationally expensive) cross-validation process for hyperparameter optimization.

Two limitations in the use of Out-of-Bag Error:

• The Out-of-Bag Error is not suitable when observations have a temporal relationship (time series). Since the selection of observations for each training is random, the temporal order is not respected and information from the future would be introduced.

• The preprocessing of the training data is done jointly, so the out-of-bag observations may suffer from data leakage). If so, the OOB-error estimates are too optimistic.

In a bootstrapping sample, if the training data size is $n$, each observation has a probability of being chosen of $\frac{1}{n}$. Therefore, the probability of not being chosen throughout the process is $(1-1/n{)}^n$, which converges to $1/ϵ$, which is approximately one third.

Predictor Importance

While it is true that the bagging (Random Forest) process improves predictive ability compared to single-tree models, this comes at a cost: model interpretability is reduced. Since it is a combination of multiple trees, it is not possible to obtain a simple graphical representation of the model, and it is not straightforward to visually identify the most important predictors. However, new strategies have been developed to quantify predictor importance, making bagging (Random Forest) models a powerful tool not only for prediction but also for exploratory analysis. Two such measures are: permutation importance and node impurity.

Permutation Importance

Identifies the influence of each predictor on a given model evaluation metric (estimated by out-of-bag error or cross-validation). The value associated with each predictor is obtained as follows:

1. Build the set of trees that make up the model.

2. Calculate a given error metric (mse, classification error, ...). This is the reference value ($erro{r}_0$).

3. For each predictor $j$:

   • Permute the values of predictor $j$ in all trees of the model, keeping the rest constant.
   • Recalculate the metric after permutation, call this ($erro{r}_j$).
   • Calculate the increase in the metric due to permutation of predictor $j$. $${\text{Increase}}_j=\frac{(erro{r}_j-erro{r}_0)}{erro{r}_0}\ast 100\%$$

If the permuted predictor was contributing to the model, it is expected that the model error will increase, as the information provided by that variable is lost. The percentage increase in error due to permutation of predictor $j$ can be interpreted as the influence of $j$ on the model. Something that often causes confusion is that this increase can be negative. If the variable does not contribute to the model, it is possible that, by shuffling it randomly, the model may improve slightly just by chance, so $(erro{r}_j-erro{r}_0)$ is negative. In general, these variables can be considered to have importance close to zero.

Although this strategy is usually the most recommended, some caution is needed in its interpretation. What it quantifies is the influence of the predictors on the model, not their relationship with the response variable. Why is this important? Suppose a scenario where this strategy is used to identify which predictors are related to a person's weight, and two of the predictors are: body mass index (BMI) and height. Since BMI and height are highly correlated (the information they provide is redundant), when one of them is permuted, the impact on the model will be minimal, as the other provides the same information. As a result, these predictors will appear as not very influential even though they are actually highly related to the response variable. One way to avoid such problems is, whenever predictors are excluded from a model, to check the impact on its predictive ability.

Increase in Node Purity

Quantifies the total increase in node purity due to splits involving the predictor (averaged over all trees). The way to calculate it is as follows: at each split in the trees, the decrease achieved in the measure used as the split criterion (Gini index, mse, entropy, ...) is recorded. For each predictor, the average decrease achieved in the set of trees that make up the ensemble is calculated. The higher this average value, the greater the contribution of the predictor to the model.

Libraries

The libraries used in this example are:

# Data processing
# ==============================================================================
import numpy as np
import pandas as pd

# Plots
# ==============================================================================
import matplotlib.pyplot as plt
plt.style.use('fivethirtyeight')
plt.rcParams['lines.linewidth'] = 1.5
plt.rcParams['font.size'] = 8

# Preprocessing and modeling
# ==============================================================================
import sklearn
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import root_mean_squared_error
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import train_test_split
from sklearn.model_selection import RepeatedKFold
from sklearn.model_selection import GridSearchCV
from sklearn.model_selection import ParameterGrid
from sklearn.model_selection import cross_val_score
from sklearn.inspection import permutation_importance
from joblib import Parallel, delayed, cpu_count
import optuna

# Warnings configuration
# ==============================================================================
import warnings
optuna.logging.set_verbosity(optuna.logging.WARNING)

print(f"scikit-learn version: {sklearn.__version__}")

scikit-learn version: 1.7.1

Data

The Boston dataset contains housing prices in the city of Boston, as well as socio-economic information about the neighborhood in which they are located. The goal is to fit a regression model to predict the median price of a house (MEDV) based on the available variables.

Number of instances: 506

Number of attributes: 13 predictive numerical/categorical.

Attribute information (in order):

• CRIM: per capita crime rate by town.
• ZN: proportion of residential land zoned for lots over 25,000 sq.ft.
• INDUS: proportion of non-retail business acres per town.
• CHAS: Charles River dummy variable (= 1 if tract bounds river; 0 otherwise).
• NOX: nitric oxides concentration (parts per 10 million).
• RM: average number of rooms per dwelling.
• AGE: proportion of owner-occupied units built prior to 1940.
• DIS: weighted distances to five Boston employment centers.
• RAD: index of accessibility to radial highways.
• TAX: full-value property-tax rate per $10,000.
• PTRATIO: pupil-teacher ratio by town.
• LSTAT: % lower status of the population.
• MEDV: Median value of owner-occupied homes in $1000's.

Missing values: None

# Data download
# ==============================================================================
url = (
    'https://raw.githubusercontent.com/JoaquinAmatRodrigo/Estadistica-machine-learning-python/'
    'master/data/Boston.csv'
    )
data = pd.read_csv(url, sep=',')
display(data.head(3))
display(data.info())

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 506 entries, 0 to 505
Data columns (total 13 columns):
 #   Column   Non-Null Count  Dtype  
---  ------   --------------  -----  
 0   CRIM     506 non-null    float64
 1   ZN       506 non-null    float64
 2   INDUS    506 non-null    float64
 3   CHAS     506 non-null    int64  
 4   NOX      506 non-null    float64
 5   RM       506 non-null    float64
 6   AGE      506 non-null    float64
 7   DIS      506 non-null    float64
 8   RAD      506 non-null    int64  
 9   TAX      506 non-null    int64  
 10  PTRATIO  506 non-null    float64
 11  LSTAT    506 non-null    float64
 12  MEDV     506 non-null    float64
dtypes: float64(10), int64(3)
memory usage: 51.5 KB

None

Model Fitting

A model is fitted using MEDV as the response variable and all other available variables as predictors.

The RandomForestRegressor class from the sklearn.ensemble module allows you to train random forest models for regression problems. The default parameters and hyperparameters used are:

• n_estimators=100
• criterion='squared_error'
• max_depth=None
• min_samples_split=2
• min_samples_leaf=1
• min_weight_fraction_leaf=0.0
• max_features=1.0
• max_leaf_nodes=None
• min_impurity_decrease=0.0
• bootstrap=True
• oob_score=False
• n_jobs=None
• random_state=None
• verbose=0
• warm_start=False
• ccp_alpha=0.0
• max_samples=None

Among all these, those that stop tree growth, control the number of trees and predictors included, and manage parallelization are particularly important:

• n_estimators: number of trees included in the model.
• max_depth: maximum depth the trees can reach.
• min_samples_split: minimum number of observations a node must have to be split. If a decimal value, it is interpreted as a fraction of the total training observations ceil(min_samples_split * n_samples).
• min_samples_leaf: minimum number of observations each child node must have for the split to occur. If a decimal value, it is interpreted as a fraction of the total training observations ceil(min_samples_split * n_samples).
• max_leaf_nodes: maximum number of terminal nodes the trees can have.

• max_features: number of predictors considered at each split. Can be:

  • An integer value
  • A fraction of the total predictors.
  • “sqrt”, square root of the total number of predictors.
  • “log2”, log2 of the total number of predictors.
  • None, uses all predictors.

• oob_score: Whether to calculate the out-of-bag R^2. By default, it is False as it increases training time.

• n_jobs: number of cores used for training. In random forest, trees are fitted independently, so parallelization significantly reduces training time. With -1, all available cores are used.
• random_state: seed for reproducibility. Must be an integer value.

As in any predictive study, it is important not only to fit the model but also to quantify its ability to predict new observations. To make this evaluation, the data is split into two groups: one for training and one for testing.

# Split data into training and test sets
# ==============================================================================
X_train, X_test, y_train, y_test = train_test_split(
    data.drop(columns="MEDV"),
    data['MEDV'],
    test_size = 0.25,
    random_state = 123
)

print(f"Training set size: {X_train.shape[0]}")
print(f"Test set size: {X_test.shape[0]}")

# Model creation
# ==============================================================================
model = RandomForestRegressor(
    n_estimators = 10,
    criterion    = 'squared_error',
    max_depth    = None,
    max_features = 1,
    oob_score    = False,
    n_jobs       = -1,
    random_state = 123
)

# Model training
# ==============================================================================
model.fit(X_train, y_train)

Training set size: 379
Test set size: 127

RandomForestRegressor(max_features=1, n_estimators=10, n_jobs=-1,
                      random_state=123)

Model Prediction and Evaluation

Once the model is trained, its predictive performance is evaluated using the test set.

# Test error of the initial model
# ==============================================================================
predictions = model.predict(X=X_test)
rmse = root_mean_squared_error(y_true=y_test, y_pred=predictions)
print(f"The test error (rmse) is: {rmse}")

The test error (rmse) is: 4.17800427058654

Hyperparameter Optimization

The initial model was trained using 10 trees (n_estimators=10) and keeping the rest of the hyperparameters at their default values. Since these are hyperparameters, the most suitable value cannot be known in advance; the way to identify them is by using validation strategies, such as cross-validation.

Random Forest models have the advantage of providing the Out-of-Bag error, which allows for an estimate of the test error without resorting to cross-validation, which is computationally expensive. In the RandomForestRegressor implementation, the metric returned as oob_score is $R^2$; if another metric is desired, the oob_decision_function_() method must be used to obtain the predictions and then calculate the metric of interest. For a more detailed explanation, see: Grid search for Random Forest models with out-of-bag error and early stopping.

Keep in mind that when searching for the optimal value of a hyperparameter with two different metrics, the result is rarely the same. The important thing is that both metrics identify the same regions of interest.

Number of Trees

In Random Forest, the number of trees is not a critical hyperparameter, since adding trees can only improve the result. Random Forest does not suffer from overfitting due to an excess of trees. However, adding trees once the improvement has stabilized is a waste of computational resources.

Hide

# Validation using Out-of-Bag error
# ==============================================================================
warnings.filterwarnings('ignore')
train_scores = []
oob_scores   = []

# Evaluated values
estimator_range = range(1, 150, 5)

# Loop to train a model with each value of n_estimators and extract its training and Out-of-Bag error.
for n_estimators in estimator_range:
    model = RandomForestRegressor(
        n_estimators = n_estimators,
        criterion    = 'squared_error',
        max_depth    = None,
        max_features = 1,
        oob_score    = True,
        n_jobs       = -1,
        random_state = 123
    )
    model.fit(X_train, y_train)
    train_scores.append(model.score(X_train, y_train))
    oob_scores.append(model.oob_score_)
    
# Plot with the evolution of errors
fig, ax = plt.subplots(figsize=(5, 3))
ax.plot(estimator_range, train_scores, label="train scores")
ax.plot(estimator_range, oob_scores, label="out-of-bag scores")
ax.plot(estimator_range[np.argmax(oob_scores)], max(oob_scores),
        marker='o', color = "red", label="max score")
ax.set_ylabel("R^2")
ax.set_xlabel("n_estimators")
ax.set_title("Out-of-bag-error evolution vs number of trees")
plt.legend();
print(f"Optimal value of n_estimators: {estimator_range[np.argmax(oob_scores)]}")
warnings.filterwarnings('default')

Optimal value of n_estimators: 66

Hide

# Validation using k-cross-validation and neg_root_mean_squared_error
# ==============================================================================
train_scores = []
cv_scores    = []

# Evaluated values
estimator_range = range(1, 150, 5)

# Loop to train a model with each value of n_estimators and extract its training and k-cross-validation error.
for n_estimators in estimator_range:
    
    model = RandomForestRegressor(
        n_estimators = n_estimators,
        criterion    = 'squared_error',
        max_depth    = None,
        max_features = 1,
        oob_score    = False,
        n_jobs       = -1,
        random_state = 123
    )
    
    # Training error
    model.fit(X_train, y_train)
    predictions = model.predict(X=X_train)
    rmse = root_mean_squared_error(
        y_true  = y_train,
        y_pred  = predictions,
    )
    train_scores.append(rmse)
    
    # Cross-validation error
    scores = cross_val_score(
        estimator = model,
        X         = X_train,
        y         = y_train,
        scoring   = 'neg_root_mean_squared_error',
        cv        = 5
    )
    # Add cross_val_score() scores and convert to positive
    cv_scores.append(-1*scores.mean())
    
# Plot with the evolution of errors
fig, ax = plt.subplots(figsize=(5, 3))
ax.plot(estimator_range, train_scores, label="train scores")
ax.plot(estimator_range, cv_scores, label="cv scores")
ax.plot(estimator_range[np.argmin(cv_scores)], min(cv_scores),
        marker='o', color = "red", label="min score")
ax.set_ylabel("root_mean_squared_error")
ax.set_xlabel("n_estimators")
ax.set_title("cv-error evolution vs number of trees")
plt.legend();
print(f"Optimal value of n_estimators: {estimator_range[np.argmin(cv_scores)]}")

Optimal value of n_estimators: 96

Although the optimal value of the metrics is reached with 66 and 96 trees, the curves indicate that, from 20 trees onward, the model's validation error stabilizes.

Max features

The value of max_features is one of the most important hyperparameters in random forest, as it controls how much the trees are decorrelated from each other.

Hide

# Validation using Out-of-Bag error
# ==============================================================================
train_scores = []
oob_scores   = []

# Evaluated values
max_features_range = range(1, X_train.shape[1] + 1, 1)

# Loop to train a model with each value of max_features and extract its training and Out-of-Bag error.
for max_features in max_features_range:
    model = RandomForestRegressor(
        n_estimators = 100,
        criterion    = 'squared_error',
        max_depth    = None,
        max_features = max_features,
        oob_score    = True,
        n_jobs       = -1,
        random_state = 123
    )
    model.fit(X_train, y_train)
    train_scores.append(model.score(X_train, y_train))
    oob_scores.append(model.oob_score_)
    
# Plot with the evolution of errors
fig, ax = plt.subplots(figsize=(5, 3))
ax.plot(max_features_range, train_scores, label="train scores")
ax.plot(max_features_range, oob_scores, label="out-of-bag scores")
ax.plot(max_features_range[np.argmax(oob_scores)], max(oob_scores),
        marker='o', color = "red")
ax.set_ylabel("R^2")
ax.set_xlabel("max_features")
ax.set_title("Out-of-bag-error evolution vs number of predictors")
plt.legend();
print(f"Optimal value of max_features: {max_features_range[np.argmax(oob_scores)]}")

Optimal value of max_features: 4

Hide

# Validation using k-cross-validation and neg_root_mean_squared_error
# ==============================================================================
train_scores = []
cv_scores    = []

# Evaluated values
max_features_range = range(1, X_train.shape[1] + 1, 1)

# Loop to train a model with each value of max_features and extract its training and k-cross-validation error.
for max_features in max_features_range:
    
    model = RandomForestRegressor(
        n_estimators = 100,
        criterion    = 'squared_error',
        max_depth    = None,
        max_features = max_features,
        oob_score    = True,
        n_jobs       = -1,
        random_state = 123
    )
    
    # Training error
    model.fit(X_train, y_train)
    predictions = model.predict(X = X_train)
    rmse = root_mean_squared_error(
        y_true  = y_train,
        y_pred  = predictions,
    )
    train_scores.append(rmse)
    
    # Cross-validation error
    scores = cross_val_score(
        estimator = model,
        X         = X_train,
        y         = y_train,
        scoring   = 'neg_root_mean_squared_error',
        cv        = 5
    )
    # Add cross_val_score() scores and convert to positive
    cv_scores.append(-1*scores.mean())
    
# Plot with the evolution of errors
fig, ax = plt.subplots(figsize=(5, 3))
ax.plot(max_features_range, train_scores, label="train scores")
ax.plot(max_features_range, cv_scores, label="cv scores")
ax.plot(max_features_range[np.argmin(cv_scores)], min(cv_scores),
        marker='o', color = "red", label="min score")
ax.set_ylabel("root_mean_squared_error")
ax.set_xlabel("max_features")
ax.set_title("cv-error evolution vs number of predictors")
plt.legend();
print(f"Optimal value of max_features: {max_features_range[np.argmin(cv_scores)]}")

Optimal value of max_features: 6

According to the two metrics used, the optimal value of max_features is between 4 and 6.

Grid search

Although individual analysis of hyperparameters is useful for understanding their impact on the model and identifying ranges of interest, the final search should not be done sequentially, as each hyperparameter interacts with the others. It is preferable to use grid search, random search, or bayesian search to analyze several combinations of hyperparameters. More information about search strategies can be found in Machine learning with Python and Scikit-learn.

# Grid of evaluated hyperparameters
# ==============================================================================
param_grid = ParameterGrid(
    {'n_estimators': [150],
     'max_features': [5, 7, 9],
     'max_depth'   : [None, 3, 10, 20]
    }
)

# Loop to fit a model with each combination of hyperparameters
# ==============================================================================
results = {'params': [], 'oob_r2': []}

for params in param_grid:
    model = RandomForestRegressor(
        oob_score    = True,
        n_jobs       = -1,
        random_state = 123,
        ** params
    )
    model.fit(X_train, y_train)
    results['params'].append(params)
    results['oob_r2'].append(model.oob_score_)
    print(f"Model: {params} ✓")

# Results
# ==============================================================================
results = pd.DataFrame(results)
results = pd.concat([results, results['params'].apply(pd.Series)], axis=1)
results = results.drop(columns = 'params')
results = results.sort_values('oob_r2', ascending=False)
results.head(4)

Model: {'max_depth': None, 'max_features': 5, 'n_estimators': 150} ✓
Model: {'max_depth': None, 'max_features': 7, 'n_estimators': 150} ✓
Model: {'max_depth': None, 'max_features': 9, 'n_estimators': 150} ✓
Model: {'max_depth': 3, 'max_features': 5, 'n_estimators': 150} ✓
Model: {'max_depth': 3, 'max_features': 7, 'n_estimators': 150} ✓
Model: {'max_depth': 3, 'max_features': 9, 'n_estimators': 150} ✓
Model: {'max_depth': 10, 'max_features': 5, 'n_estimators': 150} ✓
Model: {'max_depth': 10, 'max_features': 7, 'n_estimators': 150} ✓
Model: {'max_depth': 10, 'max_features': 9, 'n_estimators': 150} ✓
Model: {'max_depth': 20, 'max_features': 5, 'n_estimators': 150} ✓
Model: {'max_depth': 20, 'max_features': 7, 'n_estimators': 150} ✓
Model: {'max_depth': 20, 'max_features': 9, 'n_estimators': 150} ✓

# Grid of evaluated hyperparameters
# ==============================================================================
param_grid = ParameterGrid(
    {'n_estimators': [150],
     'max_features': [5, 7, 9],
     'max_depth'   : [None, 3, 10, 20]
    }
)

# Parallel loop to fit a model with each combination of hyperparameters
# ==============================================================================
def eval_oob_error(X, y, params, verbose=True):
    """
    Function to train a model using given parameters
    and return the out-of-bag error
    """
    model = RandomForestRegressor(
        oob_score    = True,
        n_jobs       = -1,
        random_state = 123,
        ** params            
    )
    model.fit(X, y)
    if verbose:
        print(f"Model: {params} ✓")
    return{'params': params, 'oob_r2': model.oob_score_}

results = Parallel(n_jobs=cpu_count()-1)(
    delayed(eval_oob_error)(X_train, y_train, params)
    for params in param_grid
)

# Results
# ==============================================================================
results = pd.DataFrame(results)
results = pd.concat([results, results['params'].apply(pd.Series)], axis=1)
results = results.drop(columns = 'params')
results = results.sort_values('oob_r2', ascending=False)
results.head(4)

Model: {'max_depth': None, 'max_features': 5, 'n_estimators': 150} ✓
Model: {'max_depth': 3, 'max_features': 5, 'n_estimators': 150} ✓
Model: {'max_depth': 3, 'max_features': 7, 'n_estimators': 150} ✓
Model: {'max_depth': 10, 'max_features': 5, 'n_estimators': 150} ✓
Model: {'max_depth': None, 'max_features': 9, 'n_estimators': 150} ✓
Model: {'max_depth': None, 'max_features': 7, 'n_estimators': 150} ✓
Model: {'max_depth': 3, 'max_features': 9, 'n_estimators': 150} ✓
Model: {'max_depth': 10, 'max_features': 9, 'n_estimators': 150} ✓
Model: {'max_depth': 10, 'max_features': 7, 'n_estimators': 150} ✓
Model: {'max_depth': 20, 'max_features': 7, 'n_estimators': 150} ✓
Model: {'max_depth': 20, 'max_features': 5, 'n_estimators': 150} ✓
Model: {'max_depth': 20, 'max_features': 9, 'n_estimators': 150} ✓

# Best hyperparameters by out-of-bag error
# ==============================================================================
print("--------------------------------------------")
print("Best hyperparameters found (oob-r2)")
print("--------------------------------------------")
print(results.iloc[0,0:])

--------------------------------------------
Best hyperparameters found (oob-r2)
--------------------------------------------
oob_r2            0.882461
max_depth        10.000000
max_features      5.000000
n_estimators    150.000000
Name: 6, dtype: float64

# Grid of evaluated hyperparameters
# ==============================================================================
param_grid = {'n_estimators': [150],
              'max_features': [5, 7, 9],
              'max_depth'   : [None, 3, 10, 20]
             }

# Grid search with cross-validation
# ==============================================================================
grid = GridSearchCV(
    estimator  = RandomForestRegressor(random_state = 123),
    param_grid = param_grid,
    scoring    = 'neg_root_mean_squared_error',
    n_jobs     = cpu_count() - 1,
    cv         = RepeatedKFold(n_splits=5, n_repeats=3, random_state=123), 
    refit      = True,
    verbose    = 0,
    return_train_score = True
)

grid.fit(X=X_train, y=y_train)

# Results
# ==============================================================================
results = pd.DataFrame(grid.cv_results_)
results.filter(regex = '(param.*|mean_t|std_t)') \
    .drop(columns = 'params') \
    .sort_values('mean_test_score', ascending = False) \
    .head(4)

# Best hyperparameters found by cross-validation
# ==============================================================================
print("----------------------------------------")
print("Best hyperparameters found (cv)")
print("----------------------------------------")
print(grid.best_params_, ":", grid.best_score_, grid.scoring)

----------------------------------------
Best hyperparameters found (cv)
----------------------------------------
{'max_depth': 10, 'max_features': 5, 'n_estimators': 150} : -3.320294830700744 neg_root_mean_squared_error

# Test error of the final model
# ==============================================================================
final_model = grid.best_estimator_
predictions = final_model.predict(X=X_test)
rmse = root_mean_squared_error(
    y_true  = y_test,
    y_pred  = predictions,
)
print(f"The test error (rmse) is: {rmse}")

The test error (rmse) is: 3.550875871674505

Bayesian search

Grid search and random search can yield good results, especially when the search range is reduced. However, neither takes into account the results obtained so far, which prevents them from focusing the search on the most interesting regions and avoiding unnecessary ones.

An alternative is to use Bayesian optimization methods to search for hyperparameters. In general terms, Bayesian hyperparameter optimization consists of creating a probabilistic model in which the objective function is the model's validation metric (RMSE, AUC, accuracy, etc.). With this strategy, the search is redirected in each iteration to the most promising regions. The ultimate goal is to reduce the number of hyperparameter combinations evaluated by the model, choosing only the best candidates. This approach is especially advantageous when the search space is very large or model evaluation is very slow.

# Bayesian hyperparameter search with optuna
# ==============================================================================
def objective(trial):
    params = {
        'n_estimators': trial.suggest_int('n_estimators', 100, 1000, step=100),
        'max_depth': trial.suggest_int('max_depth', 3, 30),
        'min_samples_split': trial.suggest_int('min_samples_split', 2, 100),
        'min_samples_leaf': trial.suggest_int('min_samples_leaf', 1, 100),
        'max_features': trial.suggest_float('max_features', 0.2, 1.0),
        'ccp_alpha': trial.suggest_float('ccp_alpha', 0.0, 1),
        # Fixed parameters
        'n_jobs': -1,
        'random_state': 4576688,
    }

    model = RandomForestRegressor(**params)
    cross_val_scores = cross_val_score(
        estimator = model,
        X = X_train,
        y = y_train,
        cv = RepeatedKFold(n_splits=5, n_repeats=3, random_state=123),
        scoring = 'neg_root_mean_squared_error',
        n_jobs=-1
    )
    score = np.mean(cross_val_scores) 
    return score

study = optuna.create_study(direction='maximize') # Maximizing because the score is negative
study.optimize(objective, n_trials=30, show_progress_bar=True, timeout=60*5)

print('Best hyperparameters:', study.best_params)
print('Best score:', study.best_value)

  0%|          | 0/30 [00:00<?, ?it/s]

Best hyperparameters: {'n_estimators': 300, 'max_depth': 20, 'min_samples_split': 2, 'min_samples_leaf': 2, 'max_features': 0.5244309455082777, 'ccp_alpha': 0.16943620121169578}
Best score: -3.492724388695251

# Test error of the final model
# ==============================================================================
final_model = RandomForestRegressor(**study.best_params)
final_model.fit(X_train, y_train)
predictions = final_model.predict(X=X_test)
rmse = root_mean_squared_error(y_true=y_test, y_pred=predictions)
print(f"The test error (rmse) is: {rmse}")

The test error (rmse) is: 3.8358066178098786

Predictor importance

# Predictor importance based on impurity reduction
# ==============================================================================
feature_importance = pd.DataFrame({
    'predictor': data.drop(columns = "MEDV").columns,
    'importance': final_model.feature_importances_
})
print("Predictor importance in the model")
print("-------------------------------------------")
feature_importance.sort_values('importance', ascending=False)

Predictor importance in the model
-------------------------------------------

Permutation

Hide

# Predictor importance based on permutation
# ==============================================================================
importance = permutation_importance(
    estimator    = final_model,
    X            = X_train,
    y            = y_train,
    n_repeats    = 5,
    scoring      = 'neg_root_mean_squared_error',
    n_jobs       = cpu_count() - 1,
    random_state = 123
)

# Store the results (mean and std) in a dataframe
df_importance = pd.DataFrame(
    {k: importance[k] for k in ['importances_mean', 'importances_std']}
    )
df_importance['feature'] = X_train.columns
df_importance.sort_values('importances_mean', ascending=False)

	importances_mean	importances_std	feature
11	4.838738	0.223503	LSTAT
5	4.209149	0.073740	RM
4	0.716986	0.051328	NOX
7	0.624620	0.043498	DIS
0	0.462018	0.037184	CRIM
10	0.372917	0.033469	PTRATIO
9	0.234956	0.017121	TAX
2	0.223568	0.028646	INDUS
6	0.179764	0.021913	AGE
8	0.048013	0.003933	RAD
1	0.008110	0.001410	ZN
3	0.003015	0.000434	CHAS

Hide

# Plot of predictor importance
# ==============================================================================
fig, ax = plt.subplots(figsize=(3.5, 4))
df_importance = df_importance.sort_values('importances_mean', ascending=True)
ax.barh(
    df_importance['feature'],
    df_importance['importances_mean'],
    xerr=df_importance['importances_std'],
    align='center',
    alpha=0
)
ax.plot(
    df_importance['importances_mean'],
    df_importance['feature'],
    marker="D",
    linestyle="",
    alpha=0.8,
    color="r"
)
ax.set_title('Predictor importance (train)')
ax.set_xlabel('Increase in error after permutation');

Both strategies identify LSTAT and RM as the most influential predictors, according to the training data.

Libraries

The libraries used in this document are:

# Data processing
# ==============================================================================
import numpy as np
import pandas as pd
import statsmodels.api as sm

# Plots
# ==============================================================================
import matplotlib.pyplot as plt
plt.style.use('fivethirtyeight')
plt.rcParams['lines.linewidth'] = 1.5
plt.rcParams['font.size'] = 8

# Preprocessing and modeling
# ==============================================================================
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import accuracy_score
from sklearn.metrics import confusion_matrix
from sklearn.metrics import ConfusionMatrixDisplay
from sklearn.metrics import classification_report
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import OneHotEncoder
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import train_test_split
from sklearn.model_selection import RepeatedKFold
from sklearn.model_selection import GridSearchCV
from sklearn.model_selection import ParameterGrid
from sklearn.inspection import permutation_importance
from joblib import Parallel, delayed
import optuna

# Warnings configuration
# ==============================================================================
import warnings

Data

The Carseats dataset, originally from the R package ISLR and accessible in Python via statsmodels.datasets.get_rdataset, contains information on child seat sales in 400 different stores. For each of the 400 stores, 11 variables have been recorded. The goal is to build a classification model to predict whether a store has high sales (Sales > 8) or low sales (Sales <= 8) based on all available variables.

# Data loading
# ==============================================================================
carseats_data = sm.datasets.get_rdataset("Carseats", "ISLR")
data = carseats_data.data
print(carseats_data.__doc__)

.. container::

   .. container::

      ======== ===============
      Carseats R Documentation
      ======== ===============

      .. rubric:: Sales of Child Car Seats
         :name: sales-of-child-car-seats

      .. rubric:: Description
         :name: description

      A simulated data set containing sales of child car seats at 400
      different stores.

      .. rubric:: Usage
         :name: usage

      .. code:: R

         Carseats

      .. rubric:: Format
         :name: format

      A data frame with 400 observations on the following 11 variables.

      ``Sales``
         Unit sales (in thousands) at each location

      ``CompPrice``
         Price charged by competitor at each location

      ``Income``
         Community income level (in thousands of dollars)

      ``Advertising``
         Local advertising budget for company at each location (in
         thousands of dollars)

      ``Population``
         Population size in region (in thousands)

      ``Price``
         Price company charges for car seats at each site

      ``ShelveLoc``
         A factor with levels ``Bad``, ``Good`` and ``Medium``
         indicating the quality of the shelving location for the car
         seats at each site

      ``Age``
         Average age of the local population

      ``Education``
         Education level at each location

      ``Urban``
         A factor with levels ``No`` and ``Yes`` to indicate whether the
         store is in an urban or rural location

      ``US``
         A factor with levels ``No`` and ``Yes`` to indicate whether the
         store is in the US or not

      .. rubric:: Source
         :name: source

      Simulated data

      .. rubric:: References
         :name: references

      James, G., Witten, D., Hastie, T., and Tibshirani, R. (2013) *An
      Introduction to Statistical Learning with applications in R*,
      https://www.statlearning.com, Springer-Verlag, New York

      .. rubric:: Examples
         :name: examples

      .. code:: R

         summary(Carseats)
         lm.fit=lm(Sales~Advertising+Price,data=Carseats)

data.head(3)

Since Sales is a continuous variable and the goal of the study is to classify stores according to whether they sell a lot or a little, a new binary variable (0, 1) called high_sales is created.

data['high_sales'] = np.where(data.Sales > 8, 1, 0)
# Once the new response variable is created, drop the original
data = data.drop(columns = 'Sales')

Model fitting and hyperparameter optimization

A classification tree is fitted using high_sales as the response variable and all available variables as predictors. First, the hyperparameters max_depth=5 and criterion='gini' are used, with the rest left at their default values. Then, pruning is applied and the results are compared to the initial model.

Unlike the previous example, these data contain categorical variables, so before training the model, it is necessary to apply one-hot encoding. A more detailed description of this process can be found in Machine learning with Python and Scikit-learn.

# Split data into train and test sets
# ==============================================================================
X_train, X_test, y_train, y_test = train_test_split(
    data.drop(columns = 'high_sales'),
    data['high_sales'],
    random_state = 123
)

# One-hot encoding of categorical variables
# ==============================================================================
# Identify the names of numeric and categorical columns
cat_cols = X_train.select_dtypes(include=['object', 'category']).columns.to_list()
numeric_cols = X_train.select_dtypes(include=['float64', 'int']).columns.to_list()

# Apply one-hot encoding only to categorical columns
preprocessor = ColumnTransformer(
    [('onehot', OneHotEncoder(handle_unknown='ignore', sparse_output=False, drop='if_binary'), cat_cols)],
    remainder='passthrough',
    verbose_feature_names_out=False
).set_output(transform="pandas")

# Once the ColumnTransformer object is defined, use fit() to learn the transformations with the training data and apply them to both sets with transform(). Both operations at once with fit_transform().
X_train_prep = preprocessor.fit_transform(X_train)
X_test_prep  = preprocessor.transform(X_test)

X_train_prep.info()

<class 'pandas.core.frame.DataFrame'>
Index: 300 entries, 170 to 365
Data columns (total 12 columns):
 #   Column            Non-Null Count  Dtype  
---  ------            --------------  -----  
 0   ShelveLoc_Bad     300 non-null    float64
 1   ShelveLoc_Good    300 non-null    float64
 2   ShelveLoc_Medium  300 non-null    float64
 3   Urban_Yes         300 non-null    float64
 4   US_Yes            300 non-null    float64
 5   CompPrice         300 non-null    int64  
 6   Income            300 non-null    int64  
 7   Advertising       300 non-null    int64  
 8   Population        300 non-null    int64  
 9   Price             300 non-null    int64  
 10  Age               300 non-null    int64  
 11  Education         300 non-null    int64  
dtypes: float64(5), int64(7)
memory usage: 30.5 KB

Although RandomForestClassifier has default values for its hyperparameters, it cannot be known in advance if these are the most suitable; the way to identify them is by using validation strategies, such as cross-validation.

Random Forest models have the advantage of providing the Out-of-Bag error, which allows for an estimate of the test error without resorting to cross-validation, which is computationally expensive. In the RandomForestClassifier implementation, the metric returned as oob_score is accuracy; if another metric is desired, the oob_decision_function_() method must be used to obtain the predictions and then calculate the metric of interest. For a more detailed explanation, see: Grid search for Random Forest models with out-of-bag error and early stopping.

Keep in mind that when searching for the optimal value of a hyperparameter with two different metrics, the result is rarely the same. The important thing is that both metrics identify the same regions of interest.

Although individual analysis of hyperparameters is useful for understanding their impact on the model and identifying ranges of interest, the final search should not be done sequentially, as each hyperparameter interacts with the others. It is preferable to use grid search or random search to analyze several combinations of hyperparameters. More information about search strategies can be found in Machine learning with Python and Scikit-learn.

# Grid of evaluated hyperparameters
# ==============================================================================
param_grid = ParameterGrid(
    {'n_estimators': [150],
     'max_features': [5, 7, 9],
     'max_depth'   : [None, 3, 10, 20],
     'criterion'   : ['gini', 'entropy']
    }
)

# Loop to fit a model with each combination of hyperparameters
# ==============================================================================
results = {'params': [], 'oob_accuracy': []}

for params in param_grid:
    model = RandomForestClassifier(
        oob_score    = True,
        n_jobs       = -1,
        random_state = 123,
        **params
    )
    model.fit(X_train_prep, y_train)
    results['params'].append(params)
    results['oob_accuracy'].append(model.oob_score_)
    print(f"Model: {params} ✓")

# Results
# ==============================================================================
results = pd.DataFrame(results)
results = pd.concat([results, results['params'].apply(pd.Series)], axis=1)
results = results.sort_values('oob_accuracy', ascending=False)
results = results.drop(columns='params')
results.head(4)

Model: {'criterion': 'gini', 'max_depth': None, 'max_features': 5, 'n_estimators': 150} ✓
Model: {'criterion': 'gini', 'max_depth': None, 'max_features': 7, 'n_estimators': 150} ✓
Model: {'criterion': 'gini', 'max_depth': None, 'max_features': 9, 'n_estimators': 150} ✓
Model: {'criterion': 'gini', 'max_depth': 3, 'max_features': 5, 'n_estimators': 150} ✓
Model: {'criterion': 'gini', 'max_depth': 3, 'max_features': 7, 'n_estimators': 150} ✓
Model: {'criterion': 'gini', 'max_depth': 3, 'max_features': 9, 'n_estimators': 150} ✓
Model: {'criterion': 'gini', 'max_depth': 10, 'max_features': 5, 'n_estimators': 150} ✓
Model: {'criterion': 'gini', 'max_depth': 10, 'max_features': 7, 'n_estimators': 150} ✓
Model: {'criterion': 'gini', 'max_depth': 10, 'max_features': 9, 'n_estimators': 150} ✓
Model: {'criterion': 'gini', 'max_depth': 20, 'max_features': 5, 'n_estimators': 150} ✓
Model: {'criterion': 'gini', 'max_depth': 20, 'max_features': 7, 'n_estimators': 150} ✓
Model: {'criterion': 'gini', 'max_depth': 20, 'max_features': 9, 'n_estimators': 150} ✓
Model: {'criterion': 'entropy', 'max_depth': None, 'max_features': 5, 'n_estimators': 150} ✓
Model: {'criterion': 'entropy', 'max_depth': None, 'max_features': 7, 'n_estimators': 150} ✓
Model: {'criterion': 'entropy', 'max_depth': None, 'max_features': 9, 'n_estimators': 150} ✓
Model: {'criterion': 'entropy', 'max_depth': 3, 'max_features': 5, 'n_estimators': 150} ✓
Model: {'criterion': 'entropy', 'max_depth': 3, 'max_features': 7, 'n_estimators': 150} ✓
Model: {'criterion': 'entropy', 'max_depth': 3, 'max_features': 9, 'n_estimators': 150} ✓
Model: {'criterion': 'entropy', 'max_depth': 10, 'max_features': 5, 'n_estimators': 150} ✓
Model: {'criterion': 'entropy', 'max_depth': 10, 'max_features': 7, 'n_estimators': 150} ✓
Model: {'criterion': 'entropy', 'max_depth': 10, 'max_features': 9, 'n_estimators': 150} ✓
Model: {'criterion': 'entropy', 'max_depth': 20, 'max_features': 5, 'n_estimators': 150} ✓
Model: {'criterion': 'entropy', 'max_depth': 20, 'max_features': 7, 'n_estimators': 150} ✓
Model: {'criterion': 'entropy', 'max_depth': 20, 'max_features': 9, 'n_estimators': 150} ✓

# Best hyperparameters by out-of-bag accuracy
# ==============================================================================
print("--------------------------------------------------")
print("Best hyperparameters found (oob-accuracy)")
print("--------------------------------------------------")
print(results.iloc[0,0:])

--------------------------------------------------
Best hyperparameters found (oob-accuracy)
--------------------------------------------------
oob_accuracy       0.82
criterion       entropy
max_depth          20.0
max_features          7
n_estimators        150
Name: 22, dtype: object

# Grid of evaluated hyperparameters
# ==============================================================================
param_grid = {
    'n_estimators': [150],
    'max_features': [5, 7, 9],
    'max_depth'   : [None, 3, 10, 20],
    'criterion'   : ['gini', 'entropy']
}

# Grid search with cross-validation
# ==============================================================================
grid = GridSearchCV(
    estimator  = RandomForestClassifier(random_state = 123),
    param_grid = param_grid,
    scoring    = 'accuracy',
    n_jobs     = cpu_count() - 1,
    cv         = RepeatedKFold(n_splits=5, n_repeats=3, random_state=123), 
    refit      = True,
    verbose    = 0,
    return_train_score = True
)

grid.fit(X=X_train_prep, y=y_train)

# Results
# ==============================================================================
results = pd.DataFrame(grid.cv_results_)
results.filter(regex='(param*|mean_t|std_t)') \
    .drop(columns='params') \
    .sort_values('mean_test_score', ascending = False) \
    .head(4)

# Best hyperparameters found by cross-validation
# ==============================================================================
print("--------------------------------------------")
print("Best hyperparameters found by (cv)")
print("--------------------------------------------")
print(grid.best_params_, ":", grid.best_score_, grid.scoring)

--------------------------------------------
Best hyperparameters found by (cv)
--------------------------------------------
{'criterion': 'entropy', 'max_depth': None, 'max_features': 5, 'n_estimators': 150} : 0.8166666666666668 accuracy

Prediction and evaluation

Finally, the predictive ability of the final model is evaluated using the test set.

# Model with the best hyperparameters
# ==============================================================================
final_model = grid.best_estimator_
final_model

RandomForestClassifier(criterion='entropy', max_features=5, n_estimators=150,
                       random_state=123)

# Test error of the final model
# ==============================================================================
predictions = final_model.predict(X=X_test_prep)
predictions[:10]

array([1, 1, 1, 1, 1, 1, 0, 0, 0, 0])

conf_matrix = confusion_matrix(y_true=y_test, y_pred=predictions)
accuracy = accuracy_score(y_true=y_test, y_pred=predictions, normalize=True)
print("Confusion matrix")
print("-------------------")
print(conf_matrix)
print("")
print(f"Test accuracy is: {100 * accuracy} % \n")
fig, ax = plt.subplots(figsize=(3, 3))
ConfusionMatrixDisplay(conf_matrix).plot(ax=ax)

print(
    classification_report(
        y_true = y_test,
        y_pred = predictions
    )
)

Confusion matrix
-------------------
[[47  3]
 [17 33]]

Test accuracy is: 80.0 % 

              precision    recall  f1-score   support

           0       0.73      0.94      0.82        50
           1       0.92      0.66      0.77        50

    accuracy                           0.80       100
   macro avg       0.83      0.80      0.80       100
weighted avg       0.83      0.80      0.80       100

After optimizing the hyperparameters, an accuracy of 78% is achieved.

Probabilistic prediction

Most implementations of Random Forest, including scikit-learn, allow probability prediction for classification problems. It is important to understand how these values are calculated to interpret and use them correctly.

In the previous example, applying .predict() returns $1$ (high sales) or $0$ (low sales) for each test observation. However, there is no information about the confidence with which the model makes this assignment. With .predict_proba(), instead of a classification, you get the probability with which the model considers that each observation may belong to each of the classes.

# Probability prediction
# ==============================================================================
predictions = final_model.predict_proba(X=X_test_prep)
predictions[:5, :]

array([[0.21333333, 0.78666667],
       [0.28      , 0.72      ],
       [0.14      , 0.86      ],
       [0.27333333, 0.72666667],
       [0.20666667, 0.79333333]])

The result of .predict_proba() is an array with one row per observation and as many columns as there are classes in the response variable. The value in the first column corresponds to the probability, according to the model, that the observation belongs to class 0, and so on. The probability value shown for each prediction corresponds to the fraction of observations of each class in the terminal nodes reached by the predicted observation across all trees.

By default, .predict() assigns each new observation to the class with the highest probability (in case of a tie, it is assigned randomly). However, this may not always be the desired behavior.

# Classification using the class with the highest probability
# ==============================================================================
df_predictions = pd.DataFrame(data=predictions, columns=['0', '1'])
df_predictions['classification_default_0.5'] = np.where(df_predictions['0'] > df_predictions['1'], 0, 1)
df_predictions.head(3)

Suppose the following scenario: the Christmas campaign is approaching and the owners of the chain want to double the stock of items in those stores expected to have high sales. Since transporting this material to the stores is costly, the manager wants to limit this strategy only to stores for which there is high confidence that they will achieve high sales.

If probabilities are available, a specific cutoff point can be set, for example, considering only as class $1$ (high sales) those stores whose prediction for this class is greater than 0.9 (90%). In this way, the final classification is better tailored to the needs of the use case.

# Final classification using a threshold of 0.8 for class 1.
# ==============================================================================
df_predictions['classification_custom_0.8'] = np.where(df_predictions['1'] > 0.9, 1, 0)
df_predictions.iloc[4:10, :]

How much should you trust these probabilities?

It is very important to keep in mind the difference between the "view" the model has of the world and the real world. Everything a model knows about the real world is what it has learned from the training data and, therefore, it has a limited "view". For example, suppose that in the training data, all stores in urban areas Urban='Yes' have high sales regardless of the value of the other predictors. When the model tries to predict a new observation, if it is in an urban area, it will classify the store as high sales with 100% confidence. However, this does not mean it is unequivocally true; there could be stores in urban areas that do not have high sales but, since they are not present in the training data, the model does not consider this possibility.

With this in mind, the probabilities generated by the model should be considered as the model's confidence, from its limited perspective, when making predictions. But not as the real-world probability that this is the case.

Predictor importance

feature_importance = pd.DataFrame(
    {'predictor': X_train_prep.columns,
     'importance': final_model.feature_importances_}
    )
print("Predictor importance in the model")
print("-------------------------------------------")
feature_importance.sort_values('importance', ascending=False)

Predictor importance in the model
-------------------------------------------

Permutation

Hide

importance = permutation_importance(
    estimator    = final_model,
    X            = X_train_prep,
    y            = y_train,
    n_repeats    = 5,
    scoring      = 'neg_root_mean_squared_error',
    n_jobs       = cpu_count() - 1,
    random_state = 123,
    )

# Store the results (mean and std) in a dataframe
df_importance = pd.DataFrame(
    {k: importance[k] for k in ['importances_mean', 'importances_std']}
    )
df_importance['feature'] = X_train_prep.columns
df_importance.sort_values('importances_mean', ascending=False)

	importances_mean	importances_std	feature
9	0.448637	0.026930	Price
1	0.338378	0.012934	ShelveLoc_Good
7	0.295299	0.011481	Advertising
5	0.257405	0.020228	CompPrice
0	0.254924	0.018636	ShelveLoc_Bad
10	0.176646	0.011373	Age
6	0.130213	0.019439	Income
8	0.084207	0.015574	Population
11	0.060537	0.031654	Education
4	0.034641	0.028284	US_Yes
3	0.000000	0.000000	Urban_Yes
2	0.000000	0.000000	ShelveLoc_Medium

Hide

# Plot
# ==============================================================================
fig, ax = plt.subplots(figsize=(3.5, 4))
df_importance = df_importance.sort_values('importances_mean', ascending=True)
ax.barh(
    df_importance['feature'],
    df_importance['importances_mean'],
    xerr=df_importance['importances_std'],
    align='center',
    alpha=0
)
ax.plot(
    df_importance['importances_mean'],
    df_importance['feature'],
    marker="D",
    linestyle="",
    alpha=0.8,
    color="r"
)
ax.set_title('Predictor importance (train)')
ax.set_xlabel('Increase in error after permutation');

Both strategies identify Price, Advertising, and ShelveLoc as the most influential predictors, according to the training data.