Feature Engineering and Feature Selection

Feature engineering changes the representation a model receives; feature selection chooses which inputs to retain. A product or distance can make a relationship accessible to a simple model. Removing inputs can reduce estimation noise or measurement costs, but selection can also discard useful information. The examples below examine these choices with fixed prediction tasks and held-out evaluation.

Interactions a linear model cannot see

A linear score adds weighted input columns, for example \(b+w_1x_1+w_2x_2\). Its weights cannot create the product \(x_1x_2\) unless that product is supplied as another column. Logistic regression applies a sigmoid to the score, so its probability is nonlinear, but its decision boundary in the two raw inputs is still a straight line. The following task requires different labels in alternating quadrants.

The examples use NumPy and scikit-learn. cross_val_score returns held-out accuracy for the classifiers and held-out R² for the regressors used here; cv=5 uses five folds. A Pipeline fits each preprocessing step with the model inside a training fold. A fixed row-wise formula such as a product or squared radius needs no fitted statistics.

import numpy as np
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures

rng = np.random.default_rng(0)
n = 4000
X = rng.normal(size=(n, 2))
y = (X[:, 0] * X[:, 1] > 0).astype(int)        # depends on the product

plain = cross_val_score(LogisticRegression(), X, y, cv=5).mean()
with_interactions = cross_val_score(
    Pipeline([("p", PolynomialFeatures(2, include_bias=False)),
              ("m", LogisticRegression())]), X, y, cv=5).mean()
print(f"raw features {plain:.4f}   with interactions {with_interactions:.4f}")
# raw features 0.6377   with interactions 0.9918

The label is 1 when the coordinates share a sign: \((2,3)\) and \((-2,-3)\) both have product 6, while \((2,-3)\) has product −6 and label 0. No single straight line separates these quadrants. The raw-input model scores 0.6377 in this run, while the polynomial version scores 0.9918. More data does not remove the straight-line restriction, although both scores and their gap can change with the sample. The engineered representation contains the correct boundary; this does not guarantee that a finite, regularized fit finds it exactly.

With two inputs, PolynomialFeatures(2, include_bias=False) returns \(x_1,x_2,x_1^2,x_1x_2,x_2^2\). Thus the experiment adds squares as well as the needed interaction. For \(p\) inputs and maximum degree \(d\), the number of output columns is \(\binom{p+d}{d}-1\) when the constant column is excluded. That gives 1,325 columns for 50 inputs at degree 2 and 23,425 at degree 3. This rapid growth affects memory and estimation; compare a broad expansion with targeted terms that the task suggests.

A radial feature for a circular boundary

The next example uses a circular decision boundary. An axis-aligned tree ensemble can approximate it with many splits. A radius-based feature gives a linear classifier a more direct representation of the same geometry.

import numpy as np
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import cross_val_score

rng = np.random.default_rng(1)
n = 3000
angle = rng.uniform(0, 2 * np.pi, n)
radius = rng.uniform(0, 2, n)
X = np.column_stack([radius * np.cos(angle), radius * np.sin(angle)])
y = (radius > 1).astype(int)                   # a circular boundary

print("raw x,y    logistic", round(cross_val_score(LogisticRegression(), X, y, cv=5).mean(), 4))
print("raw x,y    forest  ", round(cross_val_score(
    RandomForestClassifier(n_estimators=100, random_state=0), X, y, cv=5).mean(), 4))
X_plus = np.column_stack([X, (X ** 2).sum(1)])       # add the squared radius
print("+ radius^2 logistic", round(cross_val_score(LogisticRegression(), X_plus, y, cv=5).mean(), 4))
# raw x,y    logistic 0.5157
# raw x,y    forest   0.9887
# + radius^2 logistic 0.9977

Logistic regression on the raw coordinates scores 0.5157 and the random forest 0.9887. Adding \(r^2=x_1^2+x_2^2\) raises logistic regression to 0.9977 in this run. The target boundary is simply \(r^2=1\): a point at \((0.6,0.8)\) lies on it because \(0.6^2+0.8^2=1\). The code retains both original coordinates, so it fits three coefficients plus an intercept. A separate model using only squared radius would have one coefficient and an intercept, but is not the model tested here.

This target has no label noise: every label is determined by the radius. The remaining prediction errors therefore cannot be assigned to irreducible noise. A finite forest approximates the circle using axis-aligned partitions; the regularized logistic fit estimates its boundary from a finite sample. The printed accuracies do not locate the errors or separate these sources of error. This is a comparison of these configurations on one generated dataset, not a general ranking of feature engineering and forests.

