Ensemble of machine learning models in scikit-learn

There will be multiple steps to creating an ensemble model.

Start by creating the two models individually. For the first model, split the data by group and train two individual models, then join the two models together in a function. For the second model, the data can be left in its entirety (aside from removing testing data). Then, create another method to join the other two models into one ensemble model.

To demonstrate, I'll start by importing the necessary modules and loading in the dataframe:

import numpy as np
import pandas as pd
from sklearn.ensemble import RandomForestRegressor

data_str = """group,feature_1,feature_2,year,dependent_variable
group_a,12,19,2010,0.4
group_a,11,13,2011,0.9
group_a,10,5,2012,1.2
group_a,16,9,2013,3.2
group_b,8,29,2010,0.6
group_b,9,33,2011,0.1 
group_b,111,15,2012,2.1 
group_b,16,19,2013,12.2"""

data_list = [row.split(",") for row in data_str.split("\n")]
data = pd.DataFrame(data_list[1:], columns = data_list[0])

train = data.loc[data["year"] != "2013"]
test = data.loc[data["year"] == "2013"]

This will be using a RandomForestRegressor ensemble model, but any regression model can be used. In addition, it should be noted that the dataframe used here differs from the given dataframe in that this dataframe has its rows indexed from 0 rather than being indexed by group, and group is instead a column within the dataframe.

To construct the first model:

1. split the data into data for group a and for group b
2. train two independent models
3. join the models

The first two steps are done below:

# Splitting Data
train_a = train.loc[train["group"] == "group_a"]
train_b = train.loc[train["group"] == "group_b"]
test_a = test.loc[test["group"] == "group_a"]
test_b = test.loc[test["group"] == "group_b"]

# Training Two Models
model_a = RandomForestRegressor()
model_a.fit(train_a.drop(["dependent_variable", "year", "group"], axis = "columns"), train_a.dependent_variable)
model_b = RandomForestRegressor()
model_b.fit(train_b.drop(["dependent_variable", "year", "group"], axis = "columns"), train_b.dependent_variable)

Then, their predict methods can be joined together:

def individual_predictor(group, feature_1, feature_2):
    if group == "group_a": return model_a.predict([[feature_1, feature_2]])[0]
    elif group == "group_b": return model_b.predict([[feature_1, feature_2]])[0]

This will take in a group and two features individually and return the prediction. This can be adapted to whatever input and output type is necessary.

To create the second model, leave the data as whole and only train one model, which also removes the necessity to join the models:

model = RandomForestRegressor()
model.fit(train.drop(["dependent_variable", "year", "group"], axis = "columns"), train.dependent_variable)

Finally, you can join the models together into an ensemble model by averaging the result of their predict methods:

def ensemble_predict(group, feature_1, feature_2):
    return (individual_predictor(group, feature_1, feature_2) + model.predict([[feature_1, feature_2]])[0]) / 2

Again, this takes in a group and two features then returns the result. This will likely need to be adapted into another format, such as taking in a list of list of inputs and outputting a list of predictions.

This one uses 2 regressors, RandomForestRegressor, and GradientBoostingRegressor.

I Added 2013 data for r2_score calculation, it must be more than 1. Also added data from other years. Copy the text and save to txt file.

First we process the data file, separate the train and test by dataframe manipulation. We then create a model for each regressor with model 1.1 and 1.2 for group "a" and "b" respectively. Then model 2 for all data. After creating the model we then save it to disk for later processing.

After the models are created we then make predictions using all the test data and a single data. Metrics r2_square and MAE are also printed.

The last part is testing the model file by loading it and let it predict from a test. Predictions from model in memory and disk should be the same. Also there is a sample input types and how to use it in custom prediction function.

See also the docstring and comments in the code on how this works.

data.txt

group        feature_1        feature_2       year            dependent_variable
group_a         12               19           2010               0.4
group_a          7               15           2010               1.5
group_a         11               13           2011               0.9
group_a          8               8            2011               2.1
group_a         10               5            2012               1.2
group_a         11               9            2012               2.6
group_a         16               9            2013               3.2
group_a         8               10            2013               2.6
group_b         8               29            2010               0.6
group_b         11              18            2010               1.5
group_b         9               33            2011               0.1 
group_b         20              15            2011               2.8 
group_b         111             15            2012               2.1 
group_b         99              10            2012               3.6
group_b         16              19            2013               12.2
group_b         4                8            2013               5.1

