Multiple Regression
Model relationships between multiple input features and a target variable, incorporating interactions and feature engineering.
Concept Overview
Multiple regression extends linear regression to handle multiple predictor variables simultaneously. It allows us to model complex relationships where the target depends on several features and their interactions.
The multiple regression model with p features is:
ŷ = β₀ + β₁x₁ + β₂x₂ + … + βₚxₚ + ε
Where:
- ŷ (y-hat) is the predicted target value
- β₀ is the intercept (y-intercept when all features are zero)
- β₁, β₂, …, βₚ are the coefficients for each feature
- x₁, x₂, …, xₚ are the input features
- ε (epsilon) is the error term
Key Advantages:
- Can model relationships with multiple predictors
- Allows for control of confounding variables
- Enables testing of interaction effects
- More powerful than simple regression when multiple factors matter
Basic Multiple Regression
Let's start with a basic example using multiple features to predict a target variable.
import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score, mean_absolute_error
from sklearn.model_selection import train_test_split
# Generate sample data with multiple features
np.random.seed(42)
n = 200
sqft = np.random.normal(1500, 500, n)
bedrooms = np.random.randint(1, 6, n)
age = np.random.normal(10, 5, n)
bathrooms = bedrooms * 0.7 + np.random.normal(0, 0.3, n)
# Target with multiple features
price = (100 * sqft + 20000 * bedrooms + 15000 * bathrooms - 5000 * age +
np.random.normal(0, 50000, n))
df = pd.DataFrame({
'sqft': sqft,
'bedrooms': bedrooms,
'bathrooms': bathrooms,
'age': age,
'price': price
})
# Prepare features and target
X = df[['sqft', 'bedrooms', 'bathrooms', 'age']]
y = df['price']
# Split data
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
# Create and train model
model = LinearRegression()
model.fit(X_train, y_train)
# Make predictions
y_train_pred = model.predict(X_train)
y_test_pred = model.predict(X_test)
# Model parameters
print("Multiple Regression Model:")
print(f"Intercept (β₀): {model.intercept_:.2f}")
print("\nCoefficients:")
for i, feature in enumerate(X.columns):
print(f" {feature} (β{i+1}): {model.coef_[i]:.2f}")
print(f"\nModel equation:")
print(f"price = {model.intercept_:.2f}", end="")
for i, feature in enumerate(X.columns):
sign = "+" if model.coef_[i] >= 0 else ""
print(f" {sign}{model.coef_[i]:.2f} × {feature}", end="")
# Evaluate metrics
train_r2 = r2_score(y_train, y_train_pred)
test_r2 = r2_score(y_test, y_test_pred)
train_mse = mean_squared_error(y_train, y_train_pred)
test_mse = mean_squared_error(y_test, y_test_pred)
train_mae = mean_absolute_error(y_train, y_train_pred)
test_mae = mean_absolute_error(y_test, y_test_pred)
print(f"\n\nTraining Metrics:")
print(f" R²: {train_r2:.4f}")
print(f" MSE: {train_mse:.2f}")
print(f" MAE: {train_mae:.2f}")
print(f"\nTest Metrics:")
print(f" R²: {test_r2:.4f}")
print(f" MSE: {test_mse:.2f}")
print(f" MAE: {test_mae:.2f}")
# Feature importance (coefficient magnitude, standardized)
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)
model_scaled = LinearRegression()
model_scaled.fit(X_train_scaled, y_train)
print(f"\n\nStandardized Coefficients (for comparison):")
for i, feature in enumerate(X.columns):
print(f" {feature}: {model_scaled.coef_[i]:.2f}")Standardized coefficients allow comparison of feature importance when features have different scales. Larger absolute values indicate greater impact on the target variable.
Interaction Terms
Interaction terms allow features to influence each other's effect on the target. For example, the effect of square footage might depend on the number of bedrooms.
import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from sklearn.metrics import r2_score
from sklearn.model_selection import train_test_split
# Generate data with interaction effects
np.random.seed(42)
n = 200
sqft = np.random.normal(1500, 500, n)
bedrooms = np.random.randint(1, 6, n)
# Interaction: price increases more with sqft when there are more bedrooms
price = (100 * sqft + 20000 * bedrooms + 5 * sqft * bedrooms +
np.random.normal(0, 50000, n))
df = pd.DataFrame({
'sqft': sqft,
'bedrooms': bedrooms,
'price': price
})
X = df[['sqft', 'bedrooms']]
y = df['price']
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
# Model without interactions
model_no_interaction = LinearRegression()
model_no_interaction.fit(X_train, y_train)
pred_no_int = model_no_interaction.predict(X_test)
r2_no_int = r2_score(y_test, pred_no_int)
print("Model without interactions:")
print(f"Test R²: {r2_no_int:.4f}")
print(f"Coefficients: sqft={model_no_interaction.coef_[0]:.2f}, bedrooms={model_no_interaction.coef_[1]:.2f}")
# Model with interactions (degree 2 polynomial includes interactions)
from sklearn.preprocessing import PolynomialFeatures
poly_features = PolynomialFeatures(degree=2, include_bias=False, interaction_only=True)
X_train_inter = poly_features.fit_transform(X_train)
X_test_inter = poly_features.transform(X_test)
feature_names = poly_features.get_feature_names_out(X.columns)
print(f"\nInteraction features created: {feature_names}")
model_interaction = LinearRegression()
model_interaction.fit(X_train_inter, y_train)
pred_inter = model_interaction.predict(X_test_inter)
r2_inter = r2_score(y_test, pred_inter)
print(f"\nModel with interactions:")
print(f"Test R²: {r2_inter:.4f}")
print("Coefficients:")
for i, name in enumerate(feature_names):
print(f" {name}: {model_interaction.coef_[i]:.2f}")
# Manual interaction term
X_train_manual = X_train.copy()
X_train_manual['sqft_x_bedrooms'] = X_train_manual['sqft'] * X_train_manual['bedrooms']
X_test_manual = X_test.copy()
X_test_manual['sqft_x_bedrooms'] = X_test_manual['sqft'] * X_test_manual['bedrooms']
model_manual = LinearRegression()
model_manual.fit(X_train_manual, y_train)
pred_manual = model_manual.predict(X_test_manual)
r2_manual = r2_score(y_test, pred_manual)
print(f"\nModel with manual interaction term:")
print(f"Test R²: {r2_manual:.4f}")
print("Coefficients:")
for i, feature in enumerate(X_train_manual.columns):
print(f" {feature}: {model_manual.coef_[i]:.2f}")
print(f"\n\nImprovement from interactions:")
print(f"Without interactions R²: {r2_no_int:.4f}")
print(f"With interactions R²: {r2_inter:.4f}")
print(f"Improvement: {(r2_inter - r2_no_int):.4f}")Interaction terms can significantly improve model performance when features interact. Test for interactions by comparing models with and without interaction terms, or use domain knowledge to identify likely interactions.
Feature Selection
Not all features are equally important. Feature selection helps identify the most relevant features and can improve model performance and interpretability.
import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression, Ridge, Lasso
from sklearn.feature_selection import SelectKBest, f_regression, RFE
from sklearn.metrics import r2_score, mean_squared_error
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.preprocessing import StandardScaler
# Generate data with some irrelevant features
np.random.seed(42)
n = 200
sqft = np.random.normal(1500, 500, n)
bedrooms = np.random.randint(1, 6, n)
age = np.random.normal(10, 5, n)
noise1 = np.random.normal(0, 1, n) # Irrelevant
noise2 = np.random.normal(0, 1, n) # Irrelevant
noise3 = np.random.normal(0, 1, n) # Irrelevant
price = (100 * sqft + 20000 * bedrooms - 5000 * age +
np.random.normal(0, 50000, n))
df = pd.DataFrame({
'sqft': sqft,
'bedrooms': bedrooms,
'age': age,
'noise1': noise1,
'noise2': noise2,
'noise3': noise3,
'price': price
})
X = df.drop('price', axis=1)
y = df['price']
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
# 1. Univariate feature selection (F-statistic)
selector_f = SelectKBest(score_func=f_regression, k=3)
X_train_selected = selector_f.fit_transform(X_train, y_train)
X_test_selected = selector_f.transform(X_test)
selected_features = X.columns[selector_f.get_support()]
print(f"Selected features (F-test): {list(selected_features)}")
print(f"Feature scores: {dict(zip(X.columns, selector_f.scores_))}")
model_f = LinearRegression()
model_f.fit(X_train_selected, y_train)
pred_f = model_f.predict(X_test_selected)
r2_f = r2_score(y_test, pred_f)
print(f"\nModel with F-test selection - R²: {r2_f:.4f}")
# 2. Recursive Feature Elimination (RFE)
estimator = LinearRegression()
selector_rfe = RFE(estimator, n_features_to_select=3, step=1)
selector_rfe.fit(X_train, y_train)
selected_features_rfe = X.columns[selector_rfe.get_support()]
print(f"\nSelected features (RFE): {list(selected_features_rfe)}")
print(f"Feature rankings: {dict(zip(X.columns, selector_rfe.ranking_))}")
X_train_rfe = selector_rfe.transform(X_train)
X_test_rfe = selector_rfe.transform(X_test)
model_rfe = LinearRegression()
model_rfe.fit(X_train_rfe, y_train)
pred_rfe = model_rfe.predict(X_test_rfe)
r2_rfe = r2_score(y_test, pred_rfe)
print(f"Model with RFE - R²: {r2_rfe:.4f}")
# 3. Lasso for automatic feature selection
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)
lasso = Lasso(alpha=0.1, max_iter=2000)
lasso.fit(X_train_scaled, y_train)
selected_features_lasso = X.columns[np.abs(lasso.coef_) > 1e-5]
print(f"\nSelected features (Lasso): {list(selected_features_lasso)}")
print("Lasso coefficients:")
for i, feature in enumerate(X.columns):
print(f" {feature}: {lasso.coef_[i]:.4f}")
pred_lasso = lasso.predict(X_test_scaled)
r2_lasso = r2_score(y_test, pred_lasso)
print(f"Model with Lasso - R²: {r2_lasso:.4f}")
# Compare all models
print(f"\n\nModel Comparison:")
print(f"All features: R² = {r2_score(y_test, LinearRegression().fit(X_train, y_train).predict(X_test)):.4f}")
print(f"F-test selection: R² = {r2_f:.4f}")
print(f"RFE selection: R² = {r2_rfe:.4f}")
print(f"Lasso selection: R² = {r2_lasso:.4f}")Feature selection helps reduce overfitting, improve interpretability, and sometimes improve performance. Different methods work better in different scenarios—experiment with multiple approaches.
Model Diagnostics
Diagnosing multiple regression models involves checking assumptions, detecting multicollinearity, and evaluating model performance through various metrics.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score, mean_squared_error
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
# Generate data
np.random.seed(42)
n = 200
sqft = np.random.normal(1500, 500, n)
bedrooms = np.random.randint(1, 6, n)
age = np.random.normal(10, 5, n)
price = (100 * sqft + 20000 * bedrooms - 5000 * age +
np.random.normal(0, 50000, n))
df = pd.DataFrame({
'sqft': sqft,
'bedrooms': bedrooms,
'age': age,
'price': price
})
X = df[['sqft', 'bedrooms', 'age']]
y = df['price']
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
model = LinearRegression()
model.fit(X_train, y_train)
y_train_pred = model.predict(X_train)
y_test_pred = model.predict(X_test)
# Residuals
train_residuals = y_train - y_train_pred
test_residuals = y_test - y_test_pred
# Adjusted R²
n = len(y_train)
p = X_train.shape[1]
r2 = r2_score(y_train, y_train_pred)
adj_r2 = 1 - (1 - r2) * (n - 1) / (n - p - 1)
print("Model Diagnostics:")
print(f"R²: {r2:.4f}")
print(f"Adjusted R²: {adj_r2:.4f}")
print(f"Number of features (p): {p}")
print(f"Number of samples (n): {n}")
# Correlation matrix (check for multicollinearity)
correlation_matrix = X_train.corr()
print(f"\n\nCorrelation Matrix:")
print(correlation_matrix)
# VIF (Variance Inflation Factor) - requires statsmodels
try:
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor
X_with_const = sm.add_constant(X_train)
vif_data = pd.DataFrame()
vif_data["Feature"] = X_train.columns
vif_data["VIF"] = [variance_inflation_factor(X_with_const.values, i+1)
for i in range(len(X_train.columns))]
print(f"\n\nVariance Inflation Factor (VIF):")
print(vif_data)
print("\nInterpretation:")
print(" VIF < 5: Low multicollinearity (acceptable)")
print(" 5 ≤ VIF < 10: Moderate multicollinearity (caution)")
print(" VIF ≥ 10: High multicollinearity (problematic)")
except ImportError:
print("\nNote: Install statsmodels for VIF calculation")
# Residual plots
fig, axes = plt.subplots(2, 2, figsize=(14, 12))
# Residuals vs predicted
axes[0, 0].scatter(y_train_pred, train_residuals, alpha=0.6, s=50)
axes[0, 0].axhline(y=0, color='red', linestyle='--', linewidth=2)
axes[0, 0].set_xlabel('Predicted Values', fontsize=12)
axes[0, 0].set_ylabel('Residuals', fontsize=12)
axes[0, 0].set_title('Residuals vs Predicted', fontsize=14)
axes[0, 0].grid(True, alpha=0.3)
# Q-Q plot
from scipy import stats
stats.probplot(train_residuals, dist="norm", plot=axes[0, 1])
axes[0, 1].set_title('Q-Q Plot (Normality Check)', fontsize=14)
axes[0, 1].grid(True, alpha=0.3)
# Residuals histogram
axes[1, 0].hist(train_residuals, bins=20, edgecolor='black', alpha=0.7)
axes[1, 0].axvline(x=0, color='red', linestyle='--', linewidth=2)
axes[1, 0].set_xlabel('Residuals', fontsize=12)
axes[1, 0].set_ylabel('Frequency', fontsize=12)
axes[1, 0].set_title('Residuals Distribution', fontsize=14)
axes[1, 0].grid(True, alpha=0.3)
# Actual vs Predicted
axes[1, 1].scatter(y_train, y_train_pred, alpha=0.6, s=50)
min_val = min(y_train.min(), y_train_pred.min())
max_val = max(y_train.max(), y_train_pred.max())
axes[1, 1].plot([min_val, max_val], [min_val, max_val], 'r--', linewidth=2)
axes[1, 1].set_xlabel('Actual Values', fontsize=12)
axes[1, 1].set_ylabel('Predicted Values', fontsize=12)
axes[1, 1].set_title('Actual vs Predicted', fontsize=14)
axes[1, 1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()Always perform diagnostic checks on multiple regression models. Check residuals for patterns, examine multicollinearity with VIF, and use adjusted R² when comparing models with different numbers of features.
Best Practices
- Start with exploratory analysis: Understand relationships between features before modeling
- Check for multicollinearity: Use correlation matrices and VIF scores to detect highly correlated features
- Standardize when needed: Standardize features when comparing coefficient magnitudes or using regularization
- Consider interactions: Test for interaction effects when domain knowledge suggests they exist
- Feature selection: Remove irrelevant features to improve model performance and interpretability
- Use adjusted R²: Compare models with different numbers of features using adjusted R²
- Validate assumptions: Check linearity, homoscedasticity, normality, and independence
- Cross-validation: Use cross-validation to assess model performance robustly
- Regularization: Use Ridge or Lasso when features are numerous or multicollinear
- Document choices: Record feature selection and model building decisions, share insights through
info.studygrid@gmail.com
Next Steps
Learn about feature scaling to understand how to normalize features for better model performance and coefficient interpretability.
Try It Yourself
Exercise 1: Fit a model predicting CO2 from car weight and volume, then predict for weight 2300, volume 1300.
Show solution
from sklearn.linear_model import LinearRegression
X = [[790, 1000], [1160, 1200], [929, 1000], [865, 900], [1140, 1500]]
y = [99, 95, 95, 90, 105]
model = LinearRegression().fit(X, y)
print(model.predict([[2300, 1300]]).round(1))Exercise 2: Why should features be on similar scales before interpreting coefficients?
Show solution
Raw coefficients depend on each feature's units, so they are not directly comparable. Scaling (standardizing) lets you compare each feature's relative importance.
Key Takeaways
- Multiple regression uses several features to predict one target.
- Each coefficient is the effect of one feature holding others constant.
- Watch for multicollinearity between correlated features.
📘 Real-World Deep Dive
Knowing <strong>ML Multiple Regression (scikit-learn)</strong> well is what turns scikit-learn from a curiosity into a daily tool — you'll reach for it in nearly every real project.
Real-Life Scenario
An end-to-end usage of ML Multiple Regression that you'd actually see in a data pipeline or analytics notebook.
Real-Life Example
import numpy as np
from sklearn.linear_model import LinearRegression
X = np.array([[1, 2], [2, 1], [3, 4], [4, 3], [5, 5], [6, 7]])
y = np.array([5, 4, 11, 10, 15, 19])
m = LinearRegression().fit(X, y)
print("coef:", m.coef_, "intercept:", round(m.intercept_, 2))Expected Output
(see source)Common mistakes
fitexpects numeric arrays;OneHotEncoder/LabelEncoderare easy to forget for categorical features.- Calling
predicton a model trained on unscaled data and then scaling inputs at inference time silently degrades accuracy. train_test_split(X, y)requires both arrays;train_test_split(X)for unsupervised learning slips past static checkers.- Treating ML Multiple Regression as a black box without reading the docs — the API has subtle defaults that bite when you scale.
🚀 Performance & Best Practices
- Wrap preprocessing + estimator in a
Pipelinesofit/predictstay reproducible. - Use
joblib/picklefor serialising models, not the entire Python state. - For > 100 k rows, switch to
HistGradientBoostingClassifierorcuML. - When working with scikit-learn, prefer vectorised / batched operations over Python loops.
🧪 Try It Yourself
- Reproduce the snippet on a representative slice of your own data.
- Profile the snippet with
cProfileortimeitand find the single biggest improvement. - Generalise the snippet into a small, reusable function you can drop into future projects.