Feature selection and SHAP: knowing which features matter and proving it

A model can achieve excellent predictive accuracy and still be fundamentally wrong about which features are useful.

Consider a delivery-time prediction system. It identifies a particular tracking field as its most important feature. Removing that field makes predictions substantially worse. A SHAP analysis confirms that the field dominates individual predictions.

Everything appears statistically convincing.

Then someone checks the data pipeline and discovers that the tracking field is populated only after a delivery is completed.

The model has effectively learned to predict the past.

This illustrates a crucial distinction in machine learning: a feature can be statistically important, essential to a trained model, and completely unsuitable for real-world prediction.

To understand feature importance properly, we need more than one ranking algorithm. We need methods for selecting features, explaining predictions, evaluating their contribution to generalization, and auditing the underlying data.

In this tutorial, we’ll build a reproducible workflow using:

• Mutual information for filter-based feature selection.

• Recursive feature elimination (RFE) for wrapper-based selection.

• L1 regularization for embedded selection.

• XGBoost for nonlinear predictive modeling.

• SHAP for local and global model explanations.

• Feature ablation to test whether supposedly important variables actually improve predictive performance.

• A controlled data-leakage experiment in which the highest-ranked feature turns out to be unusable.

Every example uses synthetic delivery-operation data, so no external dataset or API credentials are required.

1. What does it mean for a feature to matter?

Before implementing feature selection, distinguish four questions that are frequently confused.

Statistical relevance: Does a feature contain information about the target?

Predictive utility: Does including the feature improve performance on unseen observations?

Model attribution: How much does a feature contribute to the predictions of a particular trained model?

Operational validity: Will the feature be available, correctly measured, and stable when the model makes predictions?

These are related but different properties.

A feature can have strong statistical relevance without improving a model because another feature already contains the same information.

A model can give a feature a large SHAP attribution even when that information is unavailable at prediction time.

And a feature can be causally important in the real world while having little incremental predictive value because its effects are already captured through other variables.

This is why feature selection and model explainability complement each other rather than replace each other.

Three approaches to feature selection

SHAP addresses another question: how a fitted prediction function distributes its output among input features under a particular explanation framework.

It is not, by itself, a feature-selection algorithm or a test of causality.

2. Set up a reproducible environment

We’ll use Python 3.11 or newer, with recent releases of NumPy, pandas, scikit-learn, XGBoost, SHAP, and Matplotlib.

As of October 2026, SHAP 0.53.0 is a current release, published on October 9, 2026. The modern explanation interface returns structured explanation objects and supports plotting directly from those objects. See the SHAP release notes and TreeExplainer reference.

Create a virtual environment:

python -m venv .venv
source .venv/bin/activate
python -m pip install --upgrade pip

Install the dependencies:

python -m pip install \
  "numpy>=2.0,<3" \
  "pandas>=2.2,<3" \
  "scikit-learn>=1.8,<2" \
  "xgboost>=3.1,<4" \
  "shap>=0.50,<1" \
  "matplotlib>=3.9,<4"

Verify the installed versions:

python -m pip show scikit-learn xgboost shap

The examples below are intended to run sequentially in the same Python session or notebook. Alternatively, concatenate the Python blocks into a single script in their presented order.

For production deployment, lock and test the exact dependency versions rather than relying indefinitely on version ranges.

3. Build a realistic synthetic dataset

Imagine an e-commerce business that wants to predict the total delivery duration of an order when a courier is assigned.

At that moment, the system knows information such as:

• Estimated delivery distance.

• Route-duration estimate.

• Traffic conditions.

• Rainfall.

• Warehouse workload.

• Driver experience.

• Whether the order belongs to a promotional campaign.

• Whether delivery occurs on a weekend.

The actual delivery duration becomes known only after the delivery finishes.

We’ll generate 4,000 synthetic delivery records with meaningful predictive relationships.

Import the libraries and construct the data:

import numpy as np
import pandas as pd

from sklearn.model_selection import train_test_split

rng = np.random.default_rng(2026)
n = 4000

distance = rng.uniform(1, 40, n)
traffic = rng.uniform(0, 1, n)
rain = rng.exponential(3.0, n)
backlog = rng.poisson(10, n)
experience = rng.uniform(0, 12, n)
promo = rng.binomial(1, 0.15, n)
weekend = rng.binomial(1, 2 / 7, n)

y = (
    8
    + 2.4 * distance
    + 16 * traffic
    + 0.9 * rain
    + 0.5 * backlog
    - 0.65 * experience
    + 6 * promo
    + 2.5 * weekend
    + 8 * (traffic > 0.75) * promo
    + rng.normal(0, 3, n)
)

X = pd.DataFrame({
    "distance_km": distance,
    "route_estimate_min": (
        2.4 * distance + rng.normal(0, 2, n)
    ),
    "traffic_index": traffic,
    "rain_mm": rain,
    "warehouse_backlog": backlog,
    "driver_experience_years": experience,
    "promo_day": promo,
    "weekend": weekend,
    "random_noise": rng.normal(size=n),
    "scan_to_close_min": (
        y + rng.normal(0, 0.7, n)
    ),
})