Code

myensemble.py

"""sklearn ensemble modeling.

Dependencies:
    * sklearn
    * pandas
    * numpy

References:
    * https://scikit-learn.org/stable/modules/classes.html?highlight=ensemble#module-sklearn.ensemble
    * https://pandas.pydata.org/docs/user_guide/indexing.html
"""


from typing import List, Union, Optional
import pickle  # for saving file to disk

from sklearn.ensemble import RandomForestRegressor
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.metrics import r2_score
from sklearn.metrics import mean_absolute_error
import pandas as pd
import numpy as np


def make_model(regressor, regname: str, modelfn: str, dfX: pd.DataFrame, dfy: pd.DataFrame):
    """Creates a model.

    Args:
        regressor: Can be RandomForestRegressor or GradientBoostingRegressor.
        regname: Regressor name.
        dfX: The features in pandas dataframe.
        dfy: The target in pandas dataframe.

    Returns:
        Model
    """
    X = dfX.to_numpy()
    y = dfy.to_numpy()
    model = regressor(random_state=0)
    model.fit(X, y)

    # Save model.
    with open(f'{regname}_{modelfn}', 'wb') as f:
        pickle.dump(model, f)

    return model


def get_prediction(model, test: Union[List, pd.DataFrame, np.ndarray]) -> Optional[np.ndarray]:
    """Returns prediction based on model and test input or None.
    """
    if isinstance(test, List) or isinstance(test, np.ndarray):
        return model.predict([test])
    if isinstance(test, pd.DataFrame):
        return model.predict(np.array(test))
    return None


def model_and_prediction(df: pd.DataFrame, regressor, regname: str, modelfn: str):
    """Build model and show prediction and metrics.

    To build a model we need a training data X with features
    and data y with target or dependent values.

    Args:
        df: A dataframe.
        regressor: Can be RandomForestRegressor or GradientBoostingRegressor.
        regname: The regressor name.
        modelfn: The filename where model will be saved to disk.

    Returns:
        None
    """
    features = ['feature_1', 'feature_2']

    # 1. Get the train dataframe
    train = df.loc[df.year != 2013]  # exclude 2013 in training data
    train_feature = train[features]  # select the features column
    train_target = train.dependent_variable  # select the dependent column

    model = make_model(regressor, regname, modelfn, train_feature, train_target)

    # 2. Get the test dataframe
    test = df.loc[df.year == 2013]  # only include 2013 in test data
    test_feature = test[features]
    test_target = test.dependent_variable

    # 3. Get the prediction from all rows in test feature. See step 5
    # for single data prediction.
    prediction: np.ndarray = model.predict(np.array(test_feature))

    print(f'test feature:\n{np.array(test_feature)}')
    print(f'test prediction: {prediction}')  # prediction[0] ...
    print(f'test target: {np.array(test_target)}')

    # 4. metrics
    print(f'r2_score: {r2_score(test_target, prediction)}')
    print(f'mean_absolute_error: {mean_absolute_error(test_target, prediction)}\n')

    # 5. Get prediction from the first row of test features.
    prediction_1: np.ndarray = model.predict(np.array(test_feature.iloc[[0]]))
    print(f'1st row test:\n{test_feature.iloc[[0]]}')
    print(f'1st row test prediction array: {prediction_1}')
    print(f'1st row test prediction value: {prediction_1[0]}\n')  # get the element value


