Linear Regression
Linear regression is the "Hello World" of ML — simple, interpretable, and surprisingly effective. It's the starting point for every regression task.
1. What You'll Learn
- The math behind linear regression: hypothesis function, loss function (MSE), gradient descent
- Scikit-learn implementation: LinearRegression vs SGDRegressor, regularization (Ridge/Lasso/ElasticNet)
- Multiple regression and feature interpretation: the economic meaning of coefficients
- Hypothesis testing and diagnostics: residual analysis, Q-Q plots, heteroscedasticity, multicollinearity (VIF)
- Bob's sales forecast: predicting next month's revenue using historical ad spend + traffic data
2. A Real-World E-Commerce Story
(1) The Pain Point: No Idea How to Split the Ad Budget
Bob has a monthly ad budget of 200 thousand USD to allocate across Google Ads, Facebook, and Email. His gut says "spend more, sell more," but last month an extra 50k USD on Google Ads only brought in 30k USD of additional revenue — diminishing returns, but he has no idea where the tipping point is.
(2) The Linear Regression Solution
Linear regression can quantify each channel's contribution to revenue and find the optimal budget allocation.
PYTHON
from sklearn.linear_model import LinearRegression
model = LinearRegression()
model.fit(X_train, y_train)
# Each coefficient = marginal contribution per 1k USD ad spend
for channel, coef in zip(channels, model.coef_):
print(f"{channel}: +{coef:.2f}k USD revenue per 1k USD spend")
(3) The Result: 35% ROI Improvement
Bob discovers that Google Ads has a coefficient of 0.8 (every 1k USD spent returns 0.8k USD), while Email has a coefficient of 2.5 (every 1k USD spent returns 2.5k USD). After reallocating the budget, overall ROI improves by 35%.
3. The Math Behind Linear Regression
(1) Hypothesis Function and Loss Function
Linear regression assumes the output is a linear combination of the inputs: $\hat{y} = w_1x_1 + w_2x_2 + ... + b$
The loss function uses MSE (Mean Squared Error): $L = \frac{1}{n}\sum_{i=1}^{n}(y_i - \hat{y}_i)^2$
▶ Example: Linear Regression from Scratch
PYTHON
import numpy as np
# Gradient descent for linear regression
def linear_regression_gd(X, y, lr=0.01, epochs=1000):
n_samples, n_features = X.shape
w = np.zeros(n_features)
b = 0.0
for epoch in range(epochs):
# Forward pass
y_pred = X @ w + b
# Compute gradients
dw = (2 / n_samples) * (X.T @ (y_pred - y))
db = (2 / n_samples) * np.sum(y_pred - y)
# Update parameters
w -= lr * dw
b -= lr * db
if epoch % 200 == 0:
loss = np.mean((y - y_pred) ** 2)
print(f"Epoch {epoch}: MSE = {loss:.4f}")
return w, b
# Test with simple data
rng = np.random.default_rng(42)
X = rng.uniform(0, 10, (100, 1))
y = 3 * X.squeeze() + 7 + rng.normal(0, 2, 100)
w, b = linear_regression_gd(X, y, lr=0.01, epochs=1000)
print(f"\nLearned: w={w[0]:.2f}, b={b:.2f}")
print(f"True: w=3.00, b=7.00")
Output:
TEXT
# Function defined successfully
(2) Gradient Descent Visualization
sequenceDiagram
participant Init as Parameter Init
participant Fwd as Forward Pass
participant Loss as Compute Loss
participant Grad as Compute Gradient
participant Update as Update Parameters
Init->>Fwd: w=0, b=0
loop Each Epoch
Fwd->>Loss: y_pred = Xw + b
Loss->>Grad: MSE = mean((y - y_pred)²)
Grad->>Update: dw, db = gradients
Update->>Fwd: w -= lr*dw, b -= lr*db
end
4. Scikit-learn Implementation
(1) LinearRegression vs SGDRegressor
▶ Example: Comparing Two Solvers
PYTHON
from sklearn.linear_model import LinearRegression, SGDRegressor
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
import numpy as np
rng = np.random.default_rng(42)
n = 500
X = rng.uniform(0, 100, (n, 3))
y = 50 + 0.8 * X[:, 0] + 1.2 * X[:, 1] - 0.5 * X[:, 2] + rng.normal(0, 5, n)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
# Method 1: Normal equation (closed-form, exact)
lr_normal = LinearRegression()
lr_normal.fit(X_train, y_train)
# Method 2: SGD (iterative, for large datasets)
lr_sgd = Pipeline([
("scaler", StandardScaler()),
("sgd", SGDRegressor(max_iter=1000, learning_rate="constant", eta0=0.01, random_state=42)),
])
lr_sgd.fit(X_train, y_train)
print(f"Normal Equation R²: {lr_normal.score(X_test, y_test):.4f}")
print(f"SGD R²: {lr_sgd.score(X_test, y_test):.4f}")
print(f"\nNormal coefficients: {lr_normal.coef_.round(2)}")
print(f"True coefficients: [0.80, 1.20, -0.50]")
Output:
TEXT
# Executed successfully
| Dimension | LinearRegression (Normal Equation) | SGDRegressor |
|---|---|---|
| Solver | Analytical solution $(X^TX)^{-1}X^Ty$ | Iterative gradient descent |
| Complexity | O(n³) in number of features | O(n) per step |
| Suitable data size | Small to medium (n < 100k) | Large (n > 100k) |
| Precision | Exact solution | Approximate solution |
| Requires scaling | No | Yes |
(2) Regularization: Ridge/Lasso/ElasticNet
▶ Example: Regularization Comparison
PYTHON
from sklearn.linear_model import Ridge, Lasso, ElasticNet
from sklearn.model_selection import cross_val_score
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
import numpy as np
rng = np.random.default_rng(42)
n, p = 100, 20 # p > relevant features (only 5 matter)
X = rng.standard_normal((n, p))
true_coefs = np.zeros(p)
true_coefs[:5] = [3, -2, 1.5, 0.8, -0.5]
y = X @ true_coefs + rng.normal(0, 1, n)
models = {
"LinearRegression": LinearRegression(),
"Ridge (alpha=1)": Ridge(alpha=1),
"Lasso (alpha=0.1)": Lasso(alpha=0.1),
"ElasticNet (alpha=0.1)": ElasticNet(alpha=0.1, l1_ratio=0.5),
}
for name, model in models.items():
pipe = Pipeline([("scaler", StandardScaler()), ("model", model)])
scores = cross_val_score(pipe, X, y, cv=5, scoring="r2")
pipe.fit(X, y)
nonzero = np.sum(np.abs(pipe.named_steps["model"].coef_) > 0.01)
print(f"{name:25s}: R²={scores.mean():.3f}, Non-zero coefs={nonzero}/{p}")
Output:
TEXT
# Executed successfully
| Regularization | Penalty Term | Effect | Use Case |
|---|---|---|---|
| Ridge (L2) | $\alpha \sum w_i^2$ | Shrinks coefficients toward zero | Many relevant features |
| Lasso (L1) | $\alpha \sum |w_i|$ | Sparsifies, drives some coefficients to zero | Feature selection |
| ElasticNet | L1 + L2 hybrid | Sparse + stable | Groups of correlated features |
5. Multiple Regression and Feature Interpretation
▶ Example: Bob's Ad ROI Analysis
PYTHON
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
import pandas as pd
import numpy as np
rng = np.random.default_rng(42)
n = 200
df = pd.DataFrame({
"google_ads_k": rng.uniform(20, 80, n),
"facebook_ads_k": rng.uniform(10, 50, n),
"email_marketing_k": rng.uniform(5, 30, n),
"traffic_k": rng.uniform(50, 300, n),
})
df["revenue_k"] = (
100
+ 0.8 * df["google_ads_k"]
+ 1.2 * df["facebook_ads_k"]
+ 2.5 * df["email_marketing_k"]
+ 0.15 * df["traffic_k"]
+ rng.normal(0, 10, n)
)
X = df.drop(columns=["revenue_k"])
y = df["revenue_k"]
# Train with standardized features for fair comparison
pipe = Pipeline([("scaler", StandardScaler()), ("model", LinearRegression())])
pipe.fit(X, y)
# Interpret standardized coefficients (importance ranking)
coefs = pipe.named_steps["model"].coef_
importance = pd.DataFrame({
"feature": X.columns,
"std_coef": coefs,
"abs_importance": np.abs(coefs),
}).sort_values("abs_importance", ascending=False)
print("Feature Importance (standardized coefficients):")
print(importance.to_string(index=False))
Output:
TEXT
Feature Importance (standardized coefficients):
Key takeaway: The magnitude of standardized coefficients directly reflects feature importance. The Email coefficient of 2.5 is far larger than Google's 0.8, meaning every 1k USD invested in Email returns 3x more than the same investment in Google Ads.
6. Regression Diagnostics
(1) Residual Analysis
▶ Example: 4-Panel Regression Diagnostics
PYTHON
import matplotlib.pyplot as plt
import numpy as np
from sklearn.linear_model import LinearRegression
from scipy import stats
rng = np.random.default_rng(42)
X = rng.uniform(0, 100, 200).reshape(-1, 1)
y = 50 + 0.8 * X.squeeze() + rng.normal(0, 5, 200)
model = LinearRegression().fit(X, y)
y_pred = model.predict(X)
residuals = y - y_pred
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
# (1) Residuals vs Fitted (check homoscedasticity)
axes[0, 0].scatter(y_pred, residuals, alpha=0.5, s=15)
axes[0, 0].axhline(0, color="red", linestyle="--")
axes[0, 0].set_xlabel("Fitted Values")
axes[0, 0].set_ylabel("Residuals")
axes[0, 0].set_title("Residuals vs Fitted")
# (2) Q-Q Plot (check normality)
stats.probplot(residuals, plot=axes[0, 1])
axes[0, 1].set_title("Q-Q Plot")
# (3) Scale-Location (check variance)
axes[1, 0].scatter(y_pred, np.sqrt(np.abs(residuals / np.std(residuals))), alpha=0.5, s=15)
axes[1, 0].set_xlabel("Fitted Values")
axes[1, 0].set_ylabel("Scale-Location")
axes[1, 0].set_title("Scale-Location")
# (4) Actual vs Predicted
axes[1, 1].scatter(y, y_pred, alpha=0.5, s=15)
axes[1, 1].plot([y.min(), y.max()], [y.min(), y.max()], "r--")
axes[1, 1].set_xlabel("Actual")
axes[1, 1].set_ylabel("Predicted")
axes[1, 1].set_title("Actual vs Predicted")
plt.tight_layout()
plt.savefig("regression_diagnostics.png", dpi=150)
Output:
TEXT
# Executed successfully
(2) Multicollinearity and VIF
▶ Example: VIF Detection
PYTHON
from statsmodels.stats.outliers_influence import variance_inflation_factor
import pandas as pd
import numpy as np
# VIF > 10 indicates severe multicollinearity
rng = np.random.default_rng(42)
n = 100
df = pd.DataFrame({
"ad_spend": rng.uniform(10, 100, n),
"traffic": rng.normal(0, 1, n) * 50 + 500,
"clicks": rng.normal(0, 1, n) * 100 + 1000, # Highly correlated with ad_spend
"revenue": rng.uniform(50, 500, n),
})
features = ["ad_spend", "traffic", "clicks"]
for i, col in enumerate(features):
vif = variance_inflation_factor(df[features].values, i)
print(f"{col:12s}: VIF = {vif:.2f} {'⚠️ HIGH' if vif > 10 else '✅ OK'}")
Output:
TEXT
# Executed successfully
| Diagnostic | Check Method | Normal Criteria | Fix for Violations |
|---|---|---|---|
| Linearity | Residuals vs Fitted plot | Random scatter | Add polynomial terms |
| Normality | Q-Q plot | Approximately straight line | Log transformation |
| Homoscedasticity | Scale-Location plot | Horizontal band | Weighted regression |
| Multicollinearity | VIF | VIF < 10 | Remove or merge features |
❓ FAQ
Q Does LinearRegression require feature scaling?
A It can train without scaling, but coefficients won't be comparable (different units). After standardization, coefficient magnitudes directly reflect feature importance. Ridge/Lasso/SGD require scaling.
Q Should I choose Ridge or Lasso?
A Many features that may all contribute → Ridge. Need automatic feature selection → Lasso. Groups of highly correlated features → ElasticNet. In practice, try all three and use cross-validation to pick the best.
Q How do I choose the alpha parameter?
A Use GridSearchCV or RidgeCV/LassoCV for automatic selection. Larger alpha means stronger regularization, pushing coefficients closer to zero. Typically search on a log scale: [0.001, 0.01, 0.1, 1, 10, 100].
Q How do I interpret a negative regression coefficient?
A A negative coefficient means the target decreases as that feature increases. For example, if "return rate" has a coefficient of -5, it means revenue drops by 5k USD for every 1% increase in return rate.
Q What if VIF > 10?
A Remove the feature with the highest VIF and recalculate. Alternatively, use PCA to eliminate collinearity. Or use Ridge regression (L2 regularization is naturally resistant to multicollinearity).
Q What if residuals aren't normally distributed?
A Try a log transformation on y (log(y)). If residuals show a funnel shape (heteroscedasticity), use weighted least squares or apply a Box-Cox transformation to y.
📖 Summary
- Linear regression assumes y is a linear combination of x, solved via MSE loss + gradient descent or the normal equation
- sklearn's unified API: LinearRegression (exact solution) vs SGDRegressor (large-scale data)
- The regularization trio: Ridge (L2 shrinks coefficients) / Lasso (L1 sparsifies) / ElasticNet (L1+L2 hybrid)
- Standardized coefficient magnitudes reflect feature importance and directly guide business decisions
- The four diagnostic checks: residual plot (linearity) / Q-Q plot (normality) / Scale-Location (variance) / VIF (collinearity)
- Interpreting coefficients in business terms is linear regression's greatest strength — explainability beats black-box models
📝 Exercises
- Basic (Difficulty ⭐): Use sklearn's LinearRegression to predict the California Housing dataset. Output R² and all feature coefficients. Hint:
fetch_california_housing()+model.coef_. - Intermediate (Difficulty ⭐⭐): Compare Ridge (alpha=0.1/1/10/100) on the same data using cross-validated R² and coefficient magnitudes. Observe how regularization strength affects coefficients. Hint: Use Pipeline(StandardScaler+Ridge) + GridSearchCV.
- Challenge (Difficulty ⭐⭐⭐): Implement a complete regression diagnostic pipeline — after training, plot the 4-panel diagnostic charts, compute VIF, determine whether linear regression assumptions hold, and propose improvements if they don't. Hint: Refer to the diagnostic code in Section 6.
← Previous: Comprehensive Exercise — Beginner Project | Next: Logistic Regression →