leak_column = "scan_to_close_min"

valid_columns = X.columns.drop(leak_column).tolist()

X_train, X_test, y_train, y_test = train_test_split(
    X,
    y,
    test_size=0.25,
    random_state=42,
)

X_train_valid = X_train[valid_columns]
X_test_valid = X_test[valid_columns]

print("Training rows:", len(X_train))
print("Test rows:", len(X_test))

print(
    "Distance/route correlation:",
    X_train_valid["distance_km"].corr(
        X_train_valid["route_estimate_min"]
    ),
)

There are three deliberate complications in this dataset.

Complication A: correlated predictors

The distance and route-estimate variables contain nearly identical information.

Their Pearson correlation is approximately 0.997.

Both can be useful predictors, but the model may distribute importance between them in different ways.

This will become important when interpreting SHAP values and L1 coefficients.

Complication B: nonlinear interactions

Promotional periods introduce an additional delay when traffic is particularly heavy.

Neither a simple linear coefficient nor a univariate feature score can fully represent that interaction.

A boosted tree model can learn such relationships.

Complication C: a leaked outcome proxy

The field named scan_to_close_min represents a finalized duration captured by a tracking system after delivery.

It is highly predictive because we deliberately created it from the outcome plus a small amount of measurement noise.

But it is not available at dispatch time.

We’ll initially allow it into one deliberately flawed model, then compare that model with a legitimate one.

The dataset is synthetic, but the failure mode mirrors a common production problem: a feature looks excellent in offline testing because it contains information from the future.

4. Filter selection with mutual information

A filter method scores features without fitting the final predictive model.

Mutual information measures statistical dependence between a feature and a target.

In plain language:

• Zero mutual information indicates independence at the population level.

• Higher mutual information indicates greater statistical dependence.

• Unlike ordinary correlation, mutual information can detect nonlinear relationships.

• Mutual information does not tell us whether the relationship is causal.

For regression tasks, scikit-learn provides the mutual_info_regression estimator.

We’ll use SelectKBest to retain five features.

from functools import partial

from sklearn.feature_selection import (
    SelectKBest,
    mutual_info_regression,
)

mi_selector = SelectKBest(
    score_func=partial(
        mutual_info_regression,
        random_state=42,
    ),
    k=5,
)

mi_selector.fit(X_train_valid, y_train)

mi_scores = pd.Series(
    mi_selector.scores_,
    index=valid_columns,
).sort_values(ascending=False)

mi_features = X_train_valid.columns[
    mi_selector.get_support()
].tolist()

print("Mutual information ranking:")
print(mi_scores.round(3))

print("Selected features:")
print(mi_features)

A representative result from the seeded example is:

distance_km                1.243
route_estimate_min         1.219
traffic_index              0.025
promo_day                  0.015
weekend                    0.012
driver_experience_years    0.011
random_noise               0.003
warehouse_backlog          0.000
rain_mm                    0.000

Notice something surprising.

The rainfall and backlog variables genuinely influence delivery duration in our data-generation process, yet their estimated marginal mutual information can be extremely small.

Why?

Because delivery distance dominates variation in the target, and each filter score examines a feature individually rather than conditioning on the other predictors.

A nonparametric mutual-information estimate also has finite-sample uncertainty. Small positive signals may be estimated as zero.

Consequently, a low estimated score is not a mathematical proof that a feature is useless.

Where mutual information works well

Mutual information is particularly useful when:

• There are thousands of candidate predictors.

• The target-feature relationship may be nonlinear.

• Training sophisticated models on every feature is expensive.

• We need an inexpensive preliminary ranking.

However, marginal mutual information does not directly penalize redundancy.

In our example, it ranks both distance and route duration highly, despite their near-duplicate information.

Selecting both may waste part of a limited feature budget.

For mixed categorical and numerical datasets, correctly identify discrete variables using the estimator’s discrete_features argument. Our continuous, dense example uses its default behavior; a discrete-variable-aware implementation would be preferable when discrete predictors are central to the problem.

Also remember that mutual information is measured in information units, not delivery minutes. Its numerical value cannot be interpreted as a change in the prediction.

5. Wrapper selection with recursive feature elimination

Recursive feature elimination takes a different approach.

Rather than scoring every variable independently, RFE repeatedly fits an estimator and removes the least important features.

The procedure is conceptually straightforward:

1. Fit an estimator using the current feature set.

2. Rank features using that estimator’s coefficients or importances.

3. Remove the lowest-ranked features.

4. Fit the estimator again.

5. Continue until the desired number of features remains.

We’ll use ridge regression as the estimator.

Because ridge coefficients depend on feature scale, we’ll standardize the predictors inside the estimator pipeline.

from sklearn.feature_selection import RFE
from sklearn.linear_model import Ridge
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

rfe_selector = RFE(
    estimator=make_pipeline(
        StandardScaler(),
        Ridge(alpha=10.0),
    ),
    n_features_to_select=5,
    importance_getter="named_steps.ridge.coef_",
    step=1,
)

rfe_selector.fit(X_train_valid, y_train)

rfe_features = X_train_valid.columns[
    rfe_selector.support_
].tolist()

