Multiple Regression

Multiple linear regression extends simple linear regression to use two or more independent variables to predict a continuous target. Rather than drawing a line through two dimensions, the model fits a hyperplane through as many dimensions as there are features. Understanding how the coefficients interact — and when they mislead you — is the core skill this page teaches.

This page covers:

• The multiple regression equation and what each coefficient means
• How to build a complete scikit-learn pipeline: load, preprocess, train, evaluate
• Why feature scaling matters and how to do it correctly
• How to interpret and compare scaled vs unscaled coefficients
• Diagnosing multicollinearity — the most common gotcha in multiple regression
• Residual analysis to check model assumptions
• When to choose multiple regression and what to try when it falls short

The Multiple Regression Equation

Multiple linear regression models the target y as a linear combination of n input features:

y = β₀ + β₁x₁ + β₂x₂ + ... + βₙxₙ + ε

• β₀ — the intercept: the predicted value of y when every feature equals zero
• β₁ … βₙ — the coefficients: how much y changes for a one-unit increase in each xᵢ, holding all other features constant
• ε — the error term: the part of y that the model cannot explain

The algorithm finds the coefficients by minimizing the sum of squared residuals (ordinary least squares):

SSR = Σ(yᵢ - ŷᵢ)²

This has an exact closed-form solution, so scikit-learn's LinearRegression does not need iterative gradient descent — training is nearly instant even on datasets with hundreds of thousands of rows.

Difference from Simple Linear Regression

Simple linear regression uses one feature. Multiple regression adds more features so that each coefficient captures the partial effect of that feature — its impact on the target while the remaining features are held constant. This is more powerful but introduces new risks, especially multicollinearity (see Diagnosing Multicollinearity).

The Dataset

The examples below use the California Housing dataset built into scikit-learn. It records census-block-level housing statistics for California in 1990 and has 20,640 samples across 8 features.

import pandas as pd
from sklearn.datasets import fetch_california_housing

housing = fetch_california_housing()
df = pd.DataFrame(housing.data, columns=housing.feature_names)
df['MedHouseVal'] = housing.target  # median house value in $100,000s

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

The 8 input features are:

The target MedHouseVal is the median house value in units of $100,000, so a value of 2.0 represents $200,000.

Building the Model Step by Step

Step 1 — Split the Data

Always split data before any preprocessing. Fitting a scaler on the full dataset would leak test-set statistics into training, giving you an overly optimistic evaluation. See Train/Test Split for a full explanation.

from sklearn.model_selection import train_test_split

X = df[housing.feature_names]
y = df['MedHouseVal']

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

print(f"Training samples: {len(X_train)}")  # 16512
print(f"Test samples:     {len(X_test)}")   # 4128

Step 2 — Scale the Features

Multiple regression coefficients reflect the units of each feature. MedInc is measured in tens of thousands of dollars; Population is a raw count that can reach 35,000. Without scaling, the coefficient for Population will be tiny not because population is unimportant, but because its unit is small.

StandardScaler transforms each feature to have zero mean and unit standard deviation, making coefficient magnitudes directly comparable. See Feature Scaling for details on the available scalers.

from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)   # fit ONLY on training data
X_test_scaled  = scaler.transform(X_test)         # apply the same transformation

The critical rule: call fit_transform on training data and transform (no fitting) on test data. Fitting on test data would contaminate the evaluation.

Step 3 — Train the Model

from sklearn.linear_model import LinearRegression

model = LinearRegression()
model.fit(X_train_scaled, y_train)

LinearRegression.fit() solves the OLS problem analytically. There are no hyperparameters to tune for vanilla multiple regression — what you choose is which features to include.

Step 4 — Evaluate the Model

import numpy as np
from sklearn.metrics import mean_squared_error, r2_score

y_pred = model.predict(X_test_scaled)

mse  = mean_squared_error(y_test, y_pred)
rmse = np.sqrt(mse)
r2   = r2_score(y_test, y_pred)

print(f"Mean Squared Error:      {mse:.4f}")
print(f"Root Mean Squared Error: {rmse:.4f}  (±${rmse * 100_000:,.0f})")
print(f"R-squared:               {r2:.4f}")

Expected output:

Mean Squared Error:      0.5559
Root Mean Squared Error: 0.7456  (±$74,558)
R-squared:               0.5758

What the metrics mean:

• MSE (Mean Squared Error) — average of squared differences between predictions and actual values. Squaring penalizes large errors more heavily than small ones. Units are the square of the target ($100,000s²), so it is harder to interpret directly.
• RMSE (Root Mean Squared Error) — the square root of MSE, back in the same units as the target. An RMSE of 0.75 means the model's predictions are off by about $75,000 on average.
• R² (coefficient of determination) — the fraction of variance in the target that the model explains. An R² of 0.58 means the model captures 58% of the variation in house prices. Values range from 0 (no better than predicting the mean) to 1 (perfect predictions). Negative R² is possible if the model is worse than the mean — this is a strong sign something is wrong.

An R² of ~0.58 is typical for this dataset with vanilla linear regression. The relationship between house prices and these features is partly non-linear and involves geographic clustering that a hyperplane cannot capture well. Algorithms like gradient-boosted trees routinely reach 0.80+ on this dataset.

Step 5 — Inspect the Coefficients

After scaling, coefficient magnitudes are directly comparable — they show which features most strongly influence the prediction:

coef_df = pd.DataFrame({
    'Feature':     housing.feature_names,
    'Coefficient': model.coef_
}).sort_values('Coefficient', key=abs, ascending=False)

print(coef_df.to_string(index=False))
print(f"\nIntercept: {model.intercept_:.4f}")

Expected output:

   Feature  Coefficient
  Latitude    -0.8969
 Longitude    -0.8698
    MedInc     0.8544
 AveBedrms     0.3393
  AveRooms    -0.2944
  HouseAge     0.1225
  AveOccup    -0.0408
Population    -0.0023

Intercept: 2.0719

Reading the coefficients after standard scaling:

• MedInc = 0.854 — the single strongest predictor. A one-standard-deviation increase in median income predicts a $85,400 increase in house value, holding everything else constant.
• Latitude = -0.897 and Longitude = -0.869 — the model has learned that more northern and more eastern census blocks tend to be cheaper. However, these two geographic features are highly correlated (r = -0.93), which can make their individual coefficients unstable (see Diagnosing Multicollinearity).
• AveBedrms = +0.339 vs AveRooms = -0.294 — these have opposite signs even though more rooms and more bedrooms generally mean bigger, more expensive houses. This is a classic sign of multicollinearity: AveRooms and AveBedrms are correlated (r = 0.85), so their coefficients compensate for each other. Do not interpret them in isolation.
• Population = -0.002 — very close to zero after scaling. Block population has minimal predictive power once the other features are accounted for.

The intercept (2.07) is the predicted MedHouseVal when every scaled feature is zero — that is, when all features are at their training-set mean. It equals the mean of the training targets and has no direct business meaning beyond that.

Step 6 — Make Predictions on New Data

# A new census block: high income, older house, San Francisco Bay Area
new_block = pd.DataFrame([[8.0, 41.0, 6.0, 1.0, 322, 2.5, 37.88, -122.23]],
                         columns=housing.feature_names)

new_block_scaled = scaler.transform(new_block)
prediction = model.predict(new_block_scaled)

print(f"Predicted median house value: ${prediction[0] * 100_000:,.0f}")

Expected output:

Predicted median house value: $410,895

The scaler must be the same scaler fitted on training data. Never refit the scaler on the new data — that would shift the inputs relative to what the model learned.

Diagnosing Multicollinearity

Multicollinearity occurs when two or more independent variables are highly correlated with each other. It does not prevent the model from making accurate predictions, but it makes the individual coefficients unreliable and hard to interpret. Coefficients can become large, change sign, or become statistically insignificant even for genuinely important features.

Checking with a Correlation Matrix

import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

corr = df[housing.feature_names].corr()
print(corr.round(2))

Key pairs to notice:

A correlation above 0.80 between two features is a warning sign. When you see it, consider:

1. Drop one of the correlated features. If AveRooms and AveBedrms are both in the model, try removing AveBedrms and check whether the model's predictive performance changes much.
2. Combine them. Create a derived feature (e.g., rooms_per_bedroom = AveRooms / AveBedrms) that captures the relationship without redundancy.
3. Use a regularized model. Ridge regression adds an L2 penalty that shrinks correlated coefficients toward each other, stabilizing them. Lasso (L1) can zero out redundant features entirely.

Residual Analysis

A residual is the difference between an actual value and the model's prediction: residual = y_actual - y_predicted. Plotting residuals reveals whether the model's assumptions hold.

residuals = y_test - y_pred

