LEARN · GRADIENT BOOSTING WITH XGBOOST
Why a good score is not a diagnosis
A boosted-tree model can post an attractive validation score while still being wrong for the reasons that matter: leakage, an unrepresentative split, a brittle interaction, or one small segment carrying most of the loss.
The debugging workflow below moves from training dynamics to model structure, then from global explanations to individual predictions, and finally to residual slices. The example is a synthetic delivery-time regression problem, so it is safe to run without downloading data.
Use a current XGBoost 3.x environment. The plotting helper for trees also needs Graphviz available on the operating system, not only the Python package. XGBoost’s current plotting API uses tree_idx; num_trees is deprecated in 3.x.
python -m pip install --upgrade "xgboost>=3,<4" shap scikit-learn pandas matplotlib graphviz
1. Create data with a deliberate blind spot
The target is delivery time in minutes. A diagnostic column, delivery_zone, is kept outside the feature matrix. Remote deliveries have an extra delay and higher variance, creating a realistic failure mode: the aggregate model can look useful while a minority segment is systematically underpredicted.
import numpy as np import pandas as pd import matplotlib.pyplot as plt import shap import xgboost as xgb from sklearn.inspection import PartialDependenceDisplay from sklearn.metrics import mean_absolute_error, mean_squared_error from sklearn.model_selection import train_test_split rng = np.random.default_rng(7) n = 12_000 zone = rng.choice( ["urban", "suburban", "remote"], size=n, p=[0.62, 0.28, 0.10], ) distance_km = rng.gamma(2.2, 3.0, size=n).clip(0.5, 30) traffic_index = rng.beta(2.5, 2.0, size=n) warehouse_load = rng.integers(20, 250, size=n) items = rng.poisson(3.0, size=n) + 1 hour = rng.integers(0, 24, size=n) is_weekend = rng.integers(0, 2, size=n) promo = rng.integers(0, 2, size=n) X = pd.DataFrame( { "distance_km": distance_km, "traffic_index": traffic_index, "warehouse_load": warehouse_load, "items": items, "hour_sin": np.sin(2 * np.pi * hour / 24), "hour_cos": np.cos(2 * np.pi * hour / 24), "is_weekend": is_weekend, "promo": promo, } ) zone_effect = np.select( [zone == "suburban", zone == "remote"], [5.0, 18.0 + 0.7 * distance_km * traffic_index], default=0.0, ) noise_scale = np.select( [zone == "suburban", zone == "remote"], [5.0, 12.0], default=3.5, ) y = ( 14.0 + 1.05 * distance_km + 8.0 * traffic_index + 0.03 * warehouse_load + 1.4 * np.sqrt(items) + 2.5 * is_weekend - 1.5 * promo + 3.0 * np.sin(2 * np.pi * (hour - 8) / 24) + zone_effect + rng.normal(0, noise_scale, size=n) ) metadata = pd.DataFrame({"delivery_zone": zone}, index=X.index) train_idx, valid_idx = train_test_split( X.index, test_size=0.25, random_state=7, stratify=metadata["delivery_zone"], ) X_train, X_valid = X.loc[train_idx], X.loc[valid_idx] y_train, y_valid = y[train_idx], y[valid_idx] metadata_valid = metadata.loc[valid_idx]
2. Capture learning curves during training
Put training data first and validation data last in eval_set. XGBoost names them validation_0 and validation_1; the last evaluation set is used for early stopping. The estimator’s evals_result() method returns the metric history, while prediction automatically uses best_iteration after early stopping.
model = xgb.XGBRegressor( n_estimators=1_500, learning_rate=0.03, max_depth=6, min_child_weight=5, subsample=0.85, colsample_bytree=0.85, reg_lambda=5.0, objective="reg:squarederror", eval_metric="rmse", tree_method="hist", early_stopping_rounds=75, random_state=7, n_jobs=-1, ) model.fit( X_train, y_train, eval_set=[(X_train, y_train), (X_valid, y_valid)], verbose=False, ) history = model.evals_result() train_rmse = history["validation_0"]["rmse"] valid_rmse = history["validation_1"]["rmse"] fig, ax = plt.subplots(figsize=(9, 5)) ax.plot(train_rmse, label="train") ax.plot(valid_rmse, label="validation") ax.axvline(model.best_iteration, linestyle="--", label="best iteration") ax.set( xlabel="Boosting round", ylabel="RMSE", title="Training and validation learning curves", ) ax.legend() fig.tight_layout() plt.show()
Read the shape before tuning:
-
Overfit: Training error keeps falling while validation error bottoms out and rises. Reduce depth, add regularization, increase data, or stop earlier.
-
Underfit: Both curves remain high and close. Add useful features, allow more capacity, or train longer with a suitable learning rate.
-
Broken split: The gap is extreme from the first rounds, validation is implausibly easy, or the curve is unstable. Check time ordering, entity overlap, duplicate rows, target leakage, sample size, and whether train and validation represent the same deployment population.
Early stopping limits wasted rounds; it does not repair a bad validation design.
3. Compare three incompatible meanings of importance
XGBoost defines weight as split count and gain as the average improvement from splits using a feature. Cover measures the average coverage of those splits; for this squared-error example, it is proportional to how many training rows reach them. These answer different questions, so disagreement is expected.
for importance_type in ["weight", "gain", "cover"]: fig, ax = plt.subplots(figsize=(8, 5)) xgb.plot_importance( model, ax=ax, importance_type=importance_type, max_num_features=8, show_values=False, title=f"Feature importance: {importance_type}", ) fig.tight_layout() plt.show()
Interpret the contrasts:
-
High weight but modest gain suggests a feature is a frequent fine-tuner.
-
Low weight but high gain suggests a rare, decisive splitter.
-
High cover often points to upstream splits affecting many rows.
-
A feature absent from the chart may simply be unused; XGBoost omits zero-importance features from its score dictionary.
Do not turn any of these charts into a causal ranking.
4. Read an actual tree and trace one row
A tree picture is useful for orientation, but a DataFrame dump is easier to inspect programmatically. XGBoost’s trees_to_dataframe() parses tree dumps into rows for split and leaf nodes.
booster = model.get_booster()
tree_frame = booster.trees_to_dataframe()
fig, ax = plt.subplots(figsize=(18, 10))
xgb.plot_tree(
model,
tree_idx=0,
with_stats=True,
ax=ax,
)
fig.tight_layout()
plt.show()
print(
tree_frame.loc[
tree_frame["Tree"] == 0,
["ID", "Feature", "Split", "Yes", "No", "Missing", "Gain", "Cover"],
].head(12)
)
To inspect the exact route for one row, follow the Yes or No child at every numeric split. Leaf rows store their contribution in the Gain column of this DataFrame representation.
def trace_tree( tree_frame: pd.DataFrame, row: pd.Series, tree_idx: int = 0, ) -> pd.DataFrame: nodes = tree_frame.loc[ tree_frame["Tree"] == tree_idx ].set_index("ID") node_id = f"{tree_idx}-0" path = [] while nodes.loc[node_id, "Feature"] != "Leaf": node = nodes.loc[node_id] feature = node["Feature"] threshold = float(node["Split"]) went_yes = float(row[feature]) < threshold next_id = node["Yes"] if went_yes else node["No"] path.append( { "node": node_id, "rule": f"{feature} < {threshold:.4f}", "value": float(row[feature]), "branch": "yes" if went_yes else "no", "next": next_id, } ) node_id = next_id path.append( { "node": node_id, "rule": "Leaf", "value": np.nan, "branch": "", "next": "", "leaf_value": float(nodes.loc[node_id, "Gain"]), } ) return pd.DataFrame(path)
One path explains one tree, not the ensemble. Use it to verify suspicious thresholds, missing-value routing, and features that appear unexpectedly early.
5. Use SHAP for global and local explanations
The current SHAP plotting API consumes Explanation objects. A beeswarm summarizes feature effects across rows; a waterfall shows how one prediction moves from the expected model output to the final prediction.
explainer = shap.Explainer(model)
explain_sample = X_valid.sample(1_000, random_state=7)
shap_values = explainer(explain_sample)
shap.plots.beeswarm(
shap_values,
max_display=10,
show=False,
)
plt.tight_layout()
plt.show()
In the beeswarm, horizontal position shows whether a feature pushed a prediction up or down; color shows the observed feature value. Look for mixed colors on both sides, which often indicates interactions or non-monotonic behavior.
Now explain one difficult remote delivery:
pred = model.predict(X_valid) residuals = pd.DataFrame( { "actual": y_valid, "predicted": pred, "residual": y_valid - pred, "delivery_zone": metadata_valid["delivery_zone"].to_numpy(), }, index=X_valid.index, ) worst_remote_idx = ( residuals.query("delivery_zone == 'remote'")["residual"] .abs() .idxmax() ) one_explanation = explainer(X_valid.loc[[worst_remote_idx]]) shap.plots.waterfall( one_explanation[0], max_display=10, show=False, ) plt.tight_layout() plt.show() print(trace_tree(tree_frame, X_valid.loc[worst_remote_idx], tree_idx=0))
The waterfall explains the model, not reality. A missing operational variable can produce a perfectly faithful explanation of a systematically wrong prediction.
6. Separate average behavior from individual behavior
Partial dependence shows the average prediction as one feature changes. ICE draws one curve per sampled row. Scikit-learn supports both through PartialDependenceDisplay.from_estimator() with kind="both".
PartialDependenceDisplay.from_estimator(
model,
X_valid,
features=["distance_km", "traffic_index"],
kind="both",
subsample=200,
centered=True,
random_state=7,
n_cols=2,
)
plt.tight_layout()
plt.show()
A smooth average line can hide sharply different ICE curves. That spread is a prompt to search for interactions or segments. Also remember that dependence plots are not causal, and correlated features can create unrealistic feature combinations.
7. Find the slice carrying the error
Aggregate RMSE is the start of debugging, not the end. Compute bias, absolute error, and each segment’s share of squared error.
overall_rmse = mean_squared_error(
residuals["actual"],
residuals["predicted"],
) ** 0.5
overall_mae = mean_absolute_error(
residuals["actual"],
residuals["predicted"],
)
segment_report = residuals.groupby("delivery_zone").agg(
rows=("residual", "size"),
rmse=("residual", lambda s: np.sqrt(np.mean(s**2))),
mae=("residual", lambda s: np.mean(np.abs(s))),
mean_residual=("residual", "mean"),
squared_error=("residual", lambda s: np.sum(s**2)),
)
segment_report["share_squared_error"] = (
segment_report["squared_error"]
/ segment_report["squared_error"].sum()
)
print({"rmse": overall_rmse, "mae": overall_mae})
print(segment_report.sort_values("rmse", ascending=False).round(3))
Here is the cherry on the cake: with the fixed seed above, a typical current 3.x run shows remote deliveries at about 10% of validation rows yet roughly 62% of total squared error. Their mean residual is around +15.5 minutes, so the model systematically predicts them too low. The overall RMSE, about eight minutes, hides the operational failure.
The fix is not automatically “add the segment.” First verify that the field is available, stable, lawful to use, and present at prediction time. Then compare options: include it, engineer a reliable proxy, build a segment-specific model, change the loss or weighting, collect better data, or define separate acceptance thresholds.
A repeatable debugging order
Use the tools in this sequence:
-
Validate the split and metric.
-
Inspect train-versus-validation learning curves.
-
Compare weight, gain, and cover instead of trusting one ranking.
-
Read suspicious trees and exact paths.
-
Use SHAP globally, then explain concrete failures.
-
Use PDP and ICE to expose heterogeneous responses.
-
Slice residuals by business-relevant metadata and rank slices by error contribution.
Run this toolkit on your strongest boosted model today. Do not ship until you can name where it fails, quantify who absorbs the error, and show which evidence led you there.