rfe_ranks = pd.Series(
    rfe_selector.ranking_,
    index=valid_columns,
).sort_values()

print("Selected by RFE:")
print(rfe_features)

print("RFE rankings:")
print(rfe_ranks)

The importance_getter argument tells RFE where to find the fitted ridge coefficients inside the pipeline.

For the seeded example, RFE selects:

distance_km
route_estimate_min
traffic_index
rain_mm
promo_day

This differs from the mutual-information result.

Rainfall enters the subset, but weekend status drops out.

Why?

RFE ranks features in the context of its chosen estimator. Ridge regression considers the combined linear relationship of the selected predictors with the target.

That allows features with relatively weak individual signals to become valuable alongside other features.

However, this particular RFE implementation still inherits the assumptions of ridge regression.

It does not automatically discover nonlinear interactions or optimize the predictive performance of the XGBoost model we’ll use later.

Does RFE find the optimal subset?

Not necessarily.

RFE is a greedy, backward-selection procedure. Once a feature has been removed, ordinary RFE does not reconsider it later.

It also requires choosing how many features to retain.

In production, consider RFECV, which evaluates different feature counts through cross-validation.

Use a suitable cross-validation strategy:

• Time-aware splits for chronological forecasting.

• Group-aware splits when observations share customers, accounts, or devices.

• Stratification when appropriate for classification tasks.

Do not use the test set to decide the optimal number of features.

For highly correlated predictors, RFE rankings can also be unstable. A small change in the training sample may cause one redundant variable to replace another.

6. Embedded selection with L1 regularization

Embedded feature-selection methods perform selection while fitting the predictive estimator itself.

L1 regularization is a classic example.

The idea is to penalize the absolute magnitude of model coefficients.

For linear regression, the L1-regularized objective can be expressed in plain text as:

Objective = Mean squared prediction error + α × Sum of absolute coefficient magnitudes

The exact normalization depends on the estimator implementation.

As the regularization strength α increases, coefficients are encouraged toward zero.

Some coefficients become exactly zero, allowing us to remove the corresponding features.

This is useful because selection is integrated with model optimization rather than performed through a separate ranking step.

Implement L1 selection with cross-validation

We’ll use Lasso, standardize its inputs, and tune the regularization strength through GridSearchCV.

Putting the scaler inside the searched pipeline ensures that scaling is fitted separately within each cross-validation training fold.

from sklearn.linear_model import Lasso
from sklearn.model_selection import GridSearchCV

l1_search = GridSearchCV(
    estimator=make_pipeline(
        StandardScaler(),
        Lasso(
            max_iter=20000,
            random_state=42,
        ),
    ),
    param_grid={
        "lasso__alpha": np.geomspace(
            0.01,
            3.0,
            15,
        ),
    },
    cv=5,
    scoring="neg_mean_absolute_error",
    n_jobs=-1,
)

l1_search.fit(X_train_valid, y_train)

l1_model = l1_search.best_estimator_

l1_coeffs = pd.Series(
    l1_model.named_steps["lasso"].coef_,
    index=valid_columns,
)

l1_features = l1_coeffs.index[
    l1_coeffs.abs() > 1e-6
].tolist()

print("Best alpha:")
print(l1_search.best_params_)

print("L1 coefficients:")
print(l1_coeffs.round(3))

print("Selected features:")
print(l1_features)

On the seeded data, L1 commonly retains distance, traffic, rainfall, warehouse backlog, driver experience, promotional status, and weekend status.

It sets the coefficient of the highly redundant route-estimate feature to zero.

That is a key demonstration of regularization’s behavior.

The route-estimate variable is not necessarily bad. It simply adds little unique information once distance is already present in this fitted linear model.

The correlated-feature problem with L1

When two variables are strongly correlated, L1 often retains one and suppresses the other.

Which one survives may depend on:

• Small sample differences.

• Measurement noise.

• Scaling.

• Regularization strength.

• Correlation structure.

Therefore, a zero L1 coefficient does not establish that the underlying variable is irrelevant.

For a more stable treatment of groups of correlated features, Elastic Net can be attractive because it combines L1 and L2 regularization.

One additional warning: our data contain nonlinear effects. Lasso is still fitting a linear model, so its selected subset need not be optimal for a tree ensemble.

7. Compare feature subsets using the same boosted model

We now have three candidate feature sets.

But feature rankings alone do not establish predictive usefulness.

To evaluate the selections, we’ll fit the same XGBoost regression configuration to each subset.

This ensures differences are driven primarily by the available predictors rather than by using entirely different model families.

XGBoost’s scikit-learn interface provides XGBRegressor, and its histogram tree-building method is suitable for this experiment. The current interface is described in the XGBoost Python documentation.

from sklearn.metrics import mean_absolute_error
from xgboost import XGBRegressor


def make_booster():
    return XGBRegressor(
        objective="reg:squarederror",
        n_estimators=240,
        max_depth=3,
        learning_rate=0.06,
        subsample=0.9,
        colsample_bytree=1.0,
        reg_lambda=2.0,
        tree_method="hist",
        n_jobs=4,
        random_state=42,
    )


