Feature Engineering

Raw features are rarely the best representation of your data. For this dataset, the EDA revealed three types of work to do: cap extreme outliers, create interaction features that capture relationships the raw columns miss, and engineer geographic signals that tree models can exploit.

Start with a clean copy of the data.

from sklearn.datasets import fetch_california_housing
import pandas as pd
import numpy as np
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler

housing = fetch_california_housing(as_frame=True)
df = housing.frame.copy()

print(df.shape)  # (20640, 9)

Step 1 — Handle Outliers in Room and Occupancy Columns

The EDA showed extreme outliers in AveRooms, AveBedrms, and AveOccup. These come from group-quarters blocks (prisons, dorms, military housing). Capping at the 99th percentile removes the distortion without discarding rows.

def cap_at_percentile(series, upper_pct=99):
    """Cap values at the given percentile (Winsorization)."""
    upper = series.quantile(upper_pct / 100)
    return series.clip(upper=upper)

cols_to_cap = ['AveRooms', 'AveBedrms', 'AveOccup']

for col in cols_to_cap:
    before_max = df[col].max()
    df[col] = cap_at_percentile(df[col], upper_pct=99)
    after_max = df[col].max()
    print(f"{col}: max {before_max:.1f} → {after_max:.1f}")

# AveRooms:  max 141.9 → 12.4
# AveBedrms: max  34.1 →  2.5
# AveOccup:  max 1243.3 →  8.5

Warning

Do not drop outlier rows. That reduces your dataset size for no reason. Capping preserves all 20,640 rows while removing the distortion. The 1% of extreme values are still represented — just not at their original magnitude.

Tip

Capping at the 99th percentile (also called Winsorization) is standard practice for aggregated data like census block averages. For individual transaction data, you might investigate outliers more carefully before capping.

Step 2 — Interaction Features

Individual columns tell part of the story. Ratios tell the rest.

rooms_per_person

AveRooms / AveOccup captures spaciousness — how many rooms are available per occupant. A block with 6 rooms and 2 occupants is very different from a block with 6 rooms and 6 occupants.

df['rooms_per_person'] = df['AveRooms'] / df['AveOccup']

# Sanity check
print(df['rooms_per_person'].describe().round(2))
#  count    20640.00
#  mean         2.19
#  std          0.89
#  min          0.22
#  max         14.67

bedrooms_ratio

AveBedrms / AveRooms gives the fraction of rooms that are bedrooms. Lower ratio = more common areas (living rooms, studies, kitchens) relative to sleeping space, which often signals larger, more expensive homes.

df['bedrooms_ratio'] = df['AveBedrms'] / df['AveRooms']

print(df['bedrooms_ratio'].describe().round(3))
# mean    0.213
# std     0.059
# min     0.100
# max     0.667

income_per_room

MedInc / AveRooms captures purchasing power relative to housing supply. A high-income block with few rooms available is a strong signal for high house values.

df['income_per_room'] = df['MedInc'] / df['AveRooms']

print(df['income_per_room'].describe().round(3))
# mean    0.784
# std     0.382
# min     0.043
# max     4.745

Info

Why create these features instead of just letting the model learn the relationship?

Linear models cannot learn multiplicative or ratio relationships on their own. If you give a linear model AveRooms and AveOccup separately, it can only learn a weighted sum — it cannot compute a ratio. Interaction features encode that ratio explicitly and make it available to any model.

Tree-based models can approximate ratios through sequential splits, but pre-computing the ratio makes the signal clearer and reduces the number of splits the tree needs to find it.

Step 3 — Log Transforms for Skewed Features

Population and AveOccup (after capping) are still right-skewed. Log-transforming them makes their distributions more symmetric, which helps linear models.

df['log_population'] = np.log1p(df['Population'])
df['log_ave_occup'] = np.log1p(df['AveOccup'])

print(f"Population skew:     raw={df['Population'].skew():.2f},  log={df['log_population'].skew():.2f}")
print(f"AveOccup skew:       raw={df['AveOccup'].skew():.2f},   log={df['log_ave_occup'].skew():.2f}")
# Population skew:     raw=6.17,  log=0.42
# AveOccup skew:       raw=3.64,  log=0.33

Tip

Use np.log1p(x) instead of np.log(x) when values can be zero. log1p(x) computes log(1 + x), which is defined at x=0. For this dataset values are always positive, but log1p is a good habit.

Step 4 — Geographic Features

Latitude and longitude carry strong spatial signal, but linear models cannot use them effectively as raw coordinates. Two approaches work well.

Approach A — Distance to Coast (Simplified)

A simplified coast-distance feature using a few reference points captures the coastal premium without requiring GIS libraries.

