Turning regression problem into "classification + regression"

Question

Say I have a regression problem where I'd like to predict values ranging from 0 - 100.000 based on some predictors. A single XGBoost model achieves decent overall mean performance (measured by MAE) using 5-fold CV.

Then I looked at specific sub-ranges of the [0,100.000] response interval to see where the model performed better or worse. It turned out that especially in the beginning and the end, i.e. [0, 5.000] and [80.000, 100.000], the model performed worst.

So I took a second step and partitioned the response interval into several non-overlapping sub-intervals and trained a separate regression model for each, resulting in significant better performance when compared to the predictions coming from the first model.

Now, what I have is a set of trained regression models for several sub-intervals of the response with good predictive performance for the entire range of the response. I am wondering how to apply these on unseen cases:

Obviously, for these I do not know the actual value (as this is what I want to predict), hence I do not know which of the sub-models to apply to get the prediction.

I was thinking of training another model - this time a classifier - that first predicts the interval which the case is belonging to and then secondly applying the regression model that was specifically trained on that particular interval to get the final prediction of the response.

This approach would be some kind of a "stacked" model, where the first one (the classifier) makes a rough decision on the range of the response and based on this prediction "routes" the cases to the particular regression model that predicts the exact value.

I have to questions:

1. Is this a common approach? Are there better ones? (I was searching a lot in the internet but did not find anything on that matter.)
2. I imagine training/testing this "hierarchical" model would be complicated and time-consuming to implement, especially using CV. I could imagine that it could be done using the "pipeline" framework from scikit-learn. Are there already some reference implementations I could use?

Thanks in advance for your comments/suggestions.

Cheers Chris

Answer 1

As you well noticed there is no way to know the bin in wich an unseen data's target value will be. So what you can do is to train a model that splits/clusters your data and then run a model on each cluster/group. This is possible since the first model will be able to make Inference on aun unseen x value for next running the model that corresponds to that group.

Unlike your first approach It does not take anything about your target, but is only clustering similar points so that hopefully, individual models could work better than a single model on all the dataset.

You can also try to scale the target with Standard transformation, MixMax or log so that the target features is more centered arround its mean, this in some case might help.

Below you can find an example using Boston Housing dataset:

import pandas as pd
import numpy as np
from sklearn.datasets import fetch_openml
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.cluster import KMeans
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.compose import TransformedTargetRegressor

boston = fetch_openml("boston")

X = boston["data"]
y = boston["target"]
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size= .3, random_state= 42)

kmeans = KMeans(n_clusters= 3, random_state= 42)

cluster_model = Pipeline([("scaler", StandardScaler()), ("model", kmeans)]).fit(X_train)

X_train["cluster"] = cluster_model.predict(X_train)

dict_models = {}
for k in set(cluster_model["model"].labels_):
    X_train_tmp = X_train[X_train["cluster"] == k].copy()
    X_train_tmp.drop("cluster", axis = 1, inplace = True)
    y_train_tmp = y_train[X_train_tmp.index].copy()
    n_samples = X_train_tmp.shape[0]
    
    model = GradientBoostingRegressor(random_state= 42).fit(X_train_tmp, y_train_tmp)

    score = cross_val_score(estimator= model, X = X_train_tmp, y =  y_train_tmp, scoring= "r2", cv = 3).mean()
    print(f"Model trained with: {n_samples} samples\avg cv score: {round(score,3)}")
    
    
    dict_models[k] = model
print(f"-"*50)
model = GradientBoostingRegressor(random_state= 42).fit(X_train, y_train)
score = cross_val_score(estimator= model, X = X_train, y =  y_train, scoring= "r2", cv = 3).mean()
print(f"Model trained with all the data avg cv score: {round(score,3)}")
model = TransformedTargetRegressor(regressor= GradientBoostingRegressor(), func= np.log1p, inverse_func= np.expm1).fit(X_train, y_train)
score = cross_val_score(estimator= model, X = X_train, y =  y_train, scoring= "r2", cv = 3).mean()
print(f"Model transformed target trained with all the data avg cv score: {round(score,3)}")

# evaluate on unseen data observation
sample = X_test.sample()
model_id = cluster_model.predict(sample)[0]
# this dictionary may be serialized and saved as a pickle object for production
dict_models[model_id].predict(sample)

So this is no a standard way of solving a supervised problem,so I'm not giving a recipe but an alternative that might work for your problem.

Finally consider those points:

1. It is important to take into consideration the size of each cluster since, a model trained in few observations will hardly generalize

2. Number of clusters, and cluster model itself are now hyper parameters that should be optimized via CV

3. There is no guarantee that individual models will perform better that a model with all the data

Hope it helps!