def fit_and_mae(features):
    model = make_booster()

    model.fit(
        X_train[features],
        y_train,
    )

    predictions = model.predict(
        X_test[features]
    )

    mae = mean_absolute_error(
        y_test,
        predictions,
    )

    return model, mae


feature_sets = {
    "Leaky model": X.columns.tolist(),
    "All valid features": valid_columns,
    "Mutual information": mi_features,
    "RFE": rfe_features,
    "L1": l1_features,
    "Drop distance": [
        c for c in valid_columns
        if c != "distance_km"
    ],
    "Drop route estimate": [
        c for c in valid_columns
        if c != "route_estimate_min"
    ],
    "Drop both": [
        c for c in valid_columns
        if c not in {
            "distance_km",
            "route_estimate_min",
        }
    ],
}

models = {}
results = {}

for label, columns in feature_sets.items():
    model, mae = fit_and_mae(columns)

    models[label] = model
    results[label] = mae

print("Holdout MAE in minutes:")

for label, mae in results.items():
    print(f"{label:25s} {mae:.3f}")

What the results tell us

Here are results from the seeded experiment, rounded to two decimal places (our run: scikit-learn 1.9.1, XGBoost 3.4.1, SHAP 0.54.1, NumPy 2.5.4, pandas 2.3.3). Exact results can vary slightly with dependency versions and numerical implementations.

The L1-selected subset performs comparably to using all legitimate features.

Mutual information and RFE produce smaller subsets, but predictive quality deteriorates.

This does not mean those algorithms are generally inferior.

It means that their selected subsets, under these particular settings, are less suitable for this particular boosted model.

In another dataset, the rankings might be reversed.

Why we use mean absolute error

Mean absolute error has a direct operational interpretation.

An MAE of 2.62 minutes means the average absolute difference between predicted and observed delivery duration is approximately 2.62 minutes.

That is easier to communicate to an operations team than a dimensionless score.

Other relevant regression metrics include root mean squared error for emphasizing large mistakes and percentile errors for understanding tail behavior.

In a real forecasting system, we would also examine performance during bad weather, peak traffic, and promotional campaigns.

Important evaluation distinction: This comparison is an educational diagnostic, not a procedure for repeatedly choosing a winner on the test set. In a real project, compare feature-selection strategies on validation folds, lock the final method, and evaluate it once on an untouched test set.

8. SHAP: explaining what a boosted model actually learned

Feature selection tells us which inputs a model needs to perform well.

SHAP helps explain how a fitted model uses those inputs.

SHAP builds on the idea of Shapley values from cooperative game theory.

For a particular prediction, the model’s output is decomposed into a baseline value and feature contributions.

In plain text:

Prediction = Baseline prediction + SHAP(feature 1) + SHAP(feature 2) + … + SHAP(feature n)

A positive SHAP value means that a feature pushes the modeled output above its baseline for that observation.

A negative value means that it pushes the prediction below the baseline.

For our regression task, SHAP values have the same units as the model output: minutes.

That makes individual explanations intuitive.

Suppose a delivery has:

• A baseline predicted duration of 55 minutes.

• A distance contribution of +18 minutes.

• A traffic contribution of +7 minutes.

• A driver-experience contribution of −4 minutes.

• Contributions from other features totaling +2 minutes.

The final prediction is 78 minutes.

These contributions describe the prediction function, not necessarily the causal effects of intervening on each feature.

Use the current SHAP explanation interface

For tree ensembles, shap.TreeExplainer implements efficient Tree SHAP algorithms.

We’ll explain both the deliberately leaky model and the legitimate model.

We explicitly provide a background dataset and select the interventional feature-perturbation mode. The background data come only from training observations.

import shap


def explain_model(label, columns):
    model = models[label]

    background = X_train[columns].sample(
        n=120,
        random_state=42,
    )

    observations = X_test[columns].iloc[:200]

    explainer = shap.TreeExplainer(
        model,
        data=background,
        feature_perturbation="interventional",
        model_output="raw",
    )

    explanation = explainer(observations)

    reconstructed = (
        explanation.base_values
        + explanation.values.sum(axis=1)
    )

    predictions = model.predict(observations)

    assert np.allclose(
        reconstructed,
        predictions,
        atol=1e-3,
        rtol=1e-3,
    )

    importance = pd.Series(
        np.abs(explanation.values).mean(axis=0),
        index=columns,
    ).sort_values(ascending=False)

    return explanation, importance


leaky_explanation, leaky_importance = explain_model(
    "Leaky model",
    feature_sets["Leaky model"],
)

clean_explanation, clean_importance = explain_model(
    "All valid features",
    valid_columns,
)

print("Leaky model importance:")
print(leaky_importance.round(3))

print("Clean model importance:")
print(clean_importance.round(3))

The assertion checks an essential SHAP property: the baseline plus all feature contributions should reconstruct the model’s output, within numerical tolerance.

Notice that we’re using the Explanation interface rather than older plotting patterns based on manually supplied SHAP arrays.

The modern interface retains feature values, names, baselines, and attribution values together.

Global importance: mean absolute SHAP

A common global importance score is:

Global importance(feature j) = Mean over observations of |SHAP value for feature j|

This tells us how much a feature changes model predictions on average, regardless of direction.

It does not measure prediction error improvement or causal effect.

A feature can have large positive and negative contributions that cancel in their signed average. Taking absolute values avoids that cancellation.

However, this aggregation also removes directional information.

For example, both traffic congestion and driver experience may receive large absolute contributions despite typically pushing predictions in opposite directions.

9. Visualize global and local explanations

SHAP offers three especially useful visualizations for tabular regression.

Global bar plot

The global bar plot ranks features by their average absolute attribution.

import matplotlib.pyplot as plt

shap.plots.bar(
    clean_explanation,
    max_display=10,
    show=False,
)

plt.savefig(
    "shap_global_importance.png",
    dpi=150,
    bbox_inches="tight",
)

plt.close()

This plot answers:

Which features have the largest average contributions to predictions over the observations being explained?

For our legitimate delivery-duration model, distance is dominant, followed by other operational variables such as traffic.

Beeswarm plot

A beeswarm plot provides much more information than a single importance ranking.

Each point represents an observation-feature attribution.

The horizontal position shows the SHAP contribution, while the color generally represents the underlying feature value.

shap.plots.beeswarm(
    clean_explanation,
    max_display=10,
    show=False,
)

plt.savefig(
    "shap_beeswarm.png",
    dpi=150,
    bbox_inches="tight",
)

plt.close()

A beeswarm can reveal patterns such as:

• Longer distances generally increasing predicted duration.

• Heavy traffic pushing predictions upward.

• Greater driver experience often reducing predicted duration.

• Nonlinear relationships where the effect changes abruptly after a threshold.

It also reveals whether a feature has a consistent effect or highly variable contributions.

Two features can have identical global importance but completely different attribution distributions.

One might influence almost every prediction moderately, while another matters only during rare operational events.

Waterfall plot

The waterfall visualization explains a single observation.

shap.plots.waterfall(
    clean_explanation[0],
    max_display=10,
    show=False,
)

plt.savefig(
    "shap_single_delivery.png",
    dpi=150,
    bbox_inches="tight",
)

plt.close()

This provides a feature-by-feature breakdown of how the model moved from its baseline to the final predicted duration.

It is particularly useful when investigating operational outliers.

For example, a business analyst might ask why a delivery expected to take 30 minutes received a 75-minute estimate.

A waterfall plot can identify the features responsible for the model’s unusually high prediction.

But an explanation should always be interpreted alongside prediction accuracy.

An internally consistent explanation of an inaccurate prediction is still an explanation of an inaccurate prediction.

10. The cherry on the cake: when the best feature was a leaked answer

Now return to our deliberate data-leakage experiment.

The flawed model includes scan_to_close_min, a duration recorded after delivery completion.

The model trained with this feature achieves an MAE of approximately 0.65 minutes.

That is dramatically better than the legitimate model’s 2.62-minute MAE.

SHAP also ranks the tracking field as overwhelmingly important.

On our seeded example, its mean absolute SHAP contribution is approximately 23 minutes, dwarfing the attributions of the remaining features.

At first glance, it looks like the strongest feature discovered in the project.

But the field is really a proxy for the answer.

The model is predicting delivery duration using another measurement of the completed delivery duration.

Why ordinary validation misses the problem

A random train/test split separates rows, but it does not enforce the chronology of information within each row.

If both the training and test datasets contain finalized tracking fields, the relationship between those fields and the target remains available in both splits.

The model can exploit it without overfitting in the usual sense.

Even a chronological train/test split will not fix this problem if the future test rows still include information that was recorded after their prediction timestamps.

This is a key distinction:

Correctly separating observations is not the same as correctly separating information available before and after a prediction.

How the leakage is exposed

The decisive question is not whether the feature has high SHAP importance.

It is:

At the exact moment this model must predict, does the system actually know this value?

For the tracking field, the answer is no.

A proper point-in-time feature contract would specify:

• The prediction event: courier assignment.

• The prediction timestamp: assignment time.

• The feature-generation timestamp: when the value first becomes available.

• The allowed feature cutoff: no later than assignment.

• The production source: the data system that will provide the feature during inference.

The tracking field fails this contract.

When we remove it, the model loses its unrealistic predictive advantage but becomes suitable for deployment.

That is a better engineering outcome than retaining an impressive but unattainable test score.

What SHAP actually proved here

SHAP correctly identified that the fitted model relied heavily on the tracking field.

The explanation was not defective.

The assumption that a highly attributed feature must be operationally valid was defective.

This illustrates the boundary between model interpretability and data validity.

SHAP can reveal suspicious dependence, but an operational audit is needed to establish whether that dependence constitutes leakage.

11. The correlated-feature attribution trap

Let’s investigate a subtler problem involving distance and route estimates.

These variables have a correlation of approximately 0.997.

They represent closely related information.

When training with all valid features, the boosted model may use both.

SHAP then distributes contributions between them according to the fitted trees, background distribution, and explanation assumptions.

Meanwhile, L1 may eliminate the route-estimate variable entirely.

Who is right?