def main():
    datafn = 'data.txt'
    df = pd.read_fwf(datafn)
    print(df.to_string(index=False))

    # A. Create models for each type of regressor.
    regressors = [(RandomForestRegressor, 'RandomForrest'),
                  (GradientBoostingRegressor, 'GradientBoosting')]

    for (r, name) in regressors:
        print(f'::: Regressor: {name} :::\n')

        # Model 1 using group_a
        print(':: MODEL 1.1 ::')
        grp = 'group_a'
        modelfn = f'{grp}.pkl'  # filename of model to be save to disk
        dfa = df.loc[df.group == grp]  # select group
        model_and_prediction(dfa, r, name, modelfn)

        # Model 1 using group_b
        print(':: MODEL 1.2 ::')
        grp = 'group_b'
        modelfn = f'{grp}.pkl'
        dfb = df.loc[df.group == grp]
        model_and_prediction(dfb, r, name, modelfn)

        # Model 2 using group a and b
        print(':: MODEL 2 ::')
        grp = 'group_ab'
        modelfn = f'{grp}.pkl'
        dfab = df.loc[(df.group == 'group_a') | (df.group == 'group_b')]
        model_and_prediction(dfab, r, name, modelfn)

    # B. Test saved model file prediction.
    print('::: Prediction from loaded model :::')
    mfn = 'GradientBoosting_group_ab.pkl'
    print(f'model: gradient boosting model 2, {mfn}')

    with open(mfn, 'rb') as f:
        loaded_model = pickle.load(f)

    # test: group_b  4  8  2013  5.1    
    test = [4, 8]
    prediction = loaded_model.predict([test])
    print(f'test: {test}')
    print(f'prediction: {prediction[0]}\n')

    # C. Use get_prediction().

    # input from list
    test = [4, 8]
    prediction = get_prediction(loaded_model, test)
    print(f'test from list input:\n{test}')
    print(f'prediction from get_prediction() with list input: {prediction}\n')

    # input from dataframe
    testdata = {
        'feature_1': [8, 12],
        'feature_2': [19, 15],
    }
    testdf = pd.DataFrame(testdata)
    testrow = testdf.iloc[[0]]  # first row [8, 19]
    prediction = get_prediction(loaded_model, testrow)
    print(f'test from df input:\n{testrow}')
    print(f'prediction from get_prediction() with df input: {prediction}\n')

    testrow = testdf.iloc[[1]]  # second row [12, 15]
    prediction = get_prediction(loaded_model, testrow)
    print(f'test from df input:\n{testrow}')
    print(f'prediction from get_prediction() with df input: {prediction}\n')

    # input from numpy
    test = [8, 9]
    testnp = np.array(test)
    prediction = get_prediction(loaded_model, testnp)
    print(f'test from numpy input:\n{testnp}')
    print(f'prediction from get_prediction() with numpy input: {prediction}\n')


if __name__ == '__main__':
    main()

Output

  group  feature_1  feature_2  year  dependent_variable
group_a         12         19  2010                 0.4
group_a          7         15  2010                 1.5
group_a         11         13  2011                 0.9
group_a          8          8  2011                 2.1
group_a         10          5  2012                 1.2
group_a         11          9  2012                 2.6
group_a         16          9  2013                 3.2
group_a          8         10  2013                 2.6
group_b          8         29  2010                 0.6
group_b         11         18  2010                 1.5
group_b          9         33  2011                 0.1
group_b         20         15  2011                 2.8
group_b        111         15  2012                 2.1
group_b         99         10  2012                 3.6
group_b         16         19  2013                12.2
group_b          4          8  2013                 5.1

::: Regressor: RandomForrest :::

:: MODEL 1.1 ::
test feature:
[[16  9]
 [ 8 10]]
test prediction: [1.811 2.186]
test target: [3.2 2.6]
r2_score: -10.67065000000004
mean_absolute_error: 0.9015000000000026

1st row test:
   feature_1  feature_2
6         16          9
1st row test prediction array: [1.811]
1st row test prediction value: 1.8109999999999986

:: MODEL 1.2 ::
test feature:
[[16 19]
 [ 4  8]]
test prediction: [2.116 2.408]
test target: [12.2  5.1]
r2_score: -3.3219170799444546
mean_absolute_error: 6.388

1st row test:
    feature_1  feature_2
14         16         19
1st row test prediction array: [2.116]
1st row test prediction value: 2.116000000000001