A mechanism involving a ratio, difference, distance, or elapsed time suggests a candidate feature. Define it using information available at prediction time, handle cases such as a zero denominator, and assess it on held-out data. Ordinary regression trees predict constants within leaves, so extrapolation can be difficult; an engineered feature can help represent a relevant relationship without guaranteeing accurate extrapolation. If candidate transforms are chosen using validation results, that search is part of model selection too.

The three families of selection

import numpy as np
from sklearn.feature_selection import SelectKBest, f_regression, RFE
from sklearn.linear_model import Ridge, LassoCV

rng = np.random.default_rng(2)
n, p = 400, 40
X = rng.normal(size=(n, p))
X[:, 1] = X[:, 0] + 0.01 * rng.normal(size=n)  # feature 1 is a noisy copy of feature 0
y = 2 * X[:, 0] + 1.5 * X[:, 5] + rng.normal(size=n)

k = SelectKBest(f_regression, k=5).fit(X, y)
print("filter  :", np.sort(np.where(k.get_support())[0]))
r = RFE(Ridge(), n_features_to_select=5).fit(X, y)
print("wrapper :", np.sort(np.where(r.support_)[0]))
lasso = LassoCV(cv=5, random_state=0).fit(X, y)
print("embedded:", np.sort(np.where(np.abs(lasso.coef_) > 1e-6)[0]))
# filter  : [ 0  1  5  8 37]
# wrapper : [ 0  1  5 18 30]
# embedded: [ 0  5 18 30 36]

Features 0 and 5 enter the generating equation directly. Feature 1 also predicts the target because it is a noisy copy of feature 0; it is a redundant proxy, not pure noise. All three methods retain 0 and 5 here. The filter and wrapper also retain the proxy, while this lasso fit drops it. The remaining 37 columns were generated independently of the target mechanism.

SelectKBest(f_regression) ranks each column by its marginal linear association with the target and retains five. RFE repeatedly fits Ridge and removes columns with low coefficient importance until five remain. LassoCV chooses an \(\ell_1\) penalty by internal cross-validation, then refits; coefficients near zero are treated as unselected here. Its selected set can change with the penalty, feature scales, and sample. Lasso does not invariably drop duplicates: splitting a positive coefficient between two identical columns leaves their total absolute-value penalty unchanged. Near-duplicate columns can also exchange roles across fits.

The filter and wrapper retain two independently generated noise columns, while lasso retains three in this run. The first two procedures are required to select five columns even though there are only three signal-associated candidates. Chance sample correlations and the selection rule both matter. These printed sets are descriptive fits on all rows, not estimates of predictive performance or proof of variable recovery. When evaluating a selected model, put selection inside each training fold; tune the selection settings without using the outer evaluation fold.

Cardinality on real data

Cardinality is the number of distinct levels in a categorical field. The Adult census extract has eight categorical columns, whose widths and level frequencies affect encoding choices. This descriptive audit may download data on its first run; it does not select an encoding based on test performance.

import numpy as np
from fairlearn.datasets import fetch_adult
from sklearn.preprocessing import OneHotEncoder

d = fetch_adult(as_frame=True)
X = d.data.copy()
cats = X.select_dtypes(exclude=[np.number]).columns

for c in cats:
    print(f"  {c:16s} {X.nunique(dropna=True):3d} non-missing levels")
print("non-missing category total:", sum(X.nunique(dropna=True) for c in cats))
enc = OneHotEncoder(handle_unknown="ignore").fit(X[cats])
width = sum(len(levels) for levels in enc.categories_)
print("one-hot width including missing categories:", width)
print("width with six numeric columns:", width + X.select_dtypes(include=[np.number]).shape[1])
counts = X["native-country"].value_counts()
print("countries with fewer than 100 rows:", int((counts < 100).sum()))
print("smallest country count:", int(counts.min()))
#   workclass          8 non-missing levels
#   education         16 non-missing levels
#   marital-status     7 non-missing levels
#   occupation        14 non-missing levels
#   relationship       6 non-missing levels
#   race               5 non-missing levels
#   sex                2 non-missing levels
#   native-country    41 non-missing levels
# non-missing category total: 99
# one-hot width including missing categories: 102
# width with six numeric columns: 108
# countries with fewer than 100 rows: 26
# smallest country count: 1

The categorical columns have 99 observed non-missing levels in total. With a separate missing category in workclass, occupation, and native-country, this encoder produces 102 columns. Keeping the six numeric columns gives 108 inputs altogether. These are full-dataset audit counts; a training-fold encoder can have fewer levels.

Of the 41 observed country levels, 26 have fewer than 100 rows, and the smallest has one. That makes some category-specific estimates uncertain even though the total dataset has 48,842 rows. Sparse one-hot storage can handle many zeros efficiently, but does not create information about rare levels. A rare level may be absent from a training fold; it is not guaranteed to be. Define and test an unknown-category policy for that case.

