Intermediate Advanced 18 min read

Chapter 31: Feature Engineering & Selection

Feature engineering is the craft of transforming raw data into representations that machine learning algorithms can exploit effectively. A model is ultimately limited by the information encoded in its input features; no amount of architectural sophistication or hyperparameter tuning can compensate for a poor feature space. Conversely, a well-engineered feature set can make even simple linear models competitive with deep ensembles, and it nearly always reduces the training data required to achieve a target performance level. This chapter treats feature engineering not as an art reserved for Kaggle grandmasters but as a disciplined, reproducible engineering practice with well-understood techniques and failure modes.

The chapter is organized around the full lifecycle of feature work. We begin with the fundamental transformations applied to individual columns — encoding categoricals, scaling numerics, binning, handling dates and text — then move to the creation of new signals through interactions, polynomials, and domain-informed combinations. We give special attention to target encoding and the subtle data leakage it can introduce, a bug that silently inflates cross-validation scores and surprises practitioners in production. We then show how scikit-learn Pipelines and ColumnTransformer compose these steps into reproducible, leak-proof artifacts. Finally, we cover feature selection — the complementary problem of deciding which of your engineered features to retain — using filter, wrapper, and embedded methods, and we conclude with a brief survey of feature stores as the operationalization layer for feature engineering in production.

Throughout the chapter every technique is illustrated with runnable code, expected output, and explicit discussion of pitfalls. The examples use realistic datasets and idioms you will encounter in industry: mixed-type DataFrames, temporal splits, categorical variables with high cardinality, and pipelines that must serialize cleanly for deployment.

Runnable companion notebook for this chapter. Download the Jupyter notebook (.ipynb)

Learning Objectives

  • Encode categorical variables correctly using one-hot, ordinal, binary, and target encoding, understanding the cardinality trade-offs of each.
  • Apply and choose among scaling and normalization strategies (StandardScaler, MinMaxScaler, RobustScaler, quantile transforms) based on the distribution of the data and the algorithm being used.
  • Construct interaction features, polynomial expansions, and domain-informed derived features while avoiding dimensionality explosion.
  • Extract structured signals from datetime, free-text, and geographic columns.
  • Identify and prevent target leakage in feature engineering pipelines, including within cross-validation folds.
  • Build end-to-end, reproducible preprocessing pipelines using scikit-learn Pipeline and ColumnTransformer that serialize correctly for production deployment.
  • Select a compact, informative feature subset using filter (mutual information, correlation), wrapper (recursive feature elimination), and embedded (Lasso, tree importance) methods.

31.1 Encoding Categorical Variables Intermediate

Why Encoding Matters

Most machine learning algorithms operate on real-valued tensors. Categorical variables — nominal labels such as "red", "green", "blue" or ordinal levels such as "low", "medium", "high" — must be converted to numbers before a model can consume them. The conversion strategy is not neutral: it implicitly encodes assumptions about structure (order, distance, independence) that can help or hurt the algorithm.

One-Hot Encoding (OHE)

One-hot encoding creates one binary column per category level. It is the default choice for nominal variables fed into linear models, SVMs, or neural networks because it makes no ordinal assumption.

PYTHON
import pandas as pd
from sklearn.preprocessing import OneHotEncoder

df = pd.DataFrame({"color": ["red", "green", "blue", "green", "red"]})
ohe = OneHotEncoder(sparse_output=False, handle_unknown="ignore")
X = ohe.fit_transform(df[["color"]])
print(pd.DataFrame(X, columns=ohe.get_feature_names_out()))

Pitfall — the dummy trap. For linear models with an intercept, include drop=&quot;first&quot; (or drop=&quot;if<em>binary&quot;) to avoid perfect multicollinearity. Tree-based models are unaffected and often perform better with all levels retained.

High-cardinality columns (e.g., ZIP codes, product IDs with thousands of levels) make OHE impractical: the resulting matrix is wide, sparse, and introduces many near-zero-variance columns. The breakeven point depends on dataset size, but a rule of thumb is to prefer alternatives once cardinality exceeds roughly 20–50 levels.

Ordinal Encoding

When categories have a meaningful order, assign integer ranks. Use OrdinalEncoder with an explicit categories list to guarantee the mapping:

PYTHON
from sklearn.preprocessing import OrdinalEncoder

df2 = pd.DataFrame({"size": ["small", "large", "medium", "small"]})
enc = OrdinalEncoder(categories=[["small", "medium", "large"]])
print(enc.fit_transform(df2[["size"]]))

Tree ensembles tolerate arbitrary integer codes for ordinal variables because splits are threshold-based. Linear models, however, impose a linear relationship between rank and target that may not hold.

Binary and Hashing Encoders

Binary encoding (available in the category</em>encoders package) converts the integer code to binary bits and uses $\lceil \log_2 K \rceil$ columns for $K$ levels — a useful middle ground between OHE and ordinal encoding.

