Multicollinearity

Good news: Tree-based models (Random Forest, XGBoost) are mostly immune to multicollinearity — they just pick the most useful feature at each split. But linear models (linear regression, logistic regression) are highly sensitive to it, and linear models remain extremely common in business analytics.

import pandas as pd
import numpy as np
from scipy import stats

# 12 houses, 8 features, target = price in £000s
df = pd.DataFrame({
    'price':         [285,320,195,410,265,375,220,340,290,430,255,310],
    'sq_footage':    [850,980,620,1200,780,1100,660,1020,870,1250,720,940],
    'num_rooms':     [4,  5,  3,  6,  4,  5,  3,  5,  4,  6,  3,  5],
    'num_bathrooms': [1,  2,  1,  3,  1,  2,  1,  2,  2,  3,  1,  2],
    'age_years':     [25, 8,  42, 3,  31, 12, 38, 6,  18, 2,  45, 15],
    'dist_city_km':  [4.2,2.1,8.5,1.2,5.8,1.8,7.3,2.5,3.9,0.9,9.1,3.2],
    'garden_sqm':    [45, 80, 20, 120,35, 95, 25, 75, 50, 130,15, 60],
    'crime_score':   [42, 28, 65, 18, 51, 22, 71, 31, 38, 15, 78, 35],
    'school_km':     [0.8,1.2,2.1,0.5,1.5,0.9,2.8,1.1,0.7,0.4,3.2,1.3]
})

features = ['sq_footage','num_rooms','num_bathrooms','age_years',
            'dist_city_km','garden_sqm','crime_score','school_km']

print(df[features].head(3))

What just happened?

pandas stores the data. Even before any analysis you can sense the overlaps — sq_footage, num_rooms, and num_bathrooms all probably grow together. dist_city_km and crime_score might also move in tandem (central locations typically have lower crime). Let's confirm those suspicions properly.

# Build the correlation matrix — every feature vs every other feature
corr = df[features].corr()

print("Suspicious feature pairs (|r| > 0.80):\n")

for i in range(len(features)):
    for j in range(i + 1, len(features)):   # upper triangle only — avoids listing each pair twice
        r = corr.iloc[i, j]
        if abs(r) > 0.80:
            level = "HIGH" if abs(r) > 0.90 else "MODERATE"
            print(f"  {level:9}  {features[i]}  x  {features[j]}:  r = {r:.3f}")

What just happened?

pandas' .corr() builds the correlation matrix and we loop through the upper triangle to list every suspicious pair once. 14 suspect pairs from 8 features — this dataset has a serious problem. The worst: sq_footage x num_rooms at r=0.982 (almost perfectly correlated) and dist_city_km x crime_score at r=0.975.

def calculate_vif(dataframe, feature_cols):
    """
    For each feature, measure how well other features can predict it.
    VIF formula: 1 / (1 - R2)
    When R2 is high (others predict this feature well), VIF becomes very large.
    """
    results = []
    for col in feature_cols:
        y = dataframe[col].values
        max_r2 = 0
        for other in feature_cols:
            if other == col:
                continue
            # stats.linregress fits a line from 'other' to 'col' and gives us R
            _, _, r, _, _ = stats.linregress(dataframe[other].values, y)
            if r**2 > max_r2:
                max_r2 = r**2   # keep the highest R2 found from any single predictor

        vif = round(1 / (1 - max_r2), 1) if max_r2 < 1.0 else 999
        flag = "DROP  ✗" if vif > 10 else ("REVIEW ⚠" if vif > 5 else "OK    ✓")
        results.append({'feature': col, 'VIF': vif, 'verdict': flag})

    return pd.DataFrame(results).sort_values('VIF', ascending=False)

print(calculate_vif(df, features).to_string(index=False))

What just happened?

scipy's stats.linregress() fits a line between any two columns and gives us R (correlation). Squaring it gives R2 — how much of the variation is explained. We find the highest R2 any single other feature achieves, then apply the VIF formula: 1 / (1 − R2). The closer R2 is to 1, the more the feature is redundant, and the higher the VIF.

Seven out of eight features score above 10. Only age_years (VIF=3.7) is genuinely independent. This is a severe multicollinearity problem — feeding all 8 features into a linear model would be a disaster.

# Measure each feature's correlation with the target (price)
target_corrs = {}
for feat in features:
    r, _ = stats.pearsonr(df[feat], df['price'])
    target_corrs[feat] = abs(r)   # absolute value — strength only

# For every suspicious pair, mark the weaker one (lower correlation with price) for removal
drop_set = set()
for i in range(len(features)):
    for j in range(i + 1, len(features)):
        if abs(corr.iloc[i, j]) > 0.80:
            fa, fb = features[i], features[j]
            # Drop whichever has the weaker connection to the target
            weaker = fb if target_corrs[fa] >= target_corrs[fb] else fa
            drop_set.add(weaker)

final = [f for f in features if f not in drop_set]

print("Features to KEEP:")
for f in final:
    print(f"  {f:<18}  r with price = {target_corrs[f]:.3f}")

print("\nFeatures to DROP (redundant):")
for f in sorted(drop_set):
    print(f"  {f}")

What just happened?

scipy's stats.pearsonr() gives each feature a target-correlation score. We loop through the suspicious pairs, compare those scores, and mark the weaker one for removal. From 8 features down to 3 — but these 3 features cover the most important non-overlapping signals. sq_footage (r=0.983) is the dominant predictor, garden_sqm (r=0.920) adds something else, and age_years (r=0.712, VIF=3.7) is clean and independent.

y = df['price'].values

def fit_linear(X_data, feat_names):
    """Fit OLS linear regression and print coefficients."""
    X_b = np.column_stack([np.ones(len(X_data)), X_data])   # add intercept column of 1s
    coeffs = np.linalg.lstsq(X_b, y, rcond=None)[0][1:]     # skip intercept, keep features
    for name, coef in zip(feat_names, coeffs):
        print(f"  {name:<18}  coefficient = {coef:>+.3f}")

print("=== ALL 8 FEATURES (multicollinearity present) ===")
fit_linear(df[features].values, features)

print("\n=== 3 CLEAN FEATURES (multicollinearity removed) ===")
clean = ['sq_footage', 'age_years', 'garden_sqm']
fit_linear(df[clean].values, clean)

Did you spot the problem?

Look at crime_score in the 8-feature model — it has a positive coefficient (+1.122). The model is claiming: higher crime = higher price. That is completely wrong. We know from the data that crime is negatively correlated with price (r = −0.893).

This happened because crime_score and dist_city_km are almost identical (r=0.975). The model got confused trying to split their effects and accidentally reversed the crime coefficient's sign. This is the textbook symptom of multicollinearity. The 3-feature model has sensible, interpretable coefficients throughout — age gets a negative coefficient (older houses are cheaper) and sq_footage gets a positive one (bigger houses cost more). Both make intuitive sense.

Teacher's Note

Multicollinearity doesn't always hurt predictions — but it always hurts interpretation. A model with redundant features might still make decent predictions overall. But you can't trust individual coefficients. If your manager asks "does crime rate affect house prices in our model?" — you can't answer from a multicollinear model because the coefficients are lying.

In business analytics, explanation usually matters as much as accuracy. Stakeholders make decisions based on model insights, not just predictions. Fix multicollinearity before you try to explain your model to anyone.

Up Next · Lesson 26

Visualising Distributions

Histograms, KDE curves, and box plots — build the visuals that make distribution shape instantly obvious to any audience.