Rare-level grouping and target encoding are candidates to compare with one-hot encoding. Fit any frequency threshold or grouping rule on training data. Target encoding estimates an outcome summary per category, often shrinking rare-category estimates toward a global mean. For evaluation, the mapping must exclude held-out labels. Training rows also need protection against encoding their own labels, especially for rare levels: scikit-learn TargetEncoder uses internal cross-fitting in fit_transform. Fitting a mapping on all training labels and then transforming those same rows does not provide that protection. Grouped or time-dependent data require a compatible splitting design.

Exercises

1. The cost of a feature that carries nothing. Hold three real features fixed and add 0, 10, 50, 200, and 1000 pure-noise columns. Report the cross-validated \(R^2\) at each.

Compare the feature count with the number of rows in each training fold, and keep track of which regularization setting is held fixed.

import numpy as np
from sklearn.linear_model import Ridge
from sklearn.model_selection import cross_val_score

rng = np.random.default_rng(0)
n = 300
signal = rng.normal(size=(n, 3))
y = signal @ [2.0, -1.0, 1.5] + rng.normal(size=n)

for extra in (0, 10, 50, 200, 1000):
    X = np.hstack([signal, rng.normal(size=(n, extra))]) if extra else signal
    r2 = cross_val_score(Ridge(alpha=1.0), X, y, cv=5).mean()
    print(f"{extra:5d} noise features  R2 {r2:7.4f}   total p {X.shape[1]}")
#     0 noise features  R2  0.8568   total p 3
#    10 noise features  R2  0.8466   total p 13
#    50 noise features  R2  0.8319   total p 53
#   200 noise features  R2  0.2459   total p 203
#  1000 noise features  R2  0.0664   total p 1003

2. Selection is unstable under resampling. Run a filter selection on 200 bootstrap resamples and report how often each feature is chosen. Two features carry signal and one duplicates another.

Compare the stability of the redundant proxy with that of the remaining selected columns. Does stability imply that a feature is necessary?

import numpy as np
from collections import Counter
from sklearn.feature_selection import SelectKBest, f_regression

rng = np.random.default_rng(1)
n, p = 250, 60
X = rng.normal(size=(n, p))
X[:, 1] = X[:, 0] + 0.05 * rng.normal(size=n)  # near-duplicate of feature 0
y = 2 * X[:, 0] + 1.5 * X[:, 7] + rng.normal(size=n)

counts = Counter()
for _ in range(200):
    idx = rng.integers(0, n, n)
    chosen = SelectKBest(f_regression, k=4).fit(X[idx], y[idx]).get_support()
    counts.update(np.where(chosen)[0].tolist())

for f, c in counts.most_common(8):
    print(f"feature {f:3d} selected {c / 200:.3f} of bootstraps")
# feature   0 selected 1.000 of bootstraps
# feature   1 selected 1.000 of bootstraps
# feature   7 selected 1.000 of bootstraps
# feature  19 selected 0.315 of bootstraps
# feature  31 selected 0.265 of bootstraps
# feature  38 selected 0.090 of bootstraps
# feature  18 selected 0.040 of bootstraps
# feature  11 selected 0.030 of bootstraps

3. When binning helps and when it costs. Discretize a continuous feature into 12 bins and compare against the raw feature, on a target that is linear in it and on a target that is a step function of it.

Compare the fitted prediction shapes and the held-out scores for the two targets.

import numpy as np
from sklearn.linear_model import Ridge
from sklearn.model_selection import cross_val_score
from sklearn.preprocessing import KBinsDiscretizer
from sklearn.pipeline import make_pipeline

rng = np.random.default_rng(2)
n = 3000
x = rng.uniform(-3, 3, n)
y_linear = 2 * x + 0.5 * rng.normal(size=n)
y_step = (np.floor(x) % 2 == 0).astype(float) * 2 + 0.5 * rng.normal(size=n)

for name, y in (("linear truth", y_linear), ("step truth", y_step)):
    raw = cross_val_score(Ridge(), x.reshape(-1, 1), y, cv=5).mean()
    binner = KBinsDiscretizer(n_bins=12, encode="onehot-dense", strategy="uniform")
    binned = cross_val_score(make_pipeline(binner, Ridge()), x.reshape(-1, 1), y, cv=5).mean()
    print(f"{name:14s} raw R2 {raw:7.4f}   binned R2 {binned:7.4f}")
# linear truth   raw R2  0.9799   binned R2  0.9734
# step truth     raw R2  0.0625   binned R2  0.7937