R-squared and Adjusted R²

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56

import numpy as np
 
def r_squared_from_scratch(y_true: np.ndarray, y_pred: np.ndarray) -> dict:
    """
    Compute R² and its components from scratch.
    
    Returns detailed breakdown of sum of squares.
    """
    y_true = np.asarray(y_true)
    y_pred = np.asarray(y_pred)
    n = len(y_true)
    
    # Mean of true values (baseline prediction)
    y_mean = np.mean(y_true)
    
    # Sum of Squares components
    ss_tot = np.sum((y_true - y_mean) ** 2)      # Total variance
    ss_res = np.sum((y_true - y_pred) ** 2)      # Residual (unexplained)
    ss_reg = np.sum((y_pred - y_mean) ** 2)      # Explained by model
    
    # R-squared
    r_squared = 1 - (ss_res / ss_tot)
    
    # Alternative formulation (same result for OLS)
    r_squared_alt = ss_reg / ss_tot
    
    return {
        'ss_tot': ss_tot,
        'ss_res': ss_res,
        'ss_reg': ss_reg,
        'r_squared': r_squared,
        'r_squared_alt': r_squared_alt,
        'variance_explained_pct': r_squared * 100,
        'mse': ss_res / n
    }
 
# Example: Predicting test scores
y_true = np.array([65, 72, 80, 85, 90, 75, 70, 88, 95, 78])
y_pred = np.array([68, 70, 82, 84, 88, 77, 72, 85, 92, 80])
 
# Calculate R²
result = r_squared_from_scratch(y_true, y_pred)
 