Both methods can be internally correct because they answer different questions.

L1 asks whether a coefficient is useful in a penalized linear prediction function.

SHAP explains how a particular fitted boosted model uses its inputs.

Neither necessarily measures the unique information supplied by an individual feature independently of correlated substitutes.

Test individual and grouped feature removal

Our earlier evaluation included three ablation experiments:

• Remove distance alone.

• Remove route estimate alone.

• Remove both simultaneously.

In the seeded run, the results are approximately:

Removing the route estimate barely affects performance.

Removing distance causes a modest deterioration.

Removing both destroys most of the model’s predictive capability.

This is not contradictory.

When one correlated feature disappears, the other still supplies much of the missing information.

When both disappear, the model loses the entire distance-related signal.

This is why correlated variables can appear individually expendable even when the underlying information is essential.

Group-level importance

For strongly correlated predictors, consider treating groups of variables as a unit.

For example:

• Geographic distance and route duration.

• Current price and discounted price.

• Total transaction amount and item subtotal.

• Temperature readings from adjacent sensors.

A group ablation removes the entire feature family and retrains the model.

Its performance difference answers a practical question:

How much predictive utility is lost when this source of information is unavailable?

That question is different from asking how much attribution an existing model assigns to the corresponding features.

An alternative is grouped permutation importance, where related variables are shuffled together. However, naive permutation can create unrealistic combinations and must be interpreted carefully.

Why SHAP’s dependence assumption matters

Tree SHAP needs a method for handling features considered absent from a partial coalition.

Two supported approaches are:

Interventional: Uses background data to replace hidden features, breaking dependencies between observed and background-substituted features under the chosen masking procedure.

Tree-path-dependent: Uses training-path statistics stored in the tree ensemble rather than requiring a separate background dataset.

These approaches can assign different contributions when features are strongly dependent.

In current SHAP versions, the auto setting chooses a method depending on whether background data are supplied.

The word interventional does not imply that the resulting explanations are automatically causal effects of real-world interventions.

For correlated measurements, arbitrary feature combinations can also be physically or operationally implausible.

A model may then be evaluated on combinations it almost never encounters in real life.

For deeper discussions of dependence assumptions and model support, consult the SHAP TreeExplainer documentation.

12. Go beyond explanation with feature ablation

Ablation testing supplies evidence that an explanation ranking alone cannot.

The basic procedure is:

1. Train and evaluate the reference model.

2. Remove a feature or group of features.

3. Retrain the same model family.

4. Evaluate under the same validation protocol.

5. Compare the predictive errors.

For a regression task:

Ablation cost = MAE without feature − MAE with feature

A positive ablation cost means performance got worse when the feature was removed.

A negative cost means removing the feature improved performance for that comparison.

This is a measure of the incremental predictive usefulness of the feature under the chosen training procedure, not a causal effect.

Let’s calculate the ablation costs from the models already fitted:

baseline_mae = results["All valid features"]

ablation_rows = []

for label in [
    "Drop distance",
    "Drop route estimate",
    "Drop both",
]:
    ablation_rows.append({
        "ablation": label,
        "mae": results[label],
        "additional_error": (
            results[label] - baseline_mae
        ),
    })

ablation_table = pd.DataFrame(ablation_rows)

print(
    ablation_table.sort_values(
        "additional_error",
        ascending=False,
    ).round(3)
)

For a rigorous experiment, the same procedure should be repeated across appropriate validation folds or multiple training seeds.

A single comparison can be sensitive to sampling noise and model optimization.

The difference between ablation and permutation

These two methods are often grouped together but measure different behaviors.

Retrained ablation: Removes features and allows the remaining model to adapt. It estimates the predictive information available from the retained feature set under a specified training procedure.

Permutation importance: Keeps the fitted model fixed but disrupts a feature’s association with the target by shuffling its observed values.

Permutation importance is useful for auditing an existing model, but correlated substitutes can mask the importance of individual variables.

Marginal permutation can also create unrealistic data combinations.

Neither method establishes a causal effect without additional assumptions and a suitable causal study design.

13. Feature selection must happen inside validation

One of the most common feature-selection mistakes is running a selector on the entire dataset before cross-validation.

That lets information from validation observations influence which predictors are used.

Even though the selector is not the final predictive model, supervised selection has still seen the validation targets.

This is a form of evaluation leakage.

The correct approach is to fit selection on each training fold.

A scikit-learn pipeline makes that straightforward.

For example, we can evaluate mutual-information selection and ridge regression together:

from sklearn.model_selection import cross_val_score

selection_pipeline = make_pipeline(
    SelectKBest(
        score_func=partial(
            mutual_info_regression,
            random_state=42,
        ),
        k=5,
    ),
    StandardScaler(),
    Ridge(alpha=10.0),
)

scores = cross_val_score(
    selection_pipeline,
    X_train_valid,
    y_train,
    scoring="neg_mean_absolute_error",
    cv=5,
    n_jobs=-1,
)

print("Mean validation MAE:", -scores.mean())
print("MAE standard deviation:", scores.std())

On every fold, the selector learns feature scores using only that fold’s training data.

