[Python] How to train LightGBM to directly optimize for period-over-period performance...

I'm trying to build a model to predict the outcome (5 classes at like [70%, 20%, 1%, 9%]) of individual events for different entities over several time periods (e.g., years). For each entitiy, the class mean averages vary like +/- 5%. While event-level accuracy is important, my overall goal is to create a model whose predictions are stable and predictive of future performance. The main hang-up is while the primary goal is period-to-period correlation between expected entity mean for T to real entity mean for T+1, I also want the expected outcome, on an outcome to outcome basis, to be in the right ballpark still.

I want the model's average predicted probability for a given entity in one period to have the highest possible correlation with that entity's actual performance in the next period.

Currently, I'm using hyperopt to tune a LGBM classifier. I have an objective function for hyperopt that evaluates this period-over-period correlation.

I use K-Fold cross-validation to get out-of-fold predictions for every event in my dataset.

Aggregate these predictions and the actual outcomes to the entity-period level.

For each entity, shift the aggregated predictions to align last period's prediction with this period's actual outcome.

Calculate a weighted Pearson correlation between the projected and actual outcomes. Not in love with the weighted Pearson, but I think is passes for now.

In theory, I think this works well for finding the best hyperparameters for my specific goal.

The issue is that while I'm tuning for the correlation metric, the LGBMClassifier itself is still just training to minimize standard logloss during the model.fit() process. There's a mismatch between my final evaluation metric and the model's training objective.

Is it possible to define a custom training objective for LightGBM that directly optimizes (or closely approximates) this kind of aggregate period-over-period correlation metric?

I understand that custom objectives in LightGBM typically require providing the gradient and hessian for each individual sample, which seems impossible for a metric that depends on groupby() and shift() operations across the entire dataset. Are there other approaches that would be better suited for this problem than simple logloss?

Here is my code:


X = df[['feature_1', 'feature_2', 'feature_3', 'feature_4', 'feature_5_flag', 'time_period']]
y = df['outcome']

monotone_constraints = [0, 0, 0, -1, 0, 0]

space = {
'n_estimators': hp.quniform('n_estimators', 100, 1000, 100),
'learning_rate': hp.loguniform('learning_rate', np.log(0.01), np.log(0.1)),
'max_depth': hp.quniform('max_depth', 3, 10, 1)
}

df['entity_period_id'] = df['entity_id'].astype(str) + df['time_period'].astype(str)

def objective(space):
"""
Objective function to be minimized by hyperopt.
This calculates the metric I want the model to learn directly.
"""
N_SPLITS = 5
skf = StratifiedKFold(n_splits=N_SPLITS, shuffle=True, random_state=6)
oof_preds = np.zeros(len(df))

for train_idx, val_idx in skf.split(X, y):
X_train, y_train = X.iloc[train_idx], y.iloc[train_idx]
X_val, _ = X.iloc[val_idx], y.iloc[val_idx]

model = lgb.LGBMClassifier(
n_estimators=int(space['n_estimators']),
learning_rate=space['learning_rate'],
max_depth=int(space['max_depth']),
random_state=6, n_jobs=-1,
monotone_constraints=monotone_constraints
)
model.fit(X_train, y_train)
fold_preds = model.predict_proba(X_val)[:, 1]
oof_preds[val_idx] = fold_preds

eval_df = df.copy()
eval_df['predicted_prob'] = oof_preds
entity_means = (
eval_df.groupby('entity_period_id')
.agg(
avg_predicted_prob=('predicted_prob', 'mean'),
actual_outcome_rate=('outcome', 'mean'),
event_count=('outcome', 'count')
).reset_index()
)

entity_means['entity_id'] = entity_means['entity_period_id'].str[:-4]
entity_means['period'] = entity_means['entity_period_id'].str[-4:].astype(int)
entity_means = entity_means.sort_values(by=['entity_id', 'period'])

entity_means['projected_value'] = entity_means.groupby('entity_id')['avg_predicted_prob'].shift(1)
entity_means['prev_period_event_count'] = entity_means.groupby('entity_id')['event_count'].shift(1)
entity_means.dropna(subset=['projected_value', 'actual_outcome_rate', 'prev_period_event_count'], inplace=True)

def weighted_corr(x, y, w):
if len(x) < 2: return 0
w, x, y = np.array(w), np.array(x), np.array(y)
w_mean_x, w_mean_y = np.average(x, weights=w), np.average(y, weights=w)
cov_xy = np.average((x - w_mean_x) * (y - w_mean_y), weights=w)
var_x, var_y = np.average((x - w_mean_x)**2, weights=w), np.average((y - w_mean_y)**2, weights=w)
if var_x == 0 or var_y == 0: return 0
return cov_xy / np.sqrt(var_x * var_y)

if len(entity_means) < 2:
corr = 0
else:
corr = weighted_corr(
entity_means['projected_value'],
entity_means['actual_outcome_rate'],
entity_means['prev_period_event_count']
)
if np.isnan(corr): corr = 0

return {'loss': -corr, 'status': STATUS_OK}

trials = Trials()
best = fmin(fn=objective, space=space, algo=tpe.suggest, max_evals=100, trials=trials)

print(best)

lgb_model = lgb.LGBMClassifier(
n_estimators=int(best['n_estimators']),
learning_rate=best['learning_rate'],
max_depth=int(best['max_depth']),
random_state=6,
n_jobs=-1,
monotone_constraints=monotone_constraints
)

lgb_model.fit(X, y)

Continue reading...