04 · Feature Engineering

Feature engineering is the place where domain knowledge becomes competitive advantage. Every Kaggle competition analysis that asks “why did the winner win?” arrives at the same answer: their features were better. This is the work that algorithms cannot do for you — transforming raw data into representations that expose the signal your model needs. It is neither fully automatable nor fully systematizable, but its patterns are learnable.

⚠️ What Most People Get Wrong: They treat feature engineering as something you do before modeling, once. The real workflow is iterative: model → error analysis → hypothesis about missing signal → new feature → retrain. The best features often come from studying where your model fails, not from running automated feature generation on the full dataset before you’ve fit a single model.

---

The Feature Engineering Mindset

Features are hypotheses. When you create a feature, you’re claiming: “the relationship between this transformed representation and the target is simpler (more linear, more monotone, more separable) than the relationship between the raw input and the target.” This is testable.

The budget for feature engineering effort is not uniform:

• High ROI: Encoding high-cardinality categoricals correctly, treating temporal features properly, engineering domain-specific interactions the model cannot discover implicitly

• Low ROI (often negative): Blind polynomial feature generation on high-dimensional data, automated feature crossing without domain grounding

---

1. Handling Missing Data

Missing data is not a uniform problem. The mechanism matters.

Types of Missingness

Type	Definition	Example	Implication
MCAR (Missing Completely At Random)	Probability of missing is independent of observed and unobserved values	Sensor randomly fails	Imputation safe; any method works
MAR (Missing At Random)	Probability of missing depends on observed (but not unobserved) values	Income missing more for young people, and age is observed	Imputation works if you condition on observed predictors
MNAR (Missing Not At Random)	Probability of missing depends on the missing value itself	High-income people hide their income	Imputation is biased; requires domain modeling

Imputation Strategies

import pandas as pd
import numpy as np
from sklearn.impute import SimpleImputer, KNNImputer
from sklearn.experimental import enable_iterative_imputer
from sklearn.impute import IterativeImputer

# Strategy 1: Simple imputation (fast, often fine for MCAR)
simple = SimpleImputer(strategy='median')  # or 'mean', 'most_frequent', 'constant'
X_simple = simple.fit_transform(X_train)

# Strategy 2: KNN imputation (uses similar rows to fill in)
# k=5 is a reasonable default; computationally O(n²) for large datasets
knn_imp = KNNImputer(n_neighbors=5)
X_knn = knn_imp.fit_transform(X_train)

# Strategy 3: Iterative imputation (MICE — multiple imputation by chained equations)
# Fits a model for each feature with missing values on the others
iter_imp = IterativeImputer(max_iter=10, random_state=42)
X_iter = iter_imp.fit_transform(X_train)

# Strategy 4: Add missingness indicator (critical for MNAR)
# The FACT that a value is missing is itself informative
def add_missingness_indicators(df):
    """For each column with missing values, add a binary indicator."""
    result = df.copy()
    missing_cols = df.columns[df.isnull().any()].tolist()
    for col in missing_cols:
        result[f'{col}_is_missing'] = df[col].isnull().astype(int)
    return result

df_with_indicators = add_missingness_indicators(df)

The Missingness Decision Tree

Is the feature missing > 60% of the time?
    YES → Drop the feature (likely uninformative)
    NO → Is missingness random (MCAR/MAR)?
            YES → SimpleImputer (median for numeric, mode for categorical)
            NO → Add missingness indicator + impute
Is downstream model tree-based?
    YES → Consider leaving NaN if sklearn/LightGBM/XGBoost handles it natively
    NO → Must impute before fitting

---

2. Encoding Categorical Variables

The cardinal sin: treating categorical integers as if they have ordinal meaning (encoding “Red”=1, “Green”=2, “Blue”=3 implies Red < Green < Blue).

Encoding Strategy by Cardinality

import pandas as pd
from sklearn.preprocessing import LabelEncoder, OrdinalEncoder
from category_encoders import TargetEncoder, BinaryEncoder
import numpy as np

# --- Binary categoricals (2 unique values) ---
# Just 0/1, nothing to discuss.

# --- Low-cardinality nominale (3-15 unique values) ---
# Use One-Hot Encoding
df_encoded = pd.get_dummies(df, columns=['color', 'region'], drop_first=True)
# drop_first=True to avoid perfect multicollinearity (dummy variable trap)
# sklearn equivalent for pipelines:
from sklearn.preprocessing import OneHotEncoder
ohe = OneHotEncoder(handle_unknown='ignore', sparse_output=False, drop='first')

