Why can R-squared be a misleading metric for model comparison, and how does adjusted R-squared address this?

R-squared is defined as R² = 1 - SS_res/SS_tot, where SS_res is the sum of squared residuals and SS_tot is the total sum of squares (variance of y). It represents the proportion of variance in y explained by the model. The mathematical issue is that R² is monotonically non-decreasing as you add more features to a linear model — even adding a completely random, uninformative feature cannot decrease R², because OLS will always find some coefficient for it that fits the training data at least marginally better (or, at worst, exactly zero, leaving R² unchanged).

This makes R² unsuitable for comparing models with different numbers of features, since it will always favor the more complex model even if the added complexity is just overfitting noise. Adjusted R² corrects for this by introducing a penalty for the number of predictors p: R²_adj = 1 - (1-R²)(n-1)/(n-p-1). This formula can actually decrease when an added feature doesn't improve the fit enough to outweigh the penalty for the added degree of freedom, making it a fairer basis for comparing models with different feature counts.

import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score

def adjusted_r2(r2, n, p):
    return 1 - (1 - r2) * (n - 1) / (n - p - 1)

# Demonstrating R^2 always increasing with more (even random) features
n_samples = 100
X_real = np.random.randn(n_samples, 3)
y = X_real @ [2, -1, 0.5] + np.random.randn(n_samples) * 0.5

model1 = LinearRegression().fit(X_real, y)
r2_1 = r2_score(y, model1.predict(X_real))

# Add 10 completely random, uninformative columns
X_noisy = np.column_stack([X_real, np.random.randn(n_samples, 10)])
model2 = LinearRegression().fit(X_noisy, y)
r2_2 = r2_score(y, model2.predict(X_noisy))

print(f'R^2 with 3 features:  {r2_1:.4f}')
print(f'R^2 with 13 features: {r2_2:.4f}')  # always >= r2_1, even with noise!

adj_r2_1 = adjusted_r2(r2_1, n_samples, 3)
adj_r2_2 = adjusted_r2(r2_2, n_samples, 13)
print(f'Adjusted R^2 (3 features):  {adj_r2_1:.4f}')
print(f'Adjusted R^2 (13 features): {adj_r2_2:.4f}')  # often lower!