print("=== R² Breakdown ===")
print(f"Mean of y: {np.mean(y_true):.2f}")
print(f"
Sum of Squares:")
print(f"  SS_total (variance around mean): {result['ss_tot']:.2f}")
print(f"  SS_residual (unexplained):       {result['ss_res']:.2f}")
print(f"  SS_regression (explained):       {result['ss_reg']:.2f}")
print(f"
R² = 1 - (SS_res / SS_tot)")
print(f"R² = 1 - ({result['ss_res']:.2f} / {result['ss_tot']:.2f})")
print(f"R² = {result['r_squared']:.4f}")
print(f"
→ Model explains {result['variance_explained_pct']:.1f}% of variance in test scores")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45

import numpy as np
 
def r_squared_vs_accuracy():
    """
    Demonstrate that high R² doesn't mean predictions are 'close'.
    """
    np.random.seed(42)
    
    # Scenario 1: High R² with large absolute errors
    y_true_1 = np.array([100, 200, 300, 400, 500, 600, 700, 800, 900, 1000])
    y_pred_1 = y_true_1 + np.array([-30, 20, -40, 35, -25, 30, -35, 40, -20, 25])
    
    ss_tot_1 = np.sum((y_true_1 - np.mean(y_true_1)) ** 2)
    ss_res_1 = np.sum((y_true_1 - y_pred_1) ** 2)
    r2_1 = 1 - ss_res_1 / ss_tot_1
    mae_1 = np.mean(np.abs(y_true_1 - y_pred_1))
    
    # Scenario 2: Lower R² with smaller absolute errors
    y_true_2 = np.array([100, 102, 104, 106, 108, 110, 112, 114, 116, 118])
    y_pred_2 = y_true_2 + np.array([-3, 2, -4, 3, -2, 3, -3, 4, -2, 2])
    
    ss_tot_2 = np.sum((y_true_2 - np.mean(y_true_2)) ** 2)
    ss_res_2 = np.sum((y_true_2 - y_pred_2) ** 2)
    r2_2 = 1 - ss_res_2 / ss_tot_2
    mae_2 = np.mean(np.abs(y_true_2 - y_pred_2))
    
    print("=== R² vs Prediction Accuracy ===")
    print("
Scenario 1: High variance target (100 to 1000)")
    print(f"  R²: {r2_1:.4f}")
    print(f"  MAE: {mae_1:.1f}")
    print(f"  → High R² but predictions off by ~{mae_1:.0f} on average!")
    
    print("
Scenario 2: Low variance target (100 to 118)")
    print(f"  R²: {r2_2:.4f}")
    print(f"  MAE: {mae_2:.1f}")
    print(f"  → Lower R² but predictions off by only ~{mae_2:.0f}!")
    
    print("
*** Key Insight ***")
    print("R² depends on target variance. Same MAE gives different R²")
    print("depending on how spread out the true values are.")
 
r_squared_vs_accuracy()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56

import numpy as np
 
def geometric_r_squared(y_true: np.ndarray, y_pred: np.ndarray):
    """
    Demonstrate the geometric interpretation of R².
    R² = cos²(θ) between centered y and y_hat vectors.
    """
    # Center the vectors
    y_centered = y_true - np.mean(y_true)
    y_hat_centered = y_pred - np.mean(y_pred)
    
    # Compute angle using dot product
    # cos(θ) = (a · b) / (||a|| × ||b||)
    dot_product = np.dot(y_centered, y_hat_centered)
    norm_y = np.linalg.norm(y_centered)
    norm_y_hat = np.linalg.norm(y_hat_centered)
    
    cos_theta = dot_product / (norm_y * norm_y_hat)
    theta_radians = np.arccos(np.clip(cos_theta, -1, 1))
    theta_degrees = np.degrees(theta_radians)
    
    # R² from geometry
    r_squared_geometric = cos_theta ** 2
    
    # R² from standard formula
    ss_tot = np.sum((y_true - np.mean(y_true)) ** 2)
    ss_res = np.sum((y_true - y_pred) ** 2)
    r_squared_standard = 1 - ss_res / ss_tot
    
    # Correlation coefficient
    correlation = np.corrcoef(y_true, y_pred)[0, 1]
    
    print("=== Geometric R² Interpretation ===")
    print(f"
Vector norms:")
    print(f"  ||y - ȳ||  = {norm_y:.4f}")
    print(f"  ||ŷ - ȳ|| = {norm_y_hat:.4f}")
    print(f"
Angle between centered vectors:")
    print(f"  cos(θ) = {cos_theta:.4f}")
    print(f"  θ = {theta_degrees:.2f}°")
    print(f"
R² calculations:")
    print(f"  cos²(θ) = {r_squared_geometric:.4f}")
    print(f"  1 - SS_res/SS_tot = {r_squared_standard:.4f}")
    print(f"  Corr(y, ŷ)² = {correlation**2:.4f}")
    print(f"
→ All three methods give the same R²!")
    
    return r_squared_geometric
 
# Example
y_true = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
y_pred = np.array([1.2, 2.1, 2.8, 4.2, 5.1, 5.9, 7.2, 7.8, 9.1, 10.2])
 
geometric_r_squared(y_true, y_pred)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51

import numpy as np
 
def demonstrate_negative_r_squared():
    """
    Show scenarios where R² becomes negative.
    """
    # True values
    y_true = np.array([10, 20, 30, 40, 50])
    y_mean = np.mean(y_true)  # = 30
    
    print("True values:", y_true)
    print(f"Mean: {y_mean}")
    print(f"
Baseline (mean) predictions: [{y_mean}] × 5")
    
    ss_tot = np.sum((y_true - y_mean) ** 2)
    print(f"SS_tot (variance around mean): {ss_tot}")
    
    # Scenario 1: Good predictions
    y_pred_good = np.array([12, 22, 28, 38, 52])
    ss_res_good = np.sum((y_true - y_pred_good) ** 2)
    r2_good = 1 - ss_res_good / ss_tot
    print(f"
--- Good Model ---")
    print(f"Predictions: {y_pred_good}")
    print(f"SS_res: {ss_res_good}")
    print(f"R²: {r2_good:.4f} (positive)")
    
    # Scenario 2: Bad predictions (worse than mean)
    y_pred_bad = np.array([50, 10, 50, 10, 50])  # Wrong direction!
    ss_res_bad = np.sum((y_true - y_pred_bad) ** 2)
    r2_bad = 1 - ss_res_bad / ss_tot
    print(f"
--- Bad Model ---")
    print(f"Predictions: {y_pred_bad}")
    print(f"SS_res: {ss_res_bad}")
    print(f"R²: {r2_bad:.4f} (NEGATIVE!)")
    print(f"
→ Model is worse than just predicting the mean!")
    
    # Scenario 3: Predicting exactly the mean
    y_pred_mean = np.full_like(y_true, y_mean, dtype=float)
    ss_res_mean = np.sum((y_true - y_pred_mean) ** 2)
    r2_mean = 1 - ss_res_mean / ss_tot
    print(f"
--- Baseline Model (predict mean) ---")
    print(f"Predictions: {y_pred_mean}")
    print(f"SS_res: {ss_res_mean} (equals SS_tot)")
    print(f"R²: {r2_mean:.4f} (exactly zero)")
 
demonstrate_negative_r_squared()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51

import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score
 
def demonstrate_r_squared_inflation():
    """
    Show how R² increases with model complexity,
    even for useless random features.
    """
    np.random.seed(42)
    
    # True relationship: y = 2*x1 + 3*x2 + noise
    n_samples = 100
    X_true = np.random.randn(n_samples, 2)
    y = 2 * X_true[:, 0] + 3 * X_true[:, 1] + np.random.randn(n_samples) * 0.5
    
    print("True model: y = 2*x1 + 3*x2 + noise")
    print(f"Samples: {n_samples}
")
    
    # Progressively add random (useless) features
    r_squared_values = []
    
    for n_noise_features in range(0, 80, 5):
        # Create dataset with true features + noise features
        X_noise = np.random.randn(n_samples, n_noise_features)
        X_full = np.hstack([X_true, X_noise]) if n_noise_features > 0 else X_true
        
        # Fit model and compute R²
        model = LinearRegression()
        model.fit(X_full, y)
        y_pred = model.predict(X_full)
        
        r2 = r2_score(y, y_pred)
        r_squared_values.append((2 + n_noise_features, r2))
        
        if n_noise_features in [0, 10, 30, 50, 70]:
            print(f"Features: {2 + n_noise_features:2d} ({n_noise_features} noise) → R² = {r2:.4f}")
    
    print(f"
*** Key Insight ***")
    print("R² keeps increasing even though noise features add NOTHING!")
    print("With enough features, R² → 1.0 (perfect fit to training noise)")
 
demonstrate_r_squared_inflation()
# Example output:
# Features:  2 (0 noise) → R² = 0.9721
# Features: 12 (10 noise) → R² = 0.9841
# Features: 32 (30 noise) → R² = 0.9944
# Features: 52 (50 noise) → R² = 0.9982
# Features: 72 (70 noise) → R² = 0.9996

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52

import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score
 
def adjusted_r_squared(r2: float, n: int, p: int) -> float:
    """
    Compute adjusted R² from R², sample size, and number of predictors.
    
    Parameters:
    -----------
    r2 : float - R-squared value
    n : int - Number of samples
    p : int - Number of predictors (excluding intercept)
    """
    if n - p - 1 <= 0:
        raise ValueError("n must be greater than p + 1")
    
    return 1 - ((1 - r2) * (n - 1) / (n - p - 1))
 
def compare_r2_adjusted():
    """
    Compare R² and Adjusted R² as features increase.
    """
    np.random.seed(42)
    
    n_samples = 100
    X_true = np.random.randn(n_samples, 2)
    y = 2 * X_true[:, 0] + 3 * X_true[:, 1] + np.random.randn(n_samples) * 0.5
    
    print("Effect of adding noise features on R² vs Adjusted R²")
    print(f"{'Features':^10} | {'R²':^10} | {'Adj R²':^10} | {'Difference':^10}")
    print("-" * 48)
    
    for n_noise in [0, 5, 10, 20, 30, 50, 70, 90]:
        X_noise = np.random.randn(n_samples, max(1, n_noise))
        X = np.hstack([X_true, X_noise]) if n_noise > 0 else X_true
        
        n_features = X.shape[1]
        
        model = LinearRegression()
        model.fit(X, y)
        r2 = r2_score(y, model.predict(X))
        adj_r2 = adjusted_r_squared(r2, n_samples, n_features)
        
        print(f"{n_features:^10} | {r2:^10.4f} | {adj_r2:^10.4f} | {r2 - adj_r2:^10.4f}")
    
    print("
*** Key Insight ***")
    print("R² keeps increasing, but Adjusted R² DECREASES when")
    print("useless features are added—correctly penalizing overfitting.")
 
compare_r2_adjusted()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59

import numpy as np
from sklearn.metrics import r2_score
 
def demonstrate_pitfalls():
    """
    Show scenarios where R² is misleading.
    """
    np.random.seed(42)
    
    # Pitfall 1: High R² with non-linear relationship
    print("=== Pitfall 1: Non-linearity hidden by R² ===")
    x = np.linspace(0, 10, 100)
    y_nonlinear = np.sin(x) * x + np.random.randn(100) * 0.5
    
    # Fit linear model
    from sklearn.linear_model import LinearRegression
    lr = LinearRegression()
    lr.fit(x.reshape(-1, 1), y_nonlinear)
    y_pred = lr.predict(x.reshape(-1, 1))
    
    r2 = r2_score(y_nonlinear, y_pred)
    print(f"True relationship: y = sin(x) * x")
    print(f"Linear model R²: {r2:.4f}")
    print("→ R² looks decent, but the model is fundamentally wrong!")
    
    # Pitfall 2: Range restriction
    print("
=== Pitfall 2: Range restriction ===")
    x_full = np.linspace(0, 100, 1000)
    y_full = 2 * x_full + np.random.randn(1000) * 10
    
    # Full range
    lr.fit(x_full.reshape(-1, 1), y_full)
    r2_full = r2_score(y_full, lr.predict(x_full.reshape(-1, 1)))
    
    # Restricted range (50-60 only)
    mask = (x_full >= 50) & (x_full <= 60)
    x_restricted = x_full[mask]
    y_restricted = y_full[mask]
    lr.fit(x_restricted.reshape(-1, 1), y_restricted)
    r2_restricted = r2_score(y_restricted, lr.predict(x_restricted.reshape(-1, 1)))
    
    print(f"Full range (0-100): R² = {r2_full:.4f}")
    print(f"Restricted (50-60): R² = {r2_restricted:.4f}")
    print("→ Same model quality, but range affects R² dramatically!")
    
    # Pitfall 3: R² with non-proportional variance
    print("
=== Pitfall 3: Heteroscedasticity ===")
    x = np.linspace(1, 100, 100)
    # Noise increases with x (heteroscedastic)
    y_hetero = 2 * x + np.random.randn(100) * x * 0.5
    
    lr.fit(x.reshape(-1, 1), y_hetero)
    r2 = r2_score(y_hetero, lr.predict(x.reshape(-1, 1)))
    print(f"Heteroscedastic data R²: {r2:.4f}")
    print("→ High R², but residuals aren't random—model assumptions violated!")
 
demonstrate_pitfalls()