# --- Ordinal categoricals (have natural order) ---
# "Low" < "Medium" < "High" — preserve the order
from sklearn.preprocessing import OrdinalEncoder
ord_enc = OrdinalEncoder(
    categories=[['Low', 'Medium', 'High']],  # explicit order
    handle_unknown='use_encoded_value',
    unknown_value=-1
)

# --- High-cardinality nominals (> 15 unique values, e.g., zip codes, product IDs) ---
# One-Hot Encoding explodes dimensionality: 50,000 zip codes → 50,000 columns
# Solution: Target Encoding (encode with mean of target)
# CRITICAL: Must use cross-validation to prevent target leakage

from category_encoders import TargetEncoder
target_enc = TargetEncoder(smoothing=1.0)  # smoothing regularizes rare categories
# Use only in a pipeline with CV to avoid overfitting to target

# Alternative: Frequency encoding (replace category with its frequency)
def frequency_encode(df, col):
    freq_map = df[col].value_counts(normalize=True).to_dict()
    df[f'{col}_freq'] = df[col].map(freq_map)
    return df

# Alternative: Binary encoding (good balance for 10-1000 categories)
from category_encoders import BinaryEncoder
bin_enc = BinaryEncoder()  # encodes k categories in ceil(log2(k)) binary features

Target Encoding Done Right (The Leakage-Safe Version)

from sklearn.model_selection import KFold
import numpy as np

def safe_target_encode(X_train, y_train, X_test, col, n_splits=5, smoothing=10):
    """
    Target encoding with out-of-fold computation to prevent leakage.
    The global mean acts as a prior; smoothing controls strength.
    """
    X_train = X_train.copy()
    X_test = X_test.copy()
    
    global_mean = y_train.mean()
    oof_encoded = np.zeros(len(X_train))
    
    kf = KFold(n_splits=n_splits, shuffle=True, random_state=42)
    
    for train_idx, val_idx in kf.split(X_train):
        X_tr, X_v = X_train.iloc[train_idx], X_train.iloc[val_idx]
        y_tr = y_train.iloc[train_idx]
        
        # Compute encoding on training fold
        stats = y_tr.groupby(X_tr[col]).agg(['mean', 'count'])
        # Smoothed estimate: (count * category_mean + smoothing * global_mean) / (count + smoothing)
        stats['smoothed'] = (stats['count'] * stats['mean'] + smoothing * global_mean) / \
                           (stats['count'] + smoothing)
        
        oof_encoded[val_idx] = X_v[col].map(stats['smoothed']).fillna(global_mean)
    
    X_train[f'{col}_target_enc'] = oof_encoded
    
    # For test set: use full training data statistics
    all_stats = y_train.groupby(X_train[col]).agg(['mean', 'count'])
    all_stats['smoothed'] = (all_stats['count'] * all_stats['mean'] + smoothing * global_mean) / \
                            (all_stats['count'] + smoothing)
    X_test[f'{col}_target_enc'] = X_test[col].map(all_stats['smoothed']).fillna(global_mean)
    
    return X_train, X_test

---

3. Feature Scaling

Tree-based models are invariant to monotone transformations of features. Distance-based and gradient-based models are not. Forgetting to scale before SVM, k-NN, or linear models is a correctness error, not a style preference.

from sklearn.preprocessing import StandardScaler, MinMaxScaler, RobustScaler, PowerTransformer
import numpy as np

# StandardScaler: (x - mean) / std → zero mean, unit variance
# Use for: linear models, SVMs, neural networks, PCA
std_scaler = StandardScaler()

# MinMaxScaler: (x - min) / (max - min) → [0, 1]
# Use for: image pixel values, when bounded range matters, neural nets
# Problem: sensitive to outliers
minmax = MinMaxScaler()

# RobustScaler: (x - median) / IQR → robust to outliers
# Use for: data with known outliers that you don't want to clip
robust = RobustScaler()

# PowerTransformer: Box-Cox or Yeo-Johnson → makes features more Gaussian
# Use for: skewed features before linear models or regression targets
power = PowerTransformer(method='yeo-johnson')  # yeo-johnson handles negatives