def simplified_coast_distance(lat, lon):
    """
    Approximate distance to the California coast using reference points.
    Returns distance in degrees (not exact km, but preserves the signal).
    Reference points: SF Bay (37.8, -122.4), LA (33.9, -118.4), San Diego (32.7, -117.2)
    """
    coast_points = [
        (37.8, -122.4),   # San Francisco Bay
        (36.6, -121.9),   # Monterey
        (34.4, -119.7),   # Santa Barbara
        (33.9, -118.4),   # Los Angeles
        (32.7, -117.2),   # San Diego
    ]
    min_dist = float('inf')
    for c_lat, c_lon in coast_points:
        dist = np.sqrt((lat - c_lat) ** 2 + (lon - c_lon) ** 2)
        if dist < min_dist:
            min_dist = dist
    return min_dist

df['coast_distance'] = df.apply(
    lambda row: simplified_coast_distance(row['Latitude'], row['Longitude']),
    axis=1
)

print(df['coast_distance'].describe().round(3))
# mean    1.428
# std     1.043
# min     0.001
# max     6.185

Approach B — Geographic Cluster Labels

KMeans on coordinates groups blocks into geographic regions. This gives tree-based models a categorical-style feature that encodes which part of California a block belongs to.

geo_features = df[['Latitude', 'Longitude']].values

kmeans = KMeans(n_clusters=6, random_state=42, n_init=10)
df['geo_cluster'] = kmeans.fit_predict(geo_features)

# Check cluster sizes
print(df['geo_cluster'].value_counts().sort_index())
# 0    4851
# 1    3109
# 2    2947
# 3    3762
# 4    2839
# 5    3132

# Visualize the clusters
import matplotlib.pyplot as plt

colors = ['#e41a1c','#377eb8','#4daf4a','#984ea3','#ff7f00','#a65628']
for cluster_id in range(6):
    mask = df['geo_cluster'] == cluster_id
    plt.scatter(
        df.loc[mask, 'Longitude'],
        df.loc[mask, 'Latitude'],
        c=colors[cluster_id], s=1, alpha=0.4,
        label=f'Cluster {cluster_id}'
    )
plt.legend(markerscale=5, loc='upper right')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.title('Geographic Clusters (k=6)')
plt.tight_layout()
plt.show()

Info

You do not have to choose between these two approaches. Use both. The cluster label gives a discrete region signal. The coast distance gives a continuous proximity signal. Together they let the model capture both regional and coastal effects.

Step 5 — Log-Transform the Target (Optional)

Whether to log-transform the target depends on which model you are building.

# If using linear regression: log-transform the target
df['log_target'] = np.log1p(df['MedHouseVal'])

print(f"Target skew:     raw={df['MedHouseVal'].skew():.3f}, log={df['log_target'].skew():.3f}")
# Target skew:     raw=1.079, log=0.313

For tree-based models, skip this step. Use MedHouseVal directly as the target.

Warning

If you train a model on a log-transformed target, you must invert the transform before computing MAE and RMSE in the original units. Forgetting this step makes your error metrics uninterpretable — a MAE of 0.12 in log space sounds great, but you need to convert it back to dollar units to know if it is actually good.

# Correct workflow with log target
log_preds = model.predict(X_test)
preds_original = np.expm1(log_preds)   # inverse of log1p
mae = mean_absolute_error(y_test_original, preds_original)

Step 6 — Prepare Final Feature Matrix

# Feature matrix for tree-based models (no scaling needed)
feature_cols_tree = [
    'MedInc', 'HouseAge', 'AveRooms', 'AveBedrms', 'Population',
    'AveOccup', 'Latitude', 'Longitude',
    'rooms_per_person', 'bedrooms_ratio', 'income_per_room',
    'log_population', 'log_ave_occup',
    'coast_distance', 'geo_cluster'
]

X = df[feature_cols_tree].copy()
y = df['MedHouseVal'].copy()

print(X.shape)   # (20640, 15)
print(y.shape)   # (20640,)

Step 7 — Scaling for Linear Models

Linear and regularised models (Ridge, Lasso) require features on the same scale. Tree-based models do not.

from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split

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

# Fit scaler on training data only — never on the full dataset
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)   # transform only, do not re-fit

print(X_train_scaled.mean(axis=0).round(3))  # approximately zero for all features
print(X_train_scaled.std(axis=0).round(3))   # approximately one for all features

Warning

Fit the scaler on training data only. If you fit on the full dataset before splitting, information from the test set leaks into your preprocessing step. This inflates your evaluation metrics and gives you a falsely optimistic view of model performance.

Feature Engineering Summary

Feature	Type	Rationale
rooms_per_person	Interaction	Captures spaciousness per occupant
bedrooms_ratio	Interaction	Fraction of rooms that are bedrooms; signals home type
income_per_room	Interaction	Purchasing power relative to available space
log_population	Log transform	Reduces right skew; helps linear models
log_ave_occup	Log transform	Reduces right skew after outlier capping
coast_distance	Geographic	Approximate proximity to the California coastline
geo_cluster	Geographic	Discrete geographic region label from KMeans

Success

After feature engineering you have 15 features instead of 8. The added features encode domain knowledge — that spaciousness, income-to-space ratio, and coastal proximity all matter for house prices — in a form that models can use directly.

---

← Back: EDA | Next: Model Building →