The scaler and regression estimator are also fitted using training-fold information only.

This keeps validation-fold targets out of the selection process.

If feature count, regularization strength, and model hyperparameters must all be tuned, put them inside a search procedure and use nested cross-validation when an unbiased estimate of the entire tuning process is required.

For chronological forecasting, replace ordinary shuffled or sequential k-fold designs with a time-aware evaluation strategy that respects the actual prediction horizon.

A time-aware split is still insufficient if features violate point-in-time availability.

14. Selection stability: would the same features win again?

A feature list derived from one dataset split can create a false impression of certainty.

Suppose an L1 model chooses distance and rejects route estimate.

On a different sample, it might reverse that decision.

If the predictive performance stays nearly identical, the distinction between the two selected features may have little practical meaning.

We can study this by repeatedly fitting the selector to bootstrap samples.

The following example estimates how frequently each variable receives a nonzero L1 coefficient using the regularization strength selected earlier.

from sklearn.base import clone

bootstrap_rng = np.random.default_rng(123)
n_bootstraps = 30

selection_counts = pd.Series(
    0,
    index=valid_columns,
    dtype=float,
)

base_l1_model = l1_search.best_estimator_

for _ in range(n_bootstraps):
    indices = bootstrap_rng.integers(
        0,
        len(X_train_valid),
        size=len(X_train_valid),
    )

    X_boot = X_train_valid.iloc[indices]
    y_boot = y_train[indices]

    estimator = clone(base_l1_model)
    estimator.fit(X_boot, y_boot)

    coefficients = estimator.named_steps[
        "lasso"
    ].coef_

    selection_counts += (
        np.abs(coefficients) > 1e-6
    ).astype(float)

selection_frequency = (
    selection_counts / n_bootstraps
)

print(
    selection_frequency.sort_values(
        ascending=False
    ).round(2)
)

A frequency of 1.0 means the feature was selected in every bootstrap sample.

A frequency of 0.2 means it appeared in only 20% of those samples.

These numbers describe stability under this particular resampling procedure and fixed regularization setting. They are not probabilities that a feature is genuinely causal or scientifically meaningful.

For more rigorous stability studies, repeat hyperparameter selection inside each resample and preserve the temporal or grouped sampling structure of the original data.

Interpreting unstable feature selections

Instability can indicate:

• Highly correlated predictors.

• Insufficient training data.

• Weak signal relative to noise.

• Sensitivity to regularization settings.

• Changes in the data-generating environment.

But instability is not necessarily evidence of poor predictive performance.

Models can make equally accurate predictions using alternative representations of the same underlying signal.

The operational question is whether the particular feature identities matter for collection cost, resilience, governance, or deployment.

15. Explainability is not the same as causal inference

This distinction deserves explicit attention because SHAP plots are visually persuasive.

Suppose our model assigns a positive SHAP value to promotional campaigns.

We should not conclude that canceling a campaign would reduce delivery duration by exactly that number of minutes.

The SHAP value describes a contribution to a model prediction under a particular attribution framework.

A real intervention could change other variables simultaneously.

For example, canceling a promotion might affect warehouse backlog, courier availability, traffic patterns, and order distribution.

Similarly, a model that associates driver experience with shorter deliveries does not automatically quantify the effect of giving an individual driver another year of experience.

Causal claims require an appropriate causal framework, such as randomized experimentation or carefully justified observational identification assumptions.

SHAP is useful for understanding prediction functions. It is not a replacement for those methods.

A disciplined report should distinguish:

• Predictive association: A variable helps forecast the outcome.

• Model reliance: The fitted model uses a variable to generate predictions.

• Causal effect: Intervening on the variable changes the outcome under specified assumptions.

Confusing these claims can lead to poor business decisions even when model explanations are calculated perfectly.

16. Build a production-ready feature audit

A reliable feature-importance report should include much more than a colorful bar chart.

For each important feature, document the following dimensions.

For example, a real-time traffic feed may rank highly but be expensive or occasionally unavailable.

A slightly less accurate model using reliable historical traffic estimates might be preferable in production.

The best feature set is not necessarily the one that achieves the absolute lowest offline error.

It is the one that provides an appropriate balance of predictive quality, reliability, cost, latency, and operational validity.

Monitoring after deployment

Feature importance is not necessarily stable over time.

Changes in customer behavior, logistics networks, promotions, or data collection can alter how a model uses its inputs.

A practical monitoring strategy should track:

• Input distribution drift.

• Prediction accuracy when labels become available.

• Missing-value rates.

• Feature freshness and pipeline delays.

• SHAP attribution distributions on representative monitored samples.

• Performance across relevant operational segments.

An attribution shift does not automatically prove performance degradation.

It can indicate that the model is encountering different inputs, that correlated-feature usage has changed, or that data collection is malfunctioning.

Such changes should trigger investigation rather than automatic causal conclusions.

17. Common failure modes and how to avoid them

Mistake 1: Treating high SHAP importance as proof of causality

A feature can be highly attributed because it is predictive, correlated with another variable, or leaked from the future.

Better approach: Combine explanations with data-lineage checks, ablation experiments, and causal methods when causal questions matter.