# Which models need scaling?
scaling_necessity = {
    'LinearRegression': 'Not strictly required (no regularization), but helps numerics',
    'Ridge/Lasso/ElasticNet': 'REQUIRED — regularization penalizes large weights',
    'LogisticRegression': 'REQUIRED — regularization + gradient descent convergence',
    'SVM': 'REQUIRED — kernel = dot product, scale-dependent',
    'k-NN': 'REQUIRED — distance is scale-dependent',
    'Neural Networks': 'REQUIRED — weight initialization assumes ~unit variance input',
    'PCA': 'REQUIRED — maximizes variance, large-scale features dominate',
    'DecisionTree': 'Not needed — splits are scale-invariant',
    'RandomForest': 'Not needed — ensemble of trees',
    'XGBoost/LightGBM': 'Not needed — tree-based',
    'k-Means': 'REQUIRED — uses Euclidean distance',
}
for model, need in scaling_necessity.items():
    print(f"{model:30s}: {need}")

Handling Skewed Distributions

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

def analyze_skewness(X, feature_names):
    """Identify features needing transformation."""
    results = []
    for i, name in enumerate(feature_names):
        col = X[:, i]
        skewness = stats.skew(col)
        results.append({'feature': name, 'skewness': skewness,
                       'action': 'log1p' if abs(skewness) > 1 else 'none'})
    
    return pd.DataFrame(results).sort_values('skewness', key=abs, ascending=False)

# Common transformations for right-skewed numeric features
# House prices, income, counts → often right-skewed
log_transformed = np.log1p(X_skewed)   # log(1+x), handles zeros
sqrt_transformed = np.sqrt(X_skewed)    # gentler than log
cbrt_transformed = np.cbrt(X_skewed)    # cube root, handles negatives

---

4. Feature Selection

Not all features help. Irrelevant features add noise; correlated features waste capacity and can destabilize gradient descent. Feature selection is regularization applied at the data level.

Method 1: Filter Methods (Fast, Model-Agnostic)

from sklearn.feature_selection import (
    SelectKBest, f_classif, mutual_info_classif,
    VarianceThreshold, chi2
)
import pandas as pd

# Remove near-zero variance features (almost constant — carry no information)
var_thresh = VarianceThreshold(threshold=0.01)
X_var = var_thresh.fit_transform(X_train)

# Correlation matrix — drop one of any pair with |corr| > 0.95
def remove_correlated_features(X_df, threshold=0.95):
    corr_matrix = X_df.corr().abs()
    upper = corr_matrix.where(np.triu(np.ones(corr_matrix.shape), k=1).astype(bool))
    to_drop = [col for col in upper.columns if any(upper[col] > threshold)]
    print(f"Dropping {len(to_drop)} correlated features: {to_drop}")
    return X_df.drop(columns=to_drop)

# Mutual Information (non-parametric, captures nonlinear relationships)
mi_scores = mutual_info_classif(X_train, y_train, random_state=42)
mi_df = pd.DataFrame({'feature': feature_names, 'mi_score': mi_scores})
mi_df = mi_df.sort_values('mi_score', ascending=False)
print(mi_df.head(20))

# Select top-k features by MI
selector_mi = SelectKBest(mutual_info_classif, k=20)
X_selected = selector_mi.fit_transform(X_train, y_train)
selected_features = [feature_names[i] for i in selector_mi.get_support(indices=True)]

Method 2: Wrapper Methods (Better, Expensive)

from sklearn.feature_selection import RFE, RFECV
from sklearn.ensemble import RandomForestClassifier

# Recursive Feature Elimination with CV (gold standard for small datasets)
rfecv = RFECV(
    estimator=RandomForestClassifier(n_estimators=100, random_state=42),
    step=1,
    cv=StratifiedKFold(5),
    scoring='roc_auc',
    min_features_to_select=5,
    n_jobs=-1
)
rfecv.fit(X_train, y_train)

print(f"Optimal number of features: {rfecv.n_features_}")
selected_mask = rfecv.support_
selected_features_rfe = [feature_names[i] for i, s in enumerate(selected_mask) if s]

# Plot CV scores vs number of features
plt.figure(figsize=(10, 5))
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 AUC')
plt.title('RFECV: Feature Count vs Performance')
plt.savefig('rfecv_curve.png', dpi=150)

Method 3: Embedded Methods (Best Balance of Speed and Quality)

from sklearn.linear_model import LassoCV
from sklearn.ensemble import RandomForestClassifier

# LASSO: L1 regularization drives irrelevant features to exactly zero
lasso = LassoCV(cv=5, random_state=42, max_iter=5000)
lasso.fit(X_train_scaled, y_train_reg)  # For regression targets