print(f"Mean of residuals: {residuals.mean():.4f}")  # should be close to 0
print(f"Std of residuals:  {residuals.std():.4f}")

Expected output:

Mean of residuals: 0.0035
Std of residuals:  0.7457

The mean is near zero — a sign that the model is unbiased on average. But the standard deviation of 0.75 (±$75,000) shows substantial scatter.

import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

# Predicted vs Actual
axes[0].scatter(y_test, y_pred, alpha=0.2, s=8)
axes[0].plot([0, 5], [0, 5], 'r--')
axes[0].set_xlabel('Actual')
axes[0].set_ylabel('Predicted')
axes[0].set_title('Predicted vs Actual')

# Residuals vs Predicted
axes[1].scatter(y_pred, residuals, alpha=0.2, s=8)
axes[1].axhline(0, color='r', linestyle='--')
axes[1].set_xlabel('Predicted')
axes[1].set_ylabel('Residual')
axes[1].set_title('Residuals vs Predicted')

plt.tight_layout()
plt.savefig('residuals.png', dpi=120)
print("Saved residuals.png")

What to look for:

• Predicted vs Actual — ideally, points fall along the diagonal. Systematic deviation from the diagonal (a curve, or a flattening at high values) means the linear assumption is wrong for part of the range.
• Residuals vs Predicted — ideally, residuals scatter randomly around zero at all prediction levels. A funnel shape (wider spread at higher predictions) signals heteroscedasticity — the model's error is not constant, which can make interval estimates untrustworthy.
• Clusters — distinct groups in the residual plot can reveal that the dataset contains sub-populations (e.g., urban vs rural) that warrant separate models or additional features.

Full Pipeline

Here is everything above as a single runnable script:

import numpy as np
import pandas as pd
from sklearn.datasets import fetch_california_housing
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import mean_squared_error, r2_score

# 1. Load data
housing = fetch_california_housing()
df = pd.DataFrame(housing.data, columns=housing.feature_names)
df['MedHouseVal'] = housing.target

# 2. Split — before any preprocessing
X = df[housing.feature_names]
y = df['MedHouseVal']
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

# 3. Scale
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled  = scaler.transform(X_test)

# 4. Train
model = LinearRegression()
model.fit(X_train_scaled, y_train)

# 5. Evaluate
y_pred = model.predict(X_test_scaled)
rmse = np.sqrt(mean_squared_error(y_test, y_pred))
r2   = r2_score(y_test, y_pred)
print(f"RMSE: {rmse:.4f}")
print(f"R²:   {r2:.4f}")

# 6. Coefficients (most important first)
coef_df = pd.DataFrame({
    'Feature':     housing.feature_names,
    'Coefficient': model.coef_
}).sort_values('Coefficient', key=abs, ascending=False)
print(coef_df.to_string(index=False))

When to Use Multiple Regression

Multiple regression is a strong first choice when:

• The relationship between each feature and the target is approximately linear
• Interpretability matters — each coefficient has a clear meaning
• You want a fast baseline before trying complex models
• You have enough samples relative to the number of features (a rough rule: at least 10–20 observations per feature)

Consider alternatives when:

Common Pitfalls

Fitting the scaler on all data before splitting. This leaks test-set mean and variance into training. Always split first, then fit the scaler only on the training portion.

Adding more features always helps. Adding irrelevant or redundant features can reduce interpretability, introduce multicollinearity, and hurt generalization. Use domain knowledge or a feature selection technique to choose features deliberately.

Trusting coefficients when features are correlated. When AveRooms and AveBedrms are both in the model, neither coefficient reliably reflects the true effect of that variable. Check the correlation matrix before interpreting individual coefficients.

Ignoring residual plots. An R² of 0.58 looks reasonable on paper, but residual patterns can reveal that the model is systematically wrong for high-value properties or for specific geographic regions.

Extrapolating beyond the training range. A linear model trained on houses priced between $50,000 and $500,000 should not be used to predict $5,000,000 properties. Check that new inputs fall within the range seen during training.

Next Steps

• Linear Regression — the foundation: simple (one-feature) regression and the OLS method explained in full
• Polynomial Regression — extend the linear model to capture curves by adding polynomial feature terms
• Feature Scaling — deep dive into StandardScaler, MinMaxScaler, and RobustScaler
• Train/Test Split — why proper splitting is critical and how to do it without leaking data
• Cross-Validation — a more robust alternative to a single train/test split