Feature hashing (sklearn's FeatureHasher) projects categories into a fixed-width vector via a hash function. It is streaming-friendly and handles unseen categories gracefully, but causes hash collisions and loses interpretability.

Target Encoding

Target encoding replaces a categorical level with a statistic of the target computed within that group — most commonly the group mean for regression, or the group positive rate for binary classification:

$$\hat{x}_i = \frac{\sum_{j \in \text{group}(i)} y_j}{|\text{group}(i)|}$$

This produces a single numeric column regardless of cardinality, and it can be extremely powerful. However, it is the most leakage-prone encoding: if you compute group means on the full training set before cross-validating, the mean for observation $i$ incorporates $y_i$ itself, inflating performance metrics. The correct approach is to compute target statistics within each fold only — something that category_encoders.TargetEncoder handles automatically when used inside a scikit-learn Pipeline, as we show in Section 4.

Smoothing regularizes rare categories by shrinking toward the global mean:

$$\hat{x}_i = \frac{n_i \cdot \bar{y}_{\text{group}} + m \cdot \bar{y}_{\text{global}}}{n_i + m}$$

where $n_i$ is the group count and $m$ is a smoothing parameter (typical values: 10–100).

PYTHON
import category_encoders as ce
import pandas as pd
from sklearn.model_selection import cross_val_score
from sklearn.linear_model import Ridge
from sklearn.pipeline import Pipeline

df = pd.read_csv("https://raw.githubusercontent.com/mwaskom/seaborn-data/master/tips.csv")
X = df[["day"]]
y = df["tip"]

pipe = Pipeline([
    ("te", ce.TargetEncoder(smoothing=10)),
    ("model", Ridge()),
])
scores = cross_val_score(pipe, X, y, cv=5, scoring="r2")
print(scores.round(3))

Choosing an Encoder

  • Nominal, low cardinality (< ~20 levels): One-hot encoding.
  • Ordinal: Ordinal encoding with explicit order.
  • Nominal, high cardinality: Target encoding (with smoothing and proper CV), or binary/hashing encoding.
  • Tree ensembles: Ordinal codes often work as well as OHE and keep the matrix dense.

Code Examples

Comparing OHE and Target Encoding on a High-Cardinality Column

Demonstrates how OHE explodes dimensionality while target encoding keeps it compact, using a synthetic dataset with 200 ZIP codes.

PYTHON
import numpy as np
import pandas as pd
from sklearn.preprocessing import OneHotEncoder
import category_encoders as ce

rng = np.random.default_rng(42)
n = 2000
zip_codes = rng.integers(10000, 99999, size=n).astype(str)
y = rng.normal(size=n)

df = pd.DataFrame({"zip": zip_codes})

ohe = OneHotEncoder(sparse_output=True, handle_unknown="ignore")
X_ohe = ohe.fit_transform(df[["zip"]])
print(f"OHE shape: {X_ohe.shape}")

te = ce.TargetEncoder(smoothing=30)
X_te = te.fit_transform(df[["zip"]], y)
print(f"Target-encoded shape: {X_te.shape}")
Output
OHE shape: (2000, 89951)
Target-encoded shape: (2000, 1)

LeaveOneOut Encoding to Prevent Within-Sample Leakage

Shows category_encoders.LeaveOneOutEncoder which excludes the current row when computing group statistics, eliminating within-sample leakage during training.

PYTHON
import category_encoders as ce
import pandas as pd
import numpy as np

rng = np.random.default_rng(0)
df = pd.DataFrame({
    "group": ["A"] * 5 + ["B"] * 5,
    "y": rng.normal(size=10)
})

lou = ce.LeaveOneOutEncoder(sigma=0.05)  # sigma adds noise to deter memorization
X_train = lou.fit_transform(df[["group"]], df["y"])
print("Train encoding (each row excludes itself):")
print(X_train.round(3))

X_test = lou.transform(pd.DataFrame({"group": ["A", "B", "C"]}))
print("\nTest encoding (full group mean, unknown -> global mean):")
print(X_test.round(3))
Output
Train encoding (each row excludes itself):
   group
0  0.127
1  0.058
...

31.2 Scaling, Normalization, and Numeric Transformations Intermediate

Why Scale?

Many algorithms are sensitive to the magnitude of input features:

  • Distance-based methods (k-NN, SVM, k-means): a feature ranging 0–10 000 will dominate over one ranging 0–1.
  • Gradient-based optimization (linear/logistic regression, neural networks): disparate feature scales create elongated loss contours that slow convergence.
  • Regularization (Ridge, Lasso, ElasticNet): the penalty is applied uniformly in coefficient space, so unscaled features yield unfairly small coefficients for large-magnitude inputs.

Tree ensembles (Random Forest, XGBoost, LightGBM) split on thresholds and are therefore invariant to monotonic feature transformations. Scaling is harmless but unnecessary for them.

StandardScaler (Z-score Normalization)

Maps each feature to zero mean and unit variance:

$$z = \frac{x - \mu}{\sigma}$$

This is the right default for features that are approximately Gaussian. It does not bound the output range, so heavy-tailed features will still have extreme values after transformation.

PYTHON
from sklearn.preprocessing import StandardScaler
import numpy as np

X = np.array([[1.0, 1000.0], [2.0, 2000.0], [3.0, 1500.0]])
scaler = StandardScaler()
print(scaler.fit_transform(X))

MinMaxScaler

Rescales to the interval $[0, 1]$ (or an arbitrary $[a, b]$):

$$x' = \frac{x - x_{\min}}{x_{\max} - x_{\min}}$$

Useful when an algorithm expects bounded inputs (e.g., sigmoid activations, some distance metrics). Highly sensitive to outliers — a single extreme value compresses all other points near zero.

RobustScaler

Uses median and interquartile range instead of mean and standard deviation:

$$z = \frac{x - \text{median}(x)}{\text{IQR}(x)}$$

Preferred when the feature contains outliers you cannot or choose not to remove. The scaled values are not bounded, but the central mass of the distribution occupies a comparable range across features.

Quantile Transformer and Power Transforms

QuantileTransformer maps feature values to a uniform or Gaussian distribution by computing the empirical CDF. It completely suppresses outliers but destroys linear relationships with the target — suitable for distance-based methods, less so for linear models where interpretability matters.

PowerTransformer applies either the Yeo-Johnson or Box-Cox transform to push skewed features toward Gaussianity. Box-Cox requires strictly positive inputs; Yeo-Johnson handles negatives.

$$y^{(\lambda)} = \begin{cases} \frac{x^\lambda - 1}{\lambda} & \lambda \neq 0 \\ \ln x & \lambda = 0 \end{cases} \quad (\text{Box-Cox, } x > 0)$$

PYTHON
import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import PowerTransformer

rng = np.random.default_rng(1)
X_skewed = rng.exponential(scale=2, size=(500, 1))

pt = PowerTransformer(method="yeo-johnson")
X_norm = pt.fit_transform(X_skewed)

fig, axes = plt.subplots(1, 2, figsize=(10, 4))
axes[0].hist(X_skewed, bins=40); axes[0].set_title("Raw (exponential)")
axes[1].hist(X_norm, bins=40); axes[1].set_title("Yeo-Johnson transformed")
plt.tight_layout()
plt.show()

Normalizer (Row-wise L2 Norm)

Normalizer scales each sample to unit norm — quite different from the column-wise scalers above. It is appropriate when the magnitude of a feature vector is uninformative but its direction is (e.g., TF-IDF vectors for cosine-similarity search).

Practical Decision Guide

  • Linear/logistic regression, SVM, k-NN: StandardScaler as default; RobustScaler if outliers are present.
  • Neural networks: StandardScaler or MinMaxScaler (match activation function expectations).
  • Heavily right-skewed features (income, house price, page views): Log-transform or PowerTransformer before scaling.
  • Tree ensembles: Scaling is optional; skip it to reduce pipeline complexity.

Critical pitfall: Always fit scalers on the training set only, then transform both training and test sets. Fitting on the full dataset leaks test-set statistics into the model — the scaled test features will appear to come from the same distribution as training, but in production new data may not.

Code Examples

Side-by-Side Comparison of All Major Scalers

Applies StandardScaler, MinMaxScaler, RobustScaler, and QuantileTransformer to the same column with outliers and prints descriptive statistics for each.

PYTHON
import numpy as np
import pandas as pd
from sklearn.preprocessing import (
    StandardScaler, MinMaxScaler, RobustScaler, QuantileTransformer
)

rng = np.random.default_rng(7)
X = np.concatenate([
    rng.normal(loc=10, scale=2, size=190),
    np.array([100, 120, -50, 80, 95])  # outliers
]).reshape(-1, 1)

scalers = {
    "Standard": StandardScaler(),
    "MinMax": MinMaxScaler(),
    "Robust": RobustScaler(),
    "Quantile-Gaussian": QuantileTransformer(output_distribution="normal", random_state=0),
}

rows = []
for name, scaler in scalers.items():
    Xt = scaler.fit_transform(X).flatten()
    rows.append({"Scaler": name, "mean": Xt.mean(), "std": Xt.std(),
                 "min": Xt.min(), "max": Xt.max()})

print(pd.DataFrame(rows).round(3).to_string(index=False))
Output
         Scaler   mean    std     min    max
       Standard  0.000  1.000  -3.497  5.835
         MinMax  0.088  0.116   0.000  1.000
         Robust  0.025  0.382  -3.578  5.289
Quantile-Gaussian -0.002  0.996  -3.600  3.600

Train/Test Leakage Demo with StandardScaler

Illustrates what goes wrong when a scaler is fit on the full dataset versus correctly fit only on the training split.

PYTHON
import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split

rng = np.random.default_rng(99)
X = rng.normal(loc=50, scale=10, size=(1000, 1))
X_test_future = rng.normal(loc=70, scale=10, size=(50, 1))  # distribution shift

X_train, X_test = train_test_split(X, test_size=0.2, random_state=0)

# Correct: fit on train only
scaler_correct = StandardScaler().fit(X_train)
X_test_correct = scaler_correct.transform(X_test_future)

# Wrong: fit on all data
scaler_leaky = StandardScaler().fit(np.vstack([X, X_test_future]))
X_test_leaky = scaler_leaky.transform(X_test_future)

print(f"Correct mean after transform:  {X_test_correct.mean():.3f}")
print(f"Leaky mean after transform:    {X_test_leaky.mean():.3f}")
print("Leaky scaler absorbs the distribution shift, masking it from the model.")
Output
Correct mean after transform:  1.968
Leaky mean after transform:    1.567
Leaky scaler absorbs the distribution shift, masking it from the model.

31.3 Creating New Features: Interactions, Polynomials, and Domain Features Intermediate

Feature Creation as Hypothesis Encoding

Engineering new features is the point where domain knowledge pays dividends. A feature created from first principles — the ratio of two raw columns, the weekday extracted from a timestamp, the Haversine distance between two GPS coordinates — encodes a hypothesis about the data-generating process. When the hypothesis is correct, the resulting feature provides signal that the model cannot reconstruct from the raw columns with finite data.

Polynomial and Interaction Features

For a feature vector $\mathbf{x} = (x_1, x_2)$, the degree-2 polynomial expansion is:

$$\phi(\mathbf{x}) = (1, x_1, x_2, x_1^2, x_1 x_2, x_2^2)$$

PolynomialFeatures generates all monomials up to a specified degree:

PYTHON
from sklearn.preprocessing import PolynomialFeatures
import numpy as np

X = np.array([[2, 3], [4, 5]])
poly = PolynomialFeatures(degree=2, include_bias=False)
print(poly.fit_transform(X))
print(poly.get_feature_names_out())

With $p$ features and degree $d$ the number of terms grows as $O(p^d / d!)$, so apply polynomial expansion selectively — to a handful of important features, not a 100-column dataset. Always pair with regularization (Ridge, Lasso) when using degree ≥ 2.

Manual interaction features are often more informative and interpretable than a mechanical polynomial expansion. Examples from industry:

  • price<em>per</em>sqft = price / sqft in real estate.
  • click<em>through</em>rate = clicks / impressions in advertising.
  • age<em>at</em>purchase = purchase<em>date - birth</em>date in e-commerce.
  • Binning (Discretization)

Binning converts a continuous feature into an ordinal categorical one. This can help tree models that must make many splits to approximate a non-linear relationship, or linear models where the relationship is piecewise constant.

PYTHON
import pandas as pd
from sklearn.preprocessing import KBinsDiscretizer
import numpy as np

rng = np.random.default_rng(0)
X_age = rng.integers(18, 80, size=(500, 1))

# Quantile binning: equal-frequency bins
kbd = KBinsDiscretizer(n_bins=5, encode="ordinal", strategy="quantile")
print(pd.Series(kbd.fit_transform(X_age).flatten()).value_counts().sort_index())

Three strategies are available: uniform (equal-width), quantile (equal-frequency, robust to skew), and kmeans (bin boundaries determined by k-means). For downstream linear models, use encode=&quot;onehot&quot; to avoid imposing an ordinal assumption on the bins.

Datetime Features

Timestamps contain multiple independent signals that must be extracted explicitly:

PYTHON
import pandas as pd

df = pd.DataFrame({"ts": pd.to_datetime([
    "2024-01-15 08:30", "2024-07-04 14:00", "2024-12-25 23:00"
])})

df["year"] = df["ts"].dt.year
df["month"] = df["ts"].dt.month
df["dayofweek"] = df["ts"].dt.dayofweek   # 0 = Monday
df["hour"] = df["ts"].dt.hour
df["is_weekend"] = df["ts"].dt.dayofweek >= 5
df["quarter"] = df["ts"].dt.quarter
# Cyclical encoding to preserve periodicity
import numpy as np
df["hour_sin"] = np.sin(2 * np.pi * df["hour"] / 24)
df["hour_cos"] = np.cos(2 * np.pi * df["hour"] / 24)
print(df.head())

Cyclical encoding is important for periodic features such as hour-of-day, day-of-week, and month: the raw integer treats hour 23 and hour 0 as far apart, whereas they are adjacent. Mapping to $\sin$ and $\cos$ preserves the circular topology.

Text Features

For short categorical-like text fields (job titles, product descriptions), simple engineered features are often more robust than full NLP pipelines:

  • Character count, word count.
  • Presence of specific keywords (binary flags).
  • TF-IDF on a limited vocabulary.

For longer text, bag-of-words or TF-IDF with max_features capping is the standard first step before neural approaches.

Geographic Features

For latitude/longitude pairs, useful engineered features include:

  • Distance from a meaningful anchor point (nearest city center, warehouse, airport).
  • Administrative region (neighborhood, ZIP code, county) via reverse geocoding.
  • Interaction of lat/lon with time (rush-hour patterns differ by location).

The Haversine formula gives the great-circle distance between two points:

$$d = 2r \arcsin\left(\sqrt{\sin^2\!\left(\frac{\Delta\phi}{2}\right) + \cos\phi_1\cos\phi_2\sin^2\!\left(\frac{\Delta\lambda}{2}\right)}\right)$$

where $\phi$ is latitude and $\lambda$ is longitude in radians.

Code Examples

Cyclical Encoding for Hour-of-Day

Demonstrates why raw hour values mislead linear models and how sin/cos encoding fixes it.

PYTHON
import numpy as np
import matplotlib.pyplot as plt

hours = np.arange(24)
sin_h = np.sin(2 * np.pi * hours / 24)
cos_h = np.cos(2 * np.pi * hours / 24)

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
axes[0].plot(hours, hours)
axes[0].set_title("Raw hour — 23 and 0 appear far apart")
axes[0].set_xlabel("Hour")
axes[1].scatter(sin_h, cos_h, c=hours, cmap="hsv", s=80)
axes[1].set_title("Sin/cos encoding — 23 and 0 are adjacent")
axes[1].set_xlabel("sin(hour)")
axes[1].set_ylabel("cos(hour)")
cb = plt.colorbar(axes[1].collections[0], ax=axes[1])
cb.set_label("Hour")
plt.tight_layout()
plt.show()
Output
(Two matplotlib plots displayed)

Haversine Distance Feature

Computes the great-circle distance (km) between each row's location and a fixed anchor (NYC) as a feature.

PYTHON
import numpy as np
import pandas as pd

def haversine_km(lat1, lon1, lat2, lon2):
    R = 6371.0
    phi1, phi2 = np.radians(lat1), np.radians(lat2)
    dphi = np.radians(lat2 - lat1)
    dlam = np.radians(lon2 - lon1)
    a = np.sin(dphi / 2)**2 + np.cos(phi1) * np.cos(phi2) * np.sin(dlam / 2)**2
    return 2 * R * np.arcsin(np.sqrt(a))

df = pd.DataFrame({
    "city": ["New York", "Los Angeles", "Chicago", "Houston"],
    "lat": [40.71, 34.05, 41.88, 29.76],
    "lon": [-74.01, -118.24, -87.63, -95.37],
})

NYC = (40.71, -74.01)
df["dist_from_nyc_km"] = haversine_km(df["lat"], df["lon"], *NYC)
print(df[["city", "dist_from_nyc_km"]].round(1))
Output
          city  dist_from_nyc_km
0     New York               0.0
1  Los Angeles            3940.1
2      Chicago            1144.3
3      Houston            2279.9

Selective Polynomial Interactions with Regularization

Adds degree-2 interactions for two selected features and fits Ridge regression, showing the expanded feature names.

PYTHON
import numpy as np
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.linear_model import Ridge
from sklearn.pipeline import Pipeline

rng = np.random.default_rng(42)
X = rng.standard_normal((300, 2))
y = 2 * X[:, 0] + 0.5 * X[:, 0] * X[:, 1] - X[:, 1]**2 + rng.normal(scale=0.1, size=300)

pipe = Pipeline([
    ("poly", PolynomialFeatures(degree=2, include_bias=False)),
    ("scale", StandardScaler()),
    ("ridge", Ridge(alpha=1.0)),
])
pipe.fit(X, y)

feature_names = pipe.named_steps["poly"].get_feature_names_out(["x1", "x2"])
coefs = pipe.named_steps["ridge"].coef_
for name, coef in zip(feature_names, coefs):
    print(f"{name:10s}: {coef:+.4f}")
Output
x1        : +1.9982
x2        : -0.9972
x1^2      : -0.0041
x1 x2     : +0.4998
x2^2      : -1.0018

31.4 Building Reproducible Pipelines with scikit-learn Advanced

The Case for Pipelines

In a Jupyter notebook it is easy to preprocess a DataFrame, fit a model, and evaluate it — but this workflow is a leakage minefield and is nearly impossible to deploy. The problems are:

  1. Leakage during cross-validation. If you fit scalers or encoders before calling cross<em>val</em>score, the fold held out for evaluation has already influenced the preprocessing statistics.
  2. Reproducibility. A sequence of manual df[&#039;col&#039;] = ... mutations is hard to test, version, or re-run on new data.
  3. Deployment gap. Transformations must be replicated in production using exactly the same fitted parameters. Manual code invites drift.

scikit-learn's Pipeline and ColumnTransformer solve all three problems by encapsulating the entire preprocessing-modeling workflow in a single serializable object.

Pipeline Basics

A Pipeline chains estimators sequentially. Each intermediate step must implement fit and transform; the final step only needs fit (and predict for supervised tasks):

PYTHON
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression

pipe = Pipeline([
    ("scaler", StandardScaler()),
    ("clf", LogisticRegression(max_iter=500)),
])
# pipe.fit() calls scaler.fit_transform() then clf.fit()
# pipe.predict() calls scaler.transform() then clf.predict()

Crucially, cross<em>val</em>score(pipe, X, y) re-fits the entire pipeline on each training fold, so scalers never see the validation data. This is the single most important reason to use pipelines.

ColumnTransformer: Heterogeneous Data

Real DataFrames contain a mixture of numeric, categorical, and datetime columns that require different transformations. ColumnTransformer applies different transformers to different column subsets and concatenates the results:

PYTHON
import pandas as pd
import numpy as np
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import StandardScaler, OneHotEncoder, OrdinalEncoder
from sklearn.impute import SimpleImputer
from sklearn.pipeline import Pipeline
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.model_selection import cross_val_score

df = pd.DataFrame({
    "age": [25, 32, np.nan, 45, 52],
    "income": [50000, 80000, 62000, 120000, 95000],
    "education": ["high_school", "bachelors", "masters", "bachelors", "phd"],
    "city": ["NYC", "LA", "NYC", "Chicago", "LA"],
    "churned": [0, 0, 1, 0, 1],
})
X = df.drop(columns="churned")
y = df["churned"]

numeric_pipe = Pipeline([
    ("impute", SimpleImputer(strategy="median")),
    ("scale", StandardScaler()),
])

ordinal_pipe = Pipeline([
    ("enc", OrdinalEncoder(categories=[
        ["high_school", "bachelors", "masters", "phd"]
    ])),
])

col_tf = ColumnTransformer([
    ("num",     numeric_pipe,   ["age", "income"]),
    ("ord",     ordinal_pipe,   ["education"]),
    ("ohe",     OneHotEncoder(handle_unknown="ignore", sparse_output=False), ["city"]),
], remainder="drop")

full_pipe = Pipeline([
    ("preprocess", col_tf),
    ("model", GradientBoostingClassifier(n_estimators=100, random_state=0)),
])

print("Pipeline steps:", [step[0] for step in full_pipe.steps])
print("Output feature count:", col_tf.fit_transform(X).shape[1])

Accessing Feature Names After Transformation

Post-transformation feature names are essential for interpreting model coefficients and importances:

PYTHON
col_tf.fit(X)
print(col_tf.get_feature_names_out())

The output prefixes each name with the transformer name (e.g., num<strong>age, ohe</strong>city<em>NYC).

Serialization and Deployment

A fitted pipeline serializes cleanly with joblib (preferred over pickle for large numpy arrays):

PYTHON
import joblib
full_pipe.fit(X, y)
joblib.dump(full_pipe, "churn_model_v1.joblib")

# In the serving environment:
pipe_loaded = joblib.load("churn_model_v1.joblib")
new_row = pd.DataFrame([{"age": 29, "income": 70000,
                          "education": "masters", "city": "NYC"}])
print(pipe_loaded.predict_proba(new_row))

The serialized artifact captures all fitted scaler parameters, encoder vocabularies, and model weights in a single file. Versioning this file alongside its training code is the minimum viable MLOps practice.

Custom Transformers

Business logic often requires transformations not available in sklearn. Subclass BaseEstimator and TransformerMixin:

PYTHON
from sklearn.base import BaseEstimator, TransformerMixin

class LogTransformer(BaseEstimator, TransformerMixin):
    """Log1p-transform selected columns."""
    def __init__(self, columns):
        self.columns = columns

    def fit(self, X, y=None):
        return self  # stateless

    def transform(self, X):
        X = X.copy()
        X[self.columns] = np.log1p(X[self.columns])
        return X

By implementing fit and transform separately and keeping fit stateless where appropriate, this transformer works correctly inside cross</em>val_score with no modifications.

Code Examples

Full Pipeline with Target Encoding Inside Cross-Validation

Uses category_encoders.TargetEncoder inside a Pipeline to demonstrate that target statistics are computed only on training folds, preventing leakage.

PYTHON
import pandas as pd
import numpy as np
import category_encoders as ce
from sklearn.pipeline import Pipeline
from sklearn.linear_model import Ridge
from sklearn.model_selection import cross_val_score, KFold

rng = np.random.default_rng(123)
n = 500
categories = rng.choice(["A", "B", "C", "D", "E"], size=n)
# target has a real relationship with category
group_means = {"A": 1.0, "B": 2.0, "C": 3.0, "D": 4.0, "E": 5.0}
y = np.array([group_means[c] for c in categories]) + rng.normal(scale=0.5, size=n)
X = pd.DataFrame({"cat": categories})

pipe = Pipeline([
    ("te", ce.TargetEncoder(smoothing=5)),
    ("model", Ridge()),
])

cv = KFold(n_splits=5, shuffle=True, random_state=0)
scores = cross_val_score(pipe, X, y, cv=cv, scoring="r2")
print(f"CV R² scores: {scores.round(3)}")
print(f"Mean R²:      {scores.mean():.3f}")
Output
CV R² scores: [0.921 0.934 0.911 0.918 0.927]
Mean R²:      0.922

ColumnTransformer with get_feature_names_out

Builds a realistic heterogeneous preprocessor and prints all output feature names after fitting.

PYTHON
import pandas as pd
import numpy as np
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import StandardScaler, OneHotEncoder
from sklearn.impute import SimpleImputer
from sklearn.pipeline import Pipeline

df = pd.DataFrame({
    "sqft": [1200, 900, 1500, np.nan, 2000],
    "rooms": [3, 2, 4, 3, 5],
    "neighborhood": ["downtown", "suburb", "suburb", "downtown", "rural"],
    "condition": ["good", "fair", "excellent", "good", "fair"],
})

num_pipe = Pipeline([("imp", SimpleImputer()), ("sc", StandardScaler())])
cat_pipe = OneHotEncoder(sparse_output=False, handle_unknown="ignore")

ct = ColumnTransformer([
    ("num", num_pipe, ["sqft", "rooms"]),
    ("cat", cat_pipe, ["neighborhood", "condition"]),
])

ct.fit(df)
print("Output features:")
for name in ct.get_feature_names_out():
    print(" ", name)
Output
Output features:
  num__sqft
  num__rooms
  cat__neighborhood_downtown
  cat__neighborhood_rural
  cat__neighborhood_suburb
  cat__condition_excellent
  cat__condition_fair
  cat__condition_good

31.5 Feature Selection: Filter, Wrapper, and Embedded Methods Advanced

Why Select Features?

More features are not always better. Irrelevant or redundant features can:

  • Increase training and inference time.
  • Introduce noise that harms generalization (the curse of dimensionality is especially acute for distance-based methods).
  • Make models harder to interpret and audit.
  • Increase storage and serialization costs in production.

Feature selection answers the question: given a candidate set of features, which subset should we use? The three major families — filter, wrapper, and embedded — trade off computational cost against fidelity to the model being trained.

Filter Methods

Filter methods score features independently of any model, using a statistical measure of relevance to the target.

Variance Threshold

Removes features whose variance falls below a threshold — a cheap way to discard near-constant columns:

PYTHON
from sklearn.feature_selection import VarianceThreshold
import numpy as np

X = np.array([[0, 1, 0.5], [0, 0, 0.6], [0, 1, 0.4]])
sel = VarianceThreshold(threshold=0.01)
print(sel.fit_transform(X))  # drops column 0 (all zeros)

Mutual Information

Mutual information (MI) measures the nonlinear statistical dependence between a feature $X_j$ and target $Y$:

$$I(X_j; Y) = \sum_{x,y} p(x,y) \log \frac{p(x,y)}{p(x)p(y)}$$

Unlike correlation, MI detects non-linear relationships and interactions. SelectKBest with mutual<em>info</em>classif or mutual<em>info</em>regression selects the top $k$ features by MI score:

PYTHON
from sklearn.feature_selection import SelectKBest, mutual_info_regression
from sklearn.datasets import fetch_california_housing
import pandas as pd

data = fetch_california_housing()
X, y = data.data, data.target
features = data.feature_names

sel = SelectKBest(mutual_info_regression, k=4)
sel.fit(X, y)
scores = pd.Series(sel.scores_, index=features).sort_values(ascending=False)
print(scores.round(3))
print("\nSelected:", [features[i] for i in sel.get_support(indices=True)])

Correlation-Based Filtering

Remove features that are highly correlated with each other (redundancy) or have very low correlation with the target (irrelevance). A practical implementation:

PYTHON
import pandas as pd
import numpy as np

def drop_correlated(df, threshold=0.95):
    corr = df.corr().abs()
    upper = corr.where(np.triu(np.ones(corr.shape), k=1).astype(bool))
    to_drop = [col for col in upper.columns if any(upper[col] > threshold)]
    return df.drop(columns=to_drop), to_drop

Wrapper Methods

Wrapper methods evaluate feature subsets by actually training a model and measuring validation performance. They are model-aware but expensive.

Recursive Feature Elimination (RFE)

RFE fits the model, removes the feature(s) with the lowest importance, and repeats until the desired number of features is reached:

PYTHON
from sklearn.feature_selection import RFE
from sklearn.linear_model import Ridge
from sklearn.datasets import fetch_california_housing

data = fetch_california_housing()
X, y = data.data, data.target

rfe = RFE(estimator=Ridge(), n_features_to_select=4, step=1)
rfe.fit(X, y)
print(dict(zip(data.feature_names, rfe.ranking_)))
print("Selected:", [f for f, s in zip(data.feature_names, rfe.support_) if s])

RFECV extends this by choosing the optimal feature count via cross-validation, avoiding the need to specify n<em>features</em>to<em>select manually.

Embedded Methods

Embedded methods perform feature selection as part of model training. They are typically more computationally efficient than wrappers because the model is trained once.

Lasso (L1 Regularization)

The Lasso penalty shrinks coefficients of uninformative features exactly to zero:

$$\min_{\beta} \|y - X\beta\|_2^2 + \lambda \|\beta\|_1$$

PYTHON
from sklearn.linear_model import LassoCV
from sklearn.feature_selection import SelectFromModel
from sklearn.datasets import fetch_california_housing
import numpy as np

data = fetch_california_housing()
X, y = data.data, data.target

lasso_cv = LassoCV(cv=5, random_state=0).fit(X, y)
sel = SelectFromModel(lasso_cv, prefit=True)
X_reduced = sel.transform(X)
print(f"Original features: {X.shape[1]}, After Lasso selection: {X_reduced.shape[1]}")
print("Non-zero coefs:", np.sum(lasso_cv.coef_ != 0))

Tree-Based Feature Importance

Gradient boosted trees and random forests compute feature importances as the total impurity reduction attributed to each feature across all splits. SelectFromModel with a tree estimator is a fast, non-linear alternative to Lasso:

PYTHON
from sklearn.ensemble import RandomForestRegressor
from sklearn.feature_selection import SelectFromModel
from sklearn.datasets import fetch_california_housing
import pandas as pd

data = fetch_california_housing()
X, y = data.data, data.target

rf = RandomForestRegressor(n_estimators=100, random_state=0, n_jobs=-1)
rf.fit(X, y)
importances = pd.Series(rf.feature_importances_, index=data.feature_names)
print(importances.sort_values(ascending=False).round(4))

Pitfall: The default impurity-based importance is biased toward high-cardinality continuous features. Use permutation</em>importance for a model-agnostic, bias-free estimate at the cost of additional computation.

Feature Stores: Operationalizing Feature Engineering

A feature store is a centralized repository that stores, versions, and serves computed features for both model training and online inference. It solves three problems at scale:

  1. Training-serving skew: Features are computed once and stored; the exact same values are used in training and served at prediction time.
  2. Reuse: Features computed for one model are available to any other model, avoiding duplicate computation.
  3. Point-in-time correctness: When training on historical data, a feature store ensures only features available before each label's timestamp are used — preventing temporal leakage.

Popular options include Feast (open-source), Tecton, Hopsworks, and Vertex AI Feature Store. Even without a formal feature store, the same principles apply: keep feature computation code versioned and separate from model code, and never recompute features at serving time using data not available at the prediction timestamp.

Code Examples

RFECV: Automatic Feature Count Selection

Uses RFECV with Ridge regression to automatically determine the optimal number of features on the California housing dataset.

PYTHON
import numpy as np
import matplotlib.pyplot as plt
from sklearn.feature_selection import RFECV
from sklearn.linear_model import Ridge
from sklearn.datasets import fetch_california_housing
from sklearn.preprocessing import StandardScaler

data = fetch_california_housing()
X = StandardScaler().fit_transform(data.data)
y = data.target

rfecv = RFECV(
    estimator=Ridge(),
    step=1,
    cv=5,
    scoring="r2",
    min_features_to_select=1,
    n_jobs=-1,
)
rfecv.fit(X, y)

print(f"Optimal feature count: {rfecv.n_features_}")
print(f"Selected features: {[data.feature_names[i] for i in rfecv.get_support(indices=True)]}")

plt.plot(range(1, len(rfecv.cv_results_['mean_test_score']) + 1),
         rfecv.cv_results_['mean_test_score'])
plt.xlabel("Number of features")
plt.ylabel("CV R²")
plt.title("RFECV: Performance vs Feature Count")
plt.axvline(rfecv.n_features_, color="r", linestyle="--", label=f"Optimal = {rfecv.n_features_}")
plt.legend()
plt.show()
Output
Optimal feature count: 7
Selected features: ['MedInc', 'HouseAge', 'AveRooms', 'AveBedrms', 'Population', 'AveOccup', 'Latitude']

Permutation Importance vs. Impurity Importance

Compares impurity-based and permutation-based importances to illustrate bias toward high-cardinality features in the impurity method.

PYTHON
import numpy as np
import pandas as pd
from sklearn.ensemble import RandomForestClassifier
from sklearn.inspection import permutation_importance
from sklearn.model_selection import train_test_split

rng = np.random.default_rng(7)
n = 1000
# Two informative binary features, one high-cardinality uninformative feature
X = pd.DataFrame({
    "useful_A": rng.integers(0, 2, size=n),
    "useful_B": rng.integers(0, 2, size=n),
    "noise_hicardinal": rng.uniform(0, 1, size=n),  # continuous noise
})
y = (X["useful_A"] ^ X["useful_B"]).values

X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.3, random_state=0)

rf = RandomForestClassifier(n_estimators=200, random_state=0).fit(X_tr, y_tr)

impurity = pd.Series(rf.feature_importances_, index=X.columns)
perm = permutation_importance(rf, X_te, y_te, n_repeats=30, random_state=0)
perm_means = pd.Series(perm.importances_mean, index=X.columns)

print("Impurity importance:")
print(impurity.sort_values(ascending=False).round(4))
print("\nPermutation importance:")
print(perm_means.sort_values(ascending=False).round(4))
Output
Impurity importance:
noise_hicardinal    0.4521
useful_A            0.2751
useful_B            0.2728

Permutation importance:
useful_A            0.2183
useful_B            0.2171
noise_hicardinal   -0.0011

SelectFromModel with Lasso in a Pipeline

Embeds Lasso-based feature selection inside a Pipeline so selection is re-fit on each CV fold.

PYTHON
import numpy as np
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LassoCV, Ridge
from sklearn.feature_selection import SelectFromModel
from sklearn.model_selection import cross_val_score
from sklearn.datasets import fetch_california_housing

data = fetch_california_housing()
X, y = data.data, data.target

pipe = Pipeline([
    ("scale", StandardScaler()),
    ("lasso_sel", SelectFromModel(LassoCV(cv=3, random_state=0))),
    ("ridge", Ridge()),
])

scores = cross_val_score(pipe, X, y, cv=5, scoring="r2")
print(f"CV R²: {scores.round(3)}")
print(f"Mean: {scores.mean():.3f}")
Output
CV R²: [0.614 0.601 0.589 0.633 0.617]
Mean: 0.611