# Features with non-zero coefficients
lasso_features = [feature_names[i] for i, c in enumerate(lasso.coef_) if c != 0]
print(f"LASSO selected {len(lasso_features)} of {len(feature_names)} features")

# Tree-based feature importance (fast, works for classification/regression)
rf = RandomForestClassifier(n_estimators=200, random_state=42, n_jobs=-1)
rf.fit(X_train, y_train)

importance_df = pd.DataFrame({
    'feature': feature_names,
    'importance': rf.feature_importances_
}).sort_values('importance', ascending=False)

print(importance_df.head(20))

# WARNING: RF importance is biased toward high-cardinality features
# Use permutation importance for more reliable estimates:
from sklearn.inspection import permutation_importance

perm_imp = permutation_importance(
    rf, X_val, y_val, n_repeats=20, random_state=42, n_jobs=-1
)
perm_df = pd.DataFrame({
    'feature': feature_names,
    'importance': perm_imp.importances_mean,
    'std': perm_imp.importances_std
}).sort_values('importance', ascending=False)

---

5. Feature Interactions and Creation

Domain-Driven Interactions

Before automating, think about what interactions your domain suggests.

import pandas as pd
import numpy as np

# --- Ratio features (often powerful) ---
df['price_per_sqft'] = df['price'] / (df['sqft'] + 1)
df['debt_to_income'] = df['total_debt'] / (df['income'] + 1)
df['clicks_per_impression'] = df['clicks'] / (df['impressions'] + 1)

# --- Difference features ---
df['age_at_purchase'] = df['purchase_year'] - df['birth_year']
df['days_since_last_login'] = df['current_date'] - df['last_login_date']

# --- Count features (from grouped aggregations) ---
df['user_total_orders'] = df.groupby('user_id')['order_id'].transform('count')
df['user_avg_order_value'] = df.groupby('user_id')['order_value'].transform('mean')
df['category_popularity'] = df.groupby('category')['clicks'].transform('sum')

# --- Temporal features from datetime ---
def extract_temporal_features(df, datetime_col):
    dt = pd.to_datetime(df[datetime_col])
    df[f'{datetime_col}_hour'] = dt.dt.hour
    df[f'{datetime_col}_dayofweek'] = dt.dt.dayofweek  # 0=Monday, 6=Sunday
    df[f'{datetime_col}_month'] = dt.dt.month
    df[f'{datetime_col}_quarter'] = dt.dt.quarter
    df[f'{datetime_col}_is_weekend'] = (dt.dt.dayofweek >= 5).astype(int)
    df[f'{datetime_col}_is_month_start'] = dt.dt.is_month_start.astype(int)
    df[f'{datetime_col}_is_month_end'] = dt.dt.is_month_end.astype(int)
    # Cyclical encoding for hour (23:59 is close to 00:00)
    df[f'{datetime_col}_hour_sin'] = np.sin(2 * np.pi * dt.dt.hour / 24)
    df[f'{datetime_col}_hour_cos'] = np.cos(2 * np.pi * dt.dt.hour / 24)
    return df

Automated Polynomial Features (Use Cautiously)

from sklearn.preprocessing import PolynomialFeatures

# degree=2 on 100 features → 5,051 features (n_features*(n_features+3)/2 + 1)
# degree=2 on 1000 features → 501,501 features — DO NOT DO THIS
# Limit to:
# 1. Pre-selected important features only (top 10-20)
# 2. Low-dimensional datasets (< 30 features)

top_features = importance_df.head(10)['feature'].tolist()
X_top = df[top_features]

poly = PolynomialFeatures(degree=2, interaction_only=True, include_bias=False)
X_poly = poly.fit_transform(X_top)
poly_feature_names = poly.get_feature_names_out(top_features)
print(f"Polynomial features: {X_top.shape[1]} → {X_poly.shape[1]}")

---

6. The sklearn Pipeline — Feature Engineering at Scale

from sklearn.pipeline import Pipeline, FeatureUnion
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import StandardScaler, OneHotEncoder
from sklearn.impute import SimpleImputer
from sklearn.ensemble import GradientBoostingClassifier

# Define column types
numeric_features = ['age', 'income', 'tenure_days', 'num_orders']
categorical_features = ['country', 'product_category', 'channel']

# Numeric pipeline
numeric_transformer = Pipeline(steps=[
    ('imputer', SimpleImputer(strategy='median')),
    ('scaler', StandardScaler())
])