:: MODEL 2 ::
test feature:
[[16  9]
 [ 8 10]
 [16 19]
 [ 4  8]]
test prediction: [2.425 2.145 1.01  1.958]
test target: [ 3.2  2.6 12.2  5.1]
r2_score: -1.3250936994738867
mean_absolute_error: 3.8905000000000016

1st row test:
   feature_1  feature_2
6         16          9
1st row test prediction array: [2.425]
1st row test prediction value: 2.4249999999999985

::: Regressor: GradientBoosting :::

:: MODEL 1.1 ::
test feature:
[[16  9]
 [ 8 10]]
test prediction: [2.59996945 2.21271005]
test target: [3.2 2.6]
r2_score: -1.8335008778823685
mean_absolute_error: 0.4936602458577084

1st row test:
   feature_1  feature_2
6         16          9
1st row test prediction array: [2.59996945]
1st row test prediction value: 2.59996945439128

:: MODEL 1.2 ::
test feature:
[[16 19]
 [ 4  8]]
test prediction: [1.99807124 2.63511811]
test target: [12.2  5.1]
r2_score: -3.3703627491779713
mean_absolute_error: 6.333405322236132

1st row test:
    feature_1  feature_2
14         16         19
1st row test prediction array: [1.99807124]
1st row test prediction value: 1.9980712422931164

:: MODEL 2 ::
test feature:
[[16  9]
 [ 8 10]
 [16 19]
 [ 4  8]]
test prediction: [3.60257456 2.26208935 0.402739   2.10950224]
test target: [ 3.2  2.6 12.2  5.1]
r2_score: -1.538939968014979
mean_absolute_error: 3.882060991360607

1st row test:
   feature_1  feature_2
6         16          9
1st row test prediction array: [3.60257456]
1st row test prediction value: 3.6025745572622014

::: Prediction from loaded model :::
model: gradient boosting model 2, GradientBoosting_group_ab.pkl
test: [4, 8]
prediction: 2.1095022367629728

test from list input:
[4, 8]
prediction from get_prediction() with list input: [2.10950224]

test from df input:
   feature_1  feature_2
0          8         19
prediction from get_prediction() with df input: [0.50307204]

test from df input:
   feature_1  feature_2
1         12         15
prediction from get_prediction() with df input: [1.46058714]

test from numpy input:
[8 9]
prediction from get_prediction() with numpy input: [2.30007317]

Firstly, create models using time series algorithms(use only date variable and dependent variable), fbprophet (uses features + date + dependent variable), Treebased regression algorithms like CatBoost/XGBoost/LightGBM (use features + date + dependent variable).

Using each of the mentioned algorithms create models for each group (bottom up approach). Different models will peform well for different groups. Take a weighted mean based on models' performance. Assume, group_a predictions perform best with Catboost, then with fbprophet and then with exponential moving average, use weights in proportion to accuracies derived from these models.

You can aggregate results of group level models to get aggregated results. You can also create separate models on aggregated data(summing up on year).

If I understand the last line of your question correctly, in context with the years, you are looking to capture the trend in a given calendar year via model 1 and capture trend across multiple years via model 2. And model 2 in where it could be an issue because you mentioned scikit-learn compatible models.

So I'll try to explain the approach I would take.

Model 1 is pretty straight forward, it is a regression problem so selecting the best regression model should not be an issue. You can find that by seeing results in a given calendar year.

Model 2 is where you would like to capture time-series features, kind of like YoY sort of thing. While there isn't any model in SKLearn to directly capture the time parameter like ARIMA or RNNs do, there are ways to use SKLearn Models to do the forecasting. And a lot of it depends on feature engineering. You could use features 1 and 2, sort them, shift them and then take a diff to create new features, say 1a and 2a, which could then be used with any regression model. These new features would capture the time essence. I could write a lengthy post on that here, but I think you'll find this link much better written.

Now coming to ensembling the 2 models together. As this is a regression problem, the best way I feel would be to assign weightages to the output of both these models, lets say alpha for model 1 and beta for model 2. Treat alpha and beta as hyperparameters. Tune them using the data.

This should make a pretty good ensemble with SKLearn models.