🎯 What You'll Learn
- Extend simple linear regression to multiple features with multivariate regression
- Identify when linear regression fails and a non-linear model is needed
- Transform features into polynomial features using
PolynomialFeaturesfrom sklearn - Understand interaction terms and when they add value
- Build a validation curve to select the right polynomial degree and preview the bias-variance tradeoff
pip install scikit-learn numpy matplotlib pandas. Key import: from sklearn.preprocessing import PolynomialFeatures. This lesson builds directly on the Pipeline and preprocessing concepts from Lesson 10.
1 Multivariate Linear Regression
In Lesson 11, we modeled one feature predicting one output: ŷ = β₀ + β₁x. Real datasets have many features. Multivariate linear regression extends the model to p features:
ŷ = β₀ + β₁x₁ + β₂x₂ + ... + βₚxₚ
In matrix form, where X has shape (n, p+1) with a bias column of ones prepended:
ŷ = X β where β = (XᵀX)⁻¹ Xᵀy (OLS solution)
The coefficient βⱼ is interpreted as: "a one-unit increase in xⱼ is associated with a βⱼ change in ŷ, holding all other features constant." This ceteris paribus interpretation is powerful for understanding feature effects in isolation.
import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.metrics import r2_score, mean_squared_error
np.random.seed(42)
n = 500
# Multivariate house price dataset
sqft = np.random.uniform(500, 5000, n)
bedrooms = np.random.randint(1, 7, n).astype(float)
bathrooms = np.random.uniform(1, 4, n)
age = np.random.uniform(0, 80, n)
# True relationship
price = (50000
+ 120 * sqft
+ 15000 * bedrooms
+ 25000 * bathrooms
- 500 * age
+ np.random.normal(0, 30000, n))
df = pd.DataFrame({
'sqft': sqft, 'bedrooms': bedrooms,
'bathrooms': bathrooms, 'age': age, '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)
pipeline = Pipeline([('scaler', StandardScaler()), ('model', LinearRegression())])
pipeline.fit(X_train, y_train)
print(f"Train R²: {pipeline.score(X_train, y_train):.4f}")
print(f"Test R²: {pipeline.score(X_test, y_test):.4f}")
# Interpret coefficients (unscale them)
scaler = pipeline.named_steps['scaler']
lr = pipeline.named_steps['model']
print("\nCoefficients (per unit of original feature):")
for name, coef, scale in zip(X.columns, lr.coef_, scaler.scale_):
coef_orig = coef / scale
print(f" {name:12s}: ${coef_orig:,.0f}")
# sqft: $120/sqft
# bedrooms: $15,000/bedroom
# bathrooms: $25,000/bathroom
# age: -$500/yearIf two features are highly correlated (e.g., "sqft" and "num_rooms" — larger houses have more rooms), it becomes difficult for the model to separate their individual effects. Coefficients become unstable — small changes to the data cause large changes to the coefficient values. Symptoms: absurdly large or opposite-sign coefficients. The simplest fix available to you now: remove one of the correlated features. Later lessons add two more options — Ridge regression (Lesson 21) tolerates correlated features gracefully, and PCA (Lesson 31) can merge them into one.
2 When Linearity Fails
A linear model assumes the relationship between x and y is a straight line. Many real-world relationships are not linear:
- Diminishing returns: a 500 sqft increase means more for a 1000 sqft house than for a 5000 sqft house.
- Saturation effects: adding the 10th bedroom increases price less than adding the 2nd bedroom.
- U-shaped relationships: product reviews have a U-shaped relationship with sales (very short and very long reviews both do poorly; medium-length reviews perform best).
- Threshold effects: salary jumps sharply at 5 years of experience (senior promotion).
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score
np.random.seed(42)
X = np.linspace(0, 6, 100)
y = 2 * np.sin(X) + np.random.normal(0, 0.3, 100) # non-linear: sine wave
X_2d = X.reshape(-1, 1)
lr = LinearRegression().fit(X_2d, y)
y_pred_linear = lr.predict(X_2d)
print(f"Linear regression R²: {r2_score(y, y_pred_linear):.4f}")
# R² ≈ 0.16 — very poor fit! The line can't capture the wave shape.
# Visualize the mismatch
fig, ax = plt.subplots(figsize=(8, 4))
ax.scatter(X, y, alpha=0.5, color='steelblue', label='Data')
ax.plot(X, y_pred_linear, 'r-', linewidth=2, label=f'Linear fit (R²={r2_score(y, y_pred_linear):.2f})')
ax.set_title('Linear Regression Fails on Non-Linear Data')
ax.legend()
plt.tight_layout()
plt.savefig('linearity_failure.png', dpi=100)When you plot residuals (actual − predicted) versus the fitted values and see a systematic pattern (a U-shape, wave, or funnel), that's a signal that the linearity assumption is violated and polynomial or other non-linear features are needed.
3 Polynomial Features with sklearn
PolynomialFeatures transforms a set of input features into all polynomial combinations up to a specified degree. For a single feature x with degree=2, it creates [1, x, x²]. For degree=3, it creates [1, x, x², x³].
from sklearn.preprocessing import PolynomialFeatures
import numpy as np
# Single feature
X_1d = np.array([[2], [3], [4]])
poly_d2 = PolynomialFeatures(degree=2, include_bias=True)
X_poly_d2 = poly_d2.fit_transform(X_1d)
print("Degree 2, single feature:")
print("Input: ", X_1d.ravel())
print("Output:", X_poly_d2)
# [[1. 2. 4.] → [1, x, x²]
# [1. 3. 9.]
# [1. 4. 16.]]
poly_d3 = PolynomialFeatures(degree=3, include_bias=False) # no bias column
X_poly_d3 = poly_d3.fit_transform(X_1d)
print("\nDegree 3, no bias:")
print(X_poly_d3)
# [[ 2. 4. 8.] → [x, x², x³]
# [ 3. 9. 27.]
# [ 4. 16. 64.]]
# Feature names
print("\nFeature names:", poly_d3.get_feature_names_out(['x']))Building a Polynomial Regression Pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import Pipeline
from sklearn.metrics import r2_score
import numpy as np
# Non-linear data: y = 0.5x² - 2x + 1 + noise
np.random.seed(42)
X_raw = np.sort(np.random.uniform(-3, 3, 80))
y_raw = 0.5 * X_raw**2 - 2 * X_raw + 1 + np.random.normal(0, 0.5, 80)
X_2d = X_raw.reshape(-1, 1)
# Fit with different degrees and compare
for degree in [1, 2, 3, 5]:
pipe = Pipeline([
('poly', PolynomialFeatures(degree=degree, include_bias=False)),
('scaler', StandardScaler()),
('model', LinearRegression()),
])
pipe.fit(X_2d, y_raw)
r2 = pipe.score(X_2d, y_raw)
print(f"Degree {degree}: R² = {r2:.4f}")
# Degree 1: R² = 0.6521 — underfitting
# Degree 2: R² = 0.9734 — good fit!
# Degree 3: R² = 0.9741 — marginal improvement
# Degree 5: R² = 0.9762 — very slight improvement (risk of overfitting on test)4 Underfitting vs Overfitting with Polynomial Degree
Degree is a hyperparameter controlling model complexity. Too low → underfitting. Too high → overfitting. The right degree captures the true pattern without memorizing noise.
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import Pipeline
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score
np.random.seed(0)
X_all = np.sort(np.random.uniform(-3, 3, 30))
y_all = 0.5 * X_all**2 - 2 * X_all + 1 + np.random.normal(0, 0.8, 30)
X_train, X_test, y_train, y_test = train_test_split(
X_all.reshape(-1, 1), y_all, test_size=0.33, random_state=42
)
X_plot = np.linspace(-3, 3, 300).reshape(-1, 1)
fig, axes = plt.subplots(1, 3, figsize=(14, 4), sharey=True)
degrees = [1, 2, 15]
titles = ['Degree 1 (Underfitting)', 'Degree 2 (Good Fit)', 'Degree 15 (Overfitting)']
for ax, degree, title in zip(axes, degrees, titles):
pipe = Pipeline([
('poly', PolynomialFeatures(degree=degree, include_bias=False)),
('scale', StandardScaler()),
('lr', LinearRegression()),
])
pipe.fit(X_train, y_train)
y_plot = pipe.predict(X_plot)
train_r2 = pipe.score(X_train, y_train)
test_r2 = pipe.score(X_test, y_test)
ax.scatter(X_train, y_train, color='steelblue', alpha=0.6, s=30, label='Train')
ax.scatter(X_test, y_test, color='orange', alpha=0.6, s=30, label='Test')
ax.plot(X_plot, y_plot, 'red', linewidth=2)
ax.set_title(f'{title}\nTrain R²={train_r2:.2f}, Test R²={test_r2:.2f}')
ax.set_ylim(-8, 15)
ax.legend(fontsize=8)
plt.tight_layout()
plt.savefig('poly_degrees.png', dpi=100)
print("Saved: poly_degrees.png")
# Degree 1: Train R²≈0.65, Test R²≈0.62 — both low → underfitting
# Degree 2: Train R²≈0.96, Test R²≈0.94 — both high → good fit
# Degree 15: Train R²≈0.99, Test R²≈-2.1 — gap huge → overfittingA degree-15 polynomial will look reasonable within the range of training data but can produce wildly wrong predictions just slightly outside that range. High-degree polynomial models have extreme extrapolation behavior. Always evaluate on a held-out test set, and be cautious about using high-degree polynomial models for predictions beyond the training data's range.
Try It: Drag the Degree Slider and Watch the Fit Break
Below is the same kind of noisy non-linear dataset as above (y = 0.5x² − 2x + 1 + noise, 24 points), fit live in your browser with an ordinary least-squares polynomial fit (a Vandermonde matrix solved via the normal equation — no sklearn needed). Drag the degree slider from 1 to 12:
- Degree 1–2: at degree 1 the line can't bend at all and misses the curvature (underfitting, high train MSE). Degree 2 matches the true quadratic shape almost exactly — the global minimum of training error you'd actually want.
- Degree 6–12: the curve starts threading through every individual point, looping and oscillating wildly between them — especially near the edges of the data. Train MSE keeps dropping toward zero even as the curve becomes visually absurd. That gap between "fits training points perfectly" and "represents the true relationship" is overfitting.
Degree 1 — a straight line cannot capture the curvature: high bias, high train MSE.
5 Feature Interaction Terms
When PolynomialFeatures is applied to multiple features, it also generates interaction terms — products of pairs of features. An interaction term x₁·x₂ captures the idea that the effect of x₁ on y depends on the value of x₂ (and vice versa).
For 2 features [x₁, x₂] with degree=2, PolynomialFeatures generates: [1, x₁, x₂, x₁², x₁x₂, x₂²] — 6 features total.
from sklearn.preprocessing import PolynomialFeatures
import numpy as np
# 2 features, degree=2
X_2feat = np.array([[2, 3], [4, 1], [1, 5]])
poly = PolynomialFeatures(degree=2, include_bias=True)
X_poly = poly.fit_transform(X_2feat)
print("Feature names:", poly.get_feature_names_out(['sqft', 'bedrooms']))
# ['1', 'sqft', 'bedrooms', 'sqft^2', 'sqft bedrooms', 'bedrooms^2']
print("\nTransformed matrix:")
print(X_poly)
# Interaction only (no squared terms)
poly_interact = PolynomialFeatures(degree=2, interaction_only=True, include_bias=False)
X_interact = poly_interact.fit_transform(X_2feat)
print("\nInteraction only features:", poly_interact.get_feature_names_out(['sqft', 'bedrooms']))
# ['sqft', 'bedrooms', 'sqft bedrooms']When Interactions Matter
A concrete example: in house pricing, the effect of additional square footage might be stronger when the house already has many bedrooms. The interaction term sqft × bedrooms captures this combined effect.
import numpy as np
import pandas as pd
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import Pipeline
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score
np.random.seed(42)
n = 400
sqft = np.random.uniform(500, 4000, n)
bedrooms = np.random.randint(1, 6, n).astype(float)
# True model HAS an interaction: price depends on sqft × bedrooms
price = (40000
+ 80 * sqft
+ 12000 * bedrooms
+ 0.05 * sqft * bedrooms # interaction term
+ np.random.normal(0, 25000, n))
X = np.column_stack([sqft, bedrooms])
y = price
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
# Without interaction terms
pipe_no_interact = Pipeline([
('scale', StandardScaler()),
('lr', LinearRegression()),
])
pipe_no_interact.fit(X_train, y_train)
r2_no = pipe_no_interact.score(X_test, y_test)
# With interaction terms (degree=2)
pipe_interact = Pipeline([
('poly', PolynomialFeatures(degree=2, include_bias=False)),
('scale', StandardScaler()),
('lr', LinearRegression()),
])
pipe_interact.fit(X_train, y_train)
r2_int = pipe_interact.score(X_test, y_test)
print(f"R² without interactions: {r2_no:.4f}")
print(f"R² with interactions: {r2_int:.4f}")
# R² without: 0.8914
# R² with: 0.9287 ← improvement from capturing the interaction!6 Interpreting Polynomial Coefficients
One cost of using polynomial features is reduced interpretability. In simple linear regression, β₁ = 120 means "each additional sqft adds $120." With polynomial features, the coefficient for sqft² cannot be interpreted in the same direct way — because the total effect of sqft now depends on both the linear and squared terms together.
import numpy as np
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import Pipeline
np.random.seed(42)
X = np.sort(np.random.uniform(0, 10, 100)).reshape(-1, 1)
y = 3 * X.ravel() - 0.2 * X.ravel()**2 + np.random.randn(100) * 2
pipe = Pipeline([
('poly', PolynomialFeatures(degree=2, include_bias=False)),
('scale', StandardScaler()),
('lr', LinearRegression()),
])
pipe.fit(X, y)
# Raw coefficients in standardized poly space — not directly interpretable
lr = pipe.named_steps['lr']
print("Intercept:", lr.intercept_)
print("Coefficients:", lr.coef_)
# Better approach: interpret by comparing predictions at specific x values
x_values = np.array([[1], [2], [5], [8]])
predictions = pipe.predict(x_values)
print("\nPredictions at specific x values:")
for x, p in zip(x_values.ravel(), predictions):
print(f" x={x}: ŷ={p:.2f}")
# For polynomial models, the marginal effect of x DEPENDS on the current value of x
# dy/dx = β₁ + 2β₂x (for degree-2 polynomial y = β₀ + β₁x + β₂x²)
# The effect of x is not constant — it changes as x changesDegree-2 polynomial models are usually still interpretable: you can say "price increases with sqft, but the rate of increase slows down at large sqft values (diminishing returns)." Degree-5+ models lose clear interpretability. If you need high accuracy and interpretability, consider SHAP values to explain individual predictions without relying on coefficients.
7 Choosing Polynomial Degree with Validation Curves
To find the best polynomial degree, we use a validation curve: train models with different degrees and compare training and validation performance. This directly visualises the bias-variance tradeoff (explored deeply in Lesson 14).
We need a fair way to score each degree without touching the test set. One option is to carve a chunk off the training data as a "validation set" — but then the score depends on which lucky (or unlucky) chunk you picked. K-fold cross-validation fixes this: split the training data into K equal folds (here K=5), train on 4 folds and score on the 5th, then rotate so every fold gets a turn as the scorer, and average the 5 scores. Every training sample gets used for both training and validation, and no single unlucky split can mislead you. In scikit-learn this whole procedure is one call: cross_val_score(model, X_train, y_train, cv=5).
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import Pipeline
from sklearn.model_selection import cross_val_score, train_test_split
from sklearn.metrics import r2_score
np.random.seed(42)
X_all = np.sort(np.random.uniform(-3, 3, 120))
y_all = (0.5 * X_all**2 - 2 * X_all + 1
+ np.random.normal(0, 0.8, 120))
X_2d = X_all.reshape(-1, 1)
X_train, X_test, y_train, y_test = train_test_split(X_2d, y_all, test_size=0.25, random_state=42)
degrees = range(1, 12)
train_r2s = []
val_r2s = []
for degree in degrees:
pipe = Pipeline([
('poly', PolynomialFeatures(degree=degree, include_bias=False)),
('scale', StandardScaler()),
('lr', LinearRegression()),
])
# 5-fold CV on training set for validation score
cv_scores = cross_val_score(pipe, X_train, y_train, cv=5, scoring='r2')
pipe.fit(X_train, y_train)
train_r2s.append(pipe.score(X_train, y_train))
val_r2s.append(cv_scores.mean())
best_degree = degrees[np.argmax(val_r2s)]
print(f"Best degree by CV validation R²: {best_degree}")
# Best degree: 2
# Print comparison table
print(f"\n{'Degree':>8} {'Train R²':>10} {'CV Val R²':>10}")
for d, tr, vr in zip(degrees, train_r2s, val_r2s):
marker = ' ← best' if d == best_degree else ''
print(f"{d:>8} {tr:>10.4f} {vr:>10.4f}{marker}")You might be tempted to try all degrees and pick the one with the best test R². But this turns the test set into a validation set, and the test performance becomes an optimistic estimate. Use cross-validation on the training set to pick the degree; then evaluate the chosen degree once on the held-out test set. The test set is used exactly once.
8 Practical Full Pipeline: Scale → Poly → Regress
The complete workflow for polynomial regression on real data, combining all concepts from this lesson:
import numpy as np
import pandas as pd
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import Pipeline
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.metrics import r2_score, mean_squared_error
# Build a multivariate dataset with non-linear relationships
np.random.seed(42)
n = 600
sqft = np.random.uniform(400, 5000, n)
lot = np.random.uniform(2000, 20000, n)
age = np.random.uniform(0, 100, n)
# Non-linear true model: sqft has diminishing returns (sqft²), age negatively affects price
price = (20000
+ 100 * sqft # linear sqft effect
- 0.01 * sqft**2 # diminishing returns at large sqft
+ 3 * lot
- 300 * age
+ np.random.normal(0, 40000, n))
X = np.column_stack([sqft, lot, age])
y = price
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
# Option 1: Fixed degree=2 pipeline
pipe_d2 = Pipeline([
('scaler', StandardScaler()),
('poly', PolynomialFeatures(degree=2, include_bias=False)),
('model', LinearRegression()),
])
pipe_d2.fit(X_train, y_train)
print(f"Degree-2 Test R²: {pipe_d2.score(X_test, y_test):.4f}")
# Option 2: choose the degree honestly with cross-validation
# (the Section 7 recipe, now on a realistic multivariate dataset)
best_degree, best_cv_r2 = None, -np.inf
for degree in [1, 2, 3]:
pipe = Pipeline([
('scaler', StandardScaler()),
('poly', PolynomialFeatures(degree=degree, include_bias=False)),
('model', LinearRegression()),
])
cv_r2 = cross_val_score(pipe, X_train, y_train, cv=5, scoring='r2').mean()
print(f"degree={degree}: CV R² = {cv_r2:.4f}")
if cv_r2 > best_cv_r2:
best_cv_r2, best_degree = cv_r2, degree
# Refit the winning degree on all training data; test set used exactly once
final_pipe = Pipeline([
('scaler', StandardScaler()),
('poly', PolynomialFeatures(degree=best_degree, include_bias=False)),
('model', LinearRegression()),
])
final_pipe.fit(X_train, y_train)
print(f"\nBest degree: {best_degree}")
print(f"Test R²: {final_pipe.score(X_test, y_test):.4f}")
# Looking ahead: Lesson 21 adds "regularization" to tame wild high-degree
# polynomials, and Lesson 27 introduces GridSearchCV — a tool that runs
# search loops like the one above for you, over any number of settings
# Note: number of features after PolynomialFeatures(degree=2) with 3 input features:
# 1 + 3 + 3 + 3 = (bias) + (linear) + (squared) + (interactions) = 10 features total
poly_temp = PolynomialFeatures(degree=2, include_bias=True)
print(f"\nOutput features for 3 inputs, degree=2: {poly_temp.fit_transform(X_train[:1]).shape[1]}")
print(f"Feature names: {poly_temp.get_feature_names_out(['sqft', 'lot', 'age'])}")Real-World Spotlight: House Price Prediction with Diminishing Returns
import numpy as np
import pandas as pd
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import Pipeline
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score, mean_absolute_error
np.random.seed(42)
n = 1000
sqft = np.random.uniform(500, 6000, n)
# True relationship: sqrt(sqft) captures diminishing returns more naturally
# but polynomial degree=2 is also a good approximation
price = (80000 + 200 * np.sqrt(sqft) * 30 + np.random.normal(0, 40000, n))
X = sqft.reshape(-1, 1)
y = price
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
# Compare: linear vs degree-2 polynomial
results = {}
for degree, name in [(1, 'Linear (degree 1)'), (2, 'Polynomial (degree 2)'), (3, 'Polynomial (degree 3)')]:
pipe = Pipeline([
('poly', PolynomialFeatures(degree=degree, include_bias=False)),
('scaler', StandardScaler()),
('lr', LinearRegression()),
])
pipe.fit(X_train, y_train)
y_pred = pipe.predict(X_test)
results[name] = {
'train_r2': pipe.score(X_train, y_train),
'test_r2': pipe.score(X_test, y_test),
'test_mae': mean_absolute_error(y_test, y_pred),
}
print(f"{'Model':<25} {'Train R²':>9} {'Test R²':>9} {'Test MAE':>12}")
print('-' * 60)
for name, metrics in results.items():
print(f"{name:<25} {metrics['train_r2']:>9.4f} {metrics['test_r2']:>9.4f} ${metrics['test_mae']:>10,.0f}")
# Linear (degree 1) 0.8931 0.8877 $51,234
# Polynomial (degree 2) 0.9412 0.9378 $38,201 ← significantly better
# Polynomial (degree 3) 0.9421 0.9370 $38,518 ← marginal + slight overfit riskThe degree-2 polynomial captures the diminishing returns relationship far better than the linear model, improving R² from 0.89 to 0.94 and MAE from $51K to $38K. The degree-3 model barely improves training R² and slightly degrades test R² — a warning sign of overfitting that will be explored in depth in Lesson 14.
Quick Check
✍️ Practice Exercises
- Load the Boston Housing or California Housing dataset. Plot each feature against the target variable. Identify at least 2 features that appear non-linear in their relationship with price. Confirm by comparing linear vs degree-2 polynomial R² for those features individually.
- Starting with 3 input features, compute how many features PolynomialFeatures generates for degree=2, degree=3, and degree=4 (using the formula C(p+d, d) where p=features, d=degree). Verify your answer with sklearn.
- Build pipelines for polynomial degrees 1, 2, and 3 (scaler → PolynomialFeatures → LinearRegression) and compare them with 5-fold
cross_val_scoreon the training set. Refit the winner on the full training set and report its final test R². - Create a dataset with a true interaction effect: price = 100*sqft + 50*bedrooms + 0.1*sqft*bedrooms + noise. Show that a model WITHOUT the interaction term has systematically wrong predictions for large houses with many bedrooms. Add PolynomialFeatures(degree=2) to capture the interaction and measure the improvement.
📚 Primary Source for This Lesson
sklearn User Guide: Polynomial Regression
The official documentation section on using PolynomialFeatures as a transformer with Pipeline. Also see Chapter 7 of An Introduction to Statistical Learning (free at statlearning.com) on "Moving Beyond Linearity" which covers polynomial regression, splines, and local regression in depth.