# Categorical pipeline
categorical_transformer = Pipeline(steps=[
    ('imputer', SimpleImputer(strategy='most_frequent')),
    ('encoder', OneHotEncoder(handle_unknown='ignore', sparse_output=False))
])

# Combine
preprocessor = ColumnTransformer(
    transformers=[
        ('num', numeric_transformer, numeric_features),
        ('cat', categorical_transformer, categorical_features)
    ]
)

# Full pipeline: preprocessing + model
full_pipeline = Pipeline(steps=[
    ('preprocessor', preprocessor),
    ('classifier', GradientBoostingClassifier(n_estimators=200, random_state=42))
])

# This pipeline does everything correctly:
# 1. fit() on training data only — no leakage
# 2. transform() on validation/test using training statistics
# 3. Serializable with joblib — ready for production
full_pipeline.fit(X_train_df, y_train)
val_auc = roc_auc_score(y_val, full_pipeline.predict_proba(X_val_df)[:, 1])
print(f"Validation AUC: {val_auc:.4f}")

---

7. Feature Engineering for Specific Data Types

Text Features (Quick Reference)

from sklearn.feature_extraction.text import TfidfVectorizer

# TF-IDF: term frequency * inverse document frequency
tfidf = TfidfVectorizer(
    max_features=10000,
    ngram_range=(1, 2),  # unigrams and bigrams
    min_df=5,            # ignore terms in < 5 documents
    max_df=0.95,         # ignore terms in > 95% of documents
    sublinear_tf=True    # apply log normalization to TF
)
X_text = tfidf.fit_transform(text_column)

# For production: embedding-based features beat TF-IDF for semantics
# from sentence_transformers import SentenceTransformer
# model = SentenceTransformer('all-MiniLM-L6-v2')  # 384-dim embeddings
# X_emb = model.encode(text_column.tolist())

Geographic Features

# Encode latitude/longitude with Haversine distance to reference points
from math import radians, sin, cos, sqrt, atan2

def haversine_km(lat1, lon1, lat2, lon2):
    """Distance between two lat/lon points in km."""
    R = 6371  # Earth radius in km
    dlat, dlon = radians(lat2-lat1), radians(lon2-lon1)
    a = sin(dlat/2)**2 + cos(radians(lat1)) * cos(radians(lat2)) * sin(dlon/2)**2
    return 2 * R * atan2(sqrt(a), sqrt(1-a))

# Distance to nearest city center, airport, etc. as features
df['dist_to_nyc'] = df.apply(
    lambda r: haversine_km(r['lat'], r['lon'], 40.7128, -74.0060), axis=1
)

---

Practice Problems

Problem 1 — Encoding Benchmark: On a dataset with a high-cardinality categorical column (e.g., “city” with 500+ unique values), compare: OneHotEncoding, TargetEncoding (with CV), FrequencyEncoding, and LabelEncoding. Measure AUC for a logistic regression on each. Report which wins and why, including the overfitting risk of target encoding.

Problem 2 — Missing Data Mechanism Test: Take a complete dataset. Artificially introduce MCAR, MAR, and MNAR missingness (20% each). Apply simple median imputation to all three. Compare model performance to the ground truth (no missing data). Show quantitatively why MNAR imputation fails.

Problem 3 — Feature Importance Audit: Train a Random Forest. Compare MDI (default) importance vs. permutation importance. Find at least one feature where they disagree significantly. Hypothesize why (cardinality, correlation) and verify.

Problem 4 — Pipeline Leakage Hunt: Write a deliberately broken pipeline that scales the full dataset before splitting. Write the correct version. Compare validation AUC. Report the leakage gap (typically 0.01–0.05 AUC on tabular data with a scaler, more with target encoding).

---

Resources

Resource	Type	Why	Time
sklearn Preprocessing Guide	Docs	Canonical, comprehensive	Reference
Géron Ch. 2 (End-to-End ML Project)	Book	Best practical walkthrough of FE	3 hrs
“Feature Engineering for ML” — Alice Zheng	Book	Only book dedicated to this topic	6 hrs
Kaggle Learn: Feature Engineering	Course	Practical, free, competition-tested	4 hrs
StatQuest: Missing Data video	Video	MCAR/MAR/MNAR explained clearly	20 min

---

Return to README.md · Previous: 03_model_selection_and_evaluation.md · Next: 05_classical_ml_in_production.md