Mistake 2: Selecting features before splitting data

Supervised feature selection can leak target information into evaluation.

Better approach: Fit selectors within training folds using pipelines.

Mistake 3: Ignoring feature correlations

Correlated predictors can divide attribution, suppress each other’s coefficients, and distort individual importance rankings.

Better approach: Inspect correlation structure, evaluate related features as groups, and test selection stability.

Mistake 4: Comparing importance values across different models as though they were universal

Ridge coefficients, mutual-information scores, XGBoost gain, and SHAP magnitudes do not measure the same property.

Better approach: Interpret every importance measure according to the model, scale, population, and definition used to calculate it.

Mistake 5: Assuming a smaller feature set is automatically better

Removing variables can improve simplicity but reduce predictive performance, robustness, or coverage of rare conditions.

Better approach: Compare subsets using the same model family, appropriate validation, and deployment constraints.

Mistake 6: Computing SHAP on an unrepresentative sample

Global SHAP importance depends on the observations being explained.

Explaining 200 easy deliveries might produce a different ranking from explaining 200 severe-delay cases.

The background dataset also influences attributions under background-dependent explanation modes.

Better approach: Define the population being explained, select representative samples, and document the background-data policy.

Mistake 7: Using SHAP to explain the wrong model output

SHAP contributions depend on the output scale.

For a standard regression model predicting minutes, raw SHAP values are expressed in minutes.

For binary classifiers, raw outputs from models such as XGBoost commonly represent margins or log-odds rather than probabilities.

Better approach: State the explained output explicitly and verify its additive reconstruction. Do not interpret log-odds contributions as percentage-point changes.

18. A decision framework for real projects

After working through these techniques, the choice between them becomes clearer.

These are starting points, not universal rules.

For large feature spaces, a particularly effective strategy is staged selection:

1. Remove invalid, unavailable, duplicated, and obviously unusable columns.

2. Apply an inexpensive filter to reduce the candidate space if necessary.

3. Evaluate promising subsets with wrapper or embedded methods.

4. Tune the final prediction model using validation data.

5. Compute SHAP explanations for the trained candidate.

6. Audit the most influential features for leakage, redundancy, and robustness.

7. Confirm the final model on a genuinely untouched test set.

This separates computational efficiency from predictive evaluation and operational verification.

19. What counts as convincing evidence that a feature matters?

There is no single statistic that proves every aspect of feature importance.

A defensible conclusion should be built from independent pieces of evidence.

Consider a feature that:

• Is available before prediction.

• Has a stable, documented data source.

• Improves validation performance when included.

• Retains its usefulness across time periods.

• Has meaningful SHAP contributions in the final model.

• Provides incremental information beyond correlated alternatives.

• Does not cause unacceptable operational or governance risks.

We now have a much stronger argument for keeping that feature.

Still, such evidence establishes a practical predictive case, not automatically a causal one.

This is the essential difference between producing an importance ranking and developing a trustworthy explanation.

A ranking tells us where to look.

A validation experiment tests predictive usefulness.

A data audit establishes whether the information belongs in the model.

Together, these methods create a meaningful basis for decisions.

20. Further reading and current documentation

For production implementations, use the actively maintained project documentation rather than copying old examples from notebooks that may rely on removed or superseded APIs.

• scikit-learn feature selection guide — filter, wrapper, and embedded methods.

• scikit-learn RFE API — recursive selection and importance extraction.

• scikit-learn mutual information API — continuous-target dependence estimation.

• SHAP API reference — explanation objects, explainers, and plots.

• SHAP TreeExplainer — tree-model explanation and feature-dependence assumptions.

• XGBoost Python API — current scikit-learn estimator interface.

Conclusion: move from importance rankings to evidence

The central lesson is that knowing which features a model uses is not the same as proving which features it needs or should be allowed to use.

Mutual information efficiently identifies statistical dependence.

RFE searches for useful subsets through repeated model fitting.

L1 regularization creates sparse models by penalizing coefficients.

SHAP decomposes fitted model predictions into feature attributions.

Ablation experiments test incremental predictive value.

Point-in-time data audits detect features that appear powerful only because they contain information the model would not have in production.

Our synthetic delivery experiment demonstrated all of these distinctions. The strongest feature in the flawed model was also the one that had to be removed. Correlated features produced conflicting-looking importance signals, yet grouped ablation clarified their combined predictive value.

The strongest workflow does not ask one algorithm to settle every question. It combines complementary methods and checks whether their conclusions remain valid under real operational constraints.

Your next step: Run the complete experiment, inspect the SHAP plots, then modify the synthetic dataset. Increase the correlation between distance and route estimates, strengthen the nonlinear traffic interaction, and vary the noise in the leaked tracking field. Repeat the selection and ablation analyses after each change.

Watch how the feature rankings shift, how predictive accuracy responds, and how SHAP divides credit among correlated variables.

Then apply the same process to one of your own non-sensitive datasets, beginning with a strict feature-availability audit and ending with a held-out performance comparison.

The objective is not simply to find a top-ten list of predictors. It is to build a model whose features are useful, explainable, available when needed, and supported by credible evidence.