Linear Regression: Fitting a Line to Data

Linear regression: fitting a line to minimize error

Linear regression is the simplest supervised learning method: fit a line (or hyperplane in higher dimensions) to data such that the sum of squared residuals is minimized. A residual is the vertical distance from a data point to the fitted line — the error for that point.

The optimization problem: Given inputs X (n_samples × n_features) and targets y (n_samples × 1), find weights w and bias b that minimize:

J(w, b) = (1/n) · Σᵢ(y_i - (w·xᵢ + b))² = (1/n) · ||y - (Xw + b)||²

This is the mean squared error (MSE) or least-squares objective.

Solution 1: Closed-form (normal equation)

Taking the derivative with respect to w and setting it to zero:

∂J/∂w = -(2/n) · X^T(y - Xw) = 0

Solving for w:

w = (X^T X)^(-1) X^T y

This is the least-squares solution. It computes the best-fit line in one step without iteration. Trade-off: Fast for small n and small p (features), but computing the matrix inverse is O(n·p²), which becomes slow for millions of samples or thousands of features. Also, if X^T X is ill-conditioned (features are highly correlated), the inverse is numerically unstable.

Solution 2: Gradient descent (iterative)

Start with random w; repeatedly update:

w_new = w_old - learning_rate · ∂J/∂w = w_old - learning_rate · (2/n) · X^T(Xw_old - y)

Each step moves w a small distance in the direction that reduces the loss. After many iterations, w converges to the same optimum as the normal equation (if learning rate is small enough). Trade-off: Slower per-iteration but scales better to large datasets; doesn't require matrix inversion; can be stopped early if needed.

Both find the same solution; choose based on dataset size and constraints.

Normal equation vs gradient descent: a comparison

| Aspect | Normal Equation | Gradient Descent | |--------|-----------------|------------------| | Formula | w = (X^T X)^(-1) X^T y | w ← w - lr · ∂J/∂w, repeated | | Iterations needed | 1 (analytical) | Many (100s–1000s) | | Per-iteration cost | O(n·p²) matrix inversion | O(n·p) matrix multiply | | Total time | Fast for n < 10K, p < 100 | Better for n > 100K or p > 1000 | | Numerical stability | Can be unstable if X^T X is ill-conditioned | More stable; doesn't require inversion | | Learning rate tuning | N/A | Required; affects convergence speed | | Can stop early? | No | Yes; useful if convergence is fast | | Best use case | Small datasets, offline learning | Large datasets, online/streaming learning |

For the scikit-learn LinearRegression() on small datasets (< 100K samples), the solver uses the normal equation by default. For very large datasets, use SGDRegressor (stochastic gradient descent) instead.

Multicollinearity: when features are correlated

A subtle but important problem: if two features are highly correlated (e.g., square footage and number of bedrooms), the coefficients become unstable. Small changes in the data lead to large swings in the fitted weights.

Example: Predicting house price from square footage and bedrooms. If bedrooms ≈ square_footage / 300 (bedrooms is roughly square footage divided by 300), then many combinations of (w_sqft, w_beds) fit the data equally well. The normal equation produces one arbitrary solution.

How to detect: Compute the correlation matrix of features and look for values near ±1. Or compute the condition number of X^T X; if it's > 1000, multicollinearity is severe.

How to fix:

• Use regularization: L2 (ridge regression) or L1 (lasso) adds a penalty on large weights, forcing coefficients to stay small
• Drop redundant features: keep bedrooms or square footage, but not both
• Use principal component analysis (PCA) to create uncorrelated features

Assumptions and when they break

Linear regression assumes:

• Linearity: The relationship between inputs and output is actually linear (or close to it).
• Independence: Observations are independent (no clustering, no time series autocorrelation).
• Homoscedasticity: Error variance is constant across the range of the output (not heterogeneous).
• Normality: Errors are normally distributed (less critical for prediction, more for confidence intervals).

When these fail:

• Nonlinearity: Fit is poor; consider polynomial features, splines, or tree-based methods.
• Dependence: Confidence intervals are too narrow; use methods for correlated data.
• Heteroscedasticity: Predictions are less reliable in high-variance regions; use weighted regression or robust methods.
• Non-normal errors: Confidence intervals are unreliable, but predictions can still be good.

Always plot residuals (predicted error vs predicted value) to diagnose violations:

import matplotlib.pyplot as plt
import numpy as np

# After fitting a model, compute residuals
y_pred = model.predict(X)
residuals = y_test - y_pred

# Plot residuals vs predicted values
plt.figure(figsize=(12, 4))

plt.subplot(1, 2, 1)
plt.scatter(y_pred, residuals, alpha=0.5)
plt.axhline(y=0, color='r', linestyle='--')
plt.xlabel("Predicted value")
plt.ylabel("Residual (actual - predicted)")
plt.title("Residuals vs Predicted (should be random cloud around zero)")

# Q-Q plot: are residuals normal?
from scipy import stats
plt.subplot(1, 2, 2)
stats.probplot(residuals, dist="norm", plot=plt)
plt.title("Q-Q Plot (should follow diagonal line for normality)")

plt.tight_layout()
plt.show()

Full worked example: predicting house prices

Let's fit a linear regression model on a real dataset (Boston Housing or similar) and interpret the results.

import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
import matplotlib.pyplot as plt

# Load data: synthetic house price dataset
# Features: square footage, bedrooms, bathrooms, age, etc.
np.random.seed(42)
n_samples = 200
X = np.random.randn(n_samples, 4)
X[:, 0] = np.abs(X[:, 0]) * 50 + 1000  # square footage (1000–2500)
X[:, 1] = np.round(np.abs(X[:, 1]) * 1.5 + 2)  # bedrooms (1–5)
X[:, 2] = np.round(np.abs(X[:, 2]) + 1.5)  # bathrooms (1–3)
X[:, 3] = np.abs(X[:, 3]) * 30 + 20  # age (20–80 years)

# True linear relationship (with noise)
true_weights = np.array([0.15, 30000, 20000, -500])  # coefficient for each feature
true_intercept = 50000
y = X @ true_weights + true_intercept + np.random.randn(n_samples) * 50000

# Create DataFrame for interpretability
feature_names = ["sq_ft", "bedrooms", "bathrooms", "age_years"]
df = pd.DataFrame(X, columns=feature_names)
df["price"] = y

# Split into training and test
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

# Standardize features (important for interpretation and numerical stability)
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

# Fit linear regression
model = LinearRegression()
model.fit(X_train_scaled, y_train)

# Predictions
y_train_pred = model.predict(X_train_scaled)
y_test_pred = model.predict(X_test_scaled)

# Evaluate with multiple metrics
train_r2 = r2_score(y_train, y_train_pred)
test_r2 = r2_score(y_test, y_test_pred)
train_rmse = np.sqrt(mean_squared_error(y_train, y_train_pred))
test_rmse = np.sqrt(mean_squared_error(y_test, y_test_pred))
test_mae = mean_absolute_error(y_test, y_test_pred)

print("=" * 70)
print("LINEAR REGRESSION: HOUSE PRICE PREDICTION")
print("=" * 70)
print()
print("Model coefficients (standardized features):")
for name, coef in zip(feature_names, model.coef_):
    print(f"  {name:15s}: {coef:10.1f}")
print(f"  Intercept:       {model.intercept_:10.1f}")
print()
print("Interpretation (standardized coefficients):")
print(f"  Coefficient of {feature_names[0]}: {model.coef_[0]:.1f}")
print(f"    → A 1 std-dev increase in sq_ft (from scaling) → ${model.coef_[0]:,.0f} price change")
print()
print(f"Model Performance:")
print(f"  Train R² (coefficient of determination): {train_r2:.3f}")
print(f"    → Model explains {train_r2*100:.1f}% of variance in training prices")
print(f"  Test R²:                                  {test_r2:.3f}")
print(f"    → Generalization is {'good' if test_r2 > 0.8 else 'fair' if test_r2 > 0.6 else 'poor'}")
print()
print(f"  Train RMSE (root mean squared error):    ${train_rmse:,.0f}")
print(f"  Test RMSE:                                ${test_rmse:,.0f}")
print(f"    → On average, predictions are off by ±${test_rmse:,.0f}")
print()
print(f"  Test MAE (mean absolute error):           ${test_mae:,.0f}")
print(f"    → Median error is around ${test_mae:,.0f}")
print()

# Visualize fit
plt.figure(figsize=(14, 5))

plt.subplot(1, 2, 1)
plt.scatter(y_test, y_test_pred, alpha=0.5, s=50)
# Perfect fit line
min_val, max_val = y_test.min(), y_test.max()
plt.plot([min_val, max_val], [min_val, max_val], 'r--', lw=2, label='Perfect fit')
plt.xlabel("Actual price ($)")
plt.ylabel("Predicted price ($)")
plt.title("Predicted vs Actual (Test Set)")
plt.legend()
plt.grid(True, alpha=0.3)

plt.subplot(1, 2, 2)
residuals = y_test - y_test_pred
plt.scatter(y_test_pred, residuals, alpha=0.5, s=50)
plt.axhline(y=0, color='r', linestyle='--', lw=2)
plt.xlabel("Predicted price ($)")
plt.ylabel("Residual ($)")
plt.title("Residuals (should be random around zero)")
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

Interpreting R² and other metrics

R² (coefficient of determination): Fraction of variance in y explained by the model.

• R² = 1: Perfect fit
• R² = 0.7: Good fit (explains 70% of variance)
• R² = 0: Model is no better than predicting the mean
• R² < 0: Model is worse than predicting the mean (overfitting or poor generalization)

RMSE (root mean squared error): √(mean of squared errors). Penalizes large errors more heavily.

• Interpretation: on average, predictions are off by ±RMSE

MAE (mean absolute error): Mean of absolute errors. More robust to outliers than RMSE.

• Interpretation: median error magnitude

Residuals: y_actual - y_pred. Should be randomly scattered around zero with constant variance (homoscedastic). If residuals have a pattern (e.g., they increase with predicted value), the linear model is missing structure.

Multiple linear regression and feature interaction

Linear regression extends to multiple features. Each feature gets a coefficient; the model sums their weighted contributions:

ŷ = w₀ + w₁·x₁ + w₂·x₂ + ... + wₙ·xₙ

Interpretation of coefficients: Holding all other features fixed, a unit increase in feature j is associated with a w_j increase in the output (for non-standardized features, in the original units; for standardized features, a 1 std-dev increase).

Feature interaction: If the effect of x₁ depends on x₂, linear regression misses this. You would need to add an interaction term (x₁ × x₂) as a new feature:

# Add interaction term manually
X_with_interaction = X.copy()
X_with_interaction = np.column_stack([X, X[:, 0] * X[:, 1]])  # sq_ft × bedrooms interaction

# Re-fit model with interaction
model_interaction = LinearRegression()
model_interaction.fit(X_with_interaction, y)

# Check if interaction improves R²
y_pred_interaction = model_interaction.predict(X_with_interaction)
r2_with = r2_score(y, y_pred_interaction)
print(f"R² without interaction: {r2_score(y, model.predict(X)):.3f}")
print(f"R² with interaction:    {r2_with:.3f}")

When linear regression succeeds and fails

Linear regression works well when:

• The relationship is actually linear or nearly so
• You have more samples than features (n >> p)
• Features are not highly correlated (multicollinearity is low)
• You need interpretability (coefficients tell a clear story)

Linear regression struggles when:

• The relationship is nonlinear (e.g., exponential, sinusoidal)
• Features are highly correlated (unstable coefficients)
• You have outliers (they dominate the least-squares fit)
• You have far more features than samples (overfitting, ridge/lasso regression needed)

For nonlinear relationships, consider polynomial regression (add squared/cubed features), splines, or tree-based methods (random forests, gradient boosting). For high-dimensional data, use regularization (ridge or lasso) to stabilize the fit.

Regularization: Ridge and Lasso regression for high-dimensional data

When n (samples) is much smaller than p (features), or features are correlated, vanilla linear regression overfits. Regularization adds a penalty on large weights to prevent overfitting.

Ridge regression (L2 regularization):

minimize: J(w) = (1/n) · ||y - Xw||² + λ · ||w||²

The λ term (regularization strength) penalizes large weights. The closed-form solution is:

w = (X^T X + λI)^(-1) X^T y

This is similar to the normal equation but with λI added to the diagonal, making the matrix better-conditioned (more numerically stable).

Lasso regression (L1 regularization):

minimize: J(w) = (1/n) · ||y - Xw||² + λ · ||w||

This penalizes the sum of absolute values instead of squares. A unique property: some coefficients go exactly to zero, performing automatic feature selection.

from sklearn.linear_model import Ridge, Lasso
import numpy as np

# Ridge regression
ridge = Ridge(alpha=1.0)  # alpha is the regularization strength (λ)
ridge.fit(X_train, y_train)
ridge_score = ridge.score(X_test, y_test)
print(f"Ridge R²: {ridge_score:.3f}")

# Lasso regression
lasso = Lasso(alpha=0.1, max_iter=1000)
lasso.fit(X_train, y_train)
lasso_score = lasso.score(X_test, y_test)
n_zero_coef = (lasso.coef_ == 0).sum()
print(f"Lasso R²: {lasso_score:.3f}")
print(f"  Coefficients set to zero: {n_zero_coef} / {len(lasso.coef_)}")

# Elastic Net: hybrid of Ridge and Lasso
from sklearn.linear_model import ElasticNet
elastic = ElasticNet(alpha=0.1, l1_ratio=0.5, max_iter=1000)
elastic.fit(X_train, y_train)
print(f"Elastic Net R²: {elastic.score(X_test, y_test):.3f}")

| Aspect | Ridge | Lasso | Elastic Net | |--------|-------|-------|---------| | Penalty term | λ · Σw² | λ · Σ|w| | λ · [ρ · Σ|w| + (1-ρ) · Σw²] | | Shrinks coefficients | Continuously toward zero | Some to exactly zero | Some to zero, others shrunk | | Feature selection | No (all features retained) | Yes (automatic) | Partial (depends on ρ) | | When to use | Correlated features (keep all signal) | High-dimensional (p >> n); want interpretability | Balanced; get some selection benefits | | Typical α range | 0.1 – 10 | 0.001 – 0.1 | 0.001 – 1 |

Common mistake: confusing causation with correlation

A strong coefficient does not imply causation. If house prices correlate with zip code, the zip code does not cause the price; location, school quality, and neighborhood desirability do. Always think about whether a feature is a confounder (a third variable that drives both the feature and the target) or a causal mediator (the feature itself affects the target).

Linear regression estimates correlation, not causation. Causal inference requires experimentation or causal reasoning, not just data fitting.

Real-world benchmark: predicting house prices

The California Housing dataset (20,640 houses, 8 features: latitude, longitude, housing age, total rooms, total bedrooms, population, households, median income) is a classic benchmark.

Typical results (illustrative estimates from published benchmarks):

| Model | RMSE (test) | R² (test) | Interpretation | |-------|---------|---------|---------| | Baseline (mean price) | $73,500 | 0.0 | No signal captured | | Linear regression | $73,200 | 0.58 | Explains 58% of variance | | Linear with interactions | $71,800 | 0.62 | +4% gain from feature engineering | | Ridge (α=10) | $72,500 | 0.60 | Slightly worse than OLS; regularization too strong | | Lasso (α=0.01) | $72,900 | 0.59 | Similar to Ridge; some features zeroed out | | Polynomial features (degree 2) | $68,900 | 0.71 | +13% gain from nonlinear features | | Decision tree (max_depth=10) | $71,100 | 0.65 | Competitive; captures nonlinearity |

Interpretation: Linear regression is a strong baseline, explaining 58% of price variance from just 8 raw features. Adding interactions and polynomial features boosts this to 71%, but at the cost of interpretability. A decision tree matches this performance while being more interpretable (you can trace which features matter in each split).

When linear regression fails: signs and solutions

Problem 1: Nonlinear relationship

• Sign: Residual plot shows a clear pattern (U-shape, S-curve, etc.)
• Solution: Add polynomial features (x², x³), use splines, or switch to nonlinear models (trees, neural networks)

Problem 2: Heteroscedasticity (non-constant variance)

• Sign: Residuals are small for some predicted values, large for others (e.g., scatter fan out)
• Solution: Use weighted least squares, transform the target (log, square root), or use robust regression

Problem 3: Outliers

• Sign: A few extreme residuals; coefficients change dramatically if you remove them
• Solution: Use robust regression (Huber loss), cap outliers, or investigate why they exist

Problem 4: Multicollinearity

• Sign: Coefficients have large standard errors; small data changes cause large coefficient swings
• Solution: Use Ridge/Lasso, drop redundant features, or use PCA

Problem 5: Autocorrelation (in time series)

• Sign: Residuals are correlated over time; confidence intervals are too narrow
• Solution: Use time-series models (ARIMA), add lagged features, or use GARCH for heteroscedasticity

All of these can be diagnosed visually by plotting residuals and predicted values, or statistically via residual tests.