Regression Models

Regression is one of the fundamental tasks in supervised machine learning, focused on predicting continuous numerical values rather than discrete categories. Think of it this way: while classification is like sorting mail into "spam" or "not spam" boxes, regression is like guessing the exact price tag on an item. Classification answers "what type?" questions, but regression answers "how much?" or "how many?" questions.

Here's an everyday example: when you estimate how long your commute will take based on the time of day and weather, you're doing regression in your head! Whether you're predicting house prices, forecasting tomorrow's temperature, estimating how much a customer will spend, or projecting next quarter's sales figures, regression algorithms are your go-to tools. In this comprehensive module, you'll learn the mathematical foundations behind regression (don't worry — we'll break it down step by step), implement various algorithms both from scratch and with scikit-learn, and understand how to choose and tune the right model for your specific problem.

Regression vs Classification

The key difference between regression and classification lies in the type of answer you're looking for. Classification asks "Which category does this belong to?" and outputs discrete labels like spam or not spam, cat or dog, positive or negative sentiment — essentially sorting objects into labeled bins. Regression, on the other hand, asks "How much?" or "What's the exact value?" and outputs a continuous number on a scale, like $350,000 for a house price or 72°F for tomorrow's temperature.

Here's a practical example that illustrates the distinction clearly. If you're a bank, asking "Will this customer pay back their loan?" is a classification problem with a Yes/No answer. But asking "How much interest should we charge this customer?" is a regression problem that outputs a percentage like 4.5% or 7.2%. This distinction matters because it affects the algorithms we use, how we measure success, and how we interpret results. Learning to recognize whether your problem is regression or classification is one of the most important skills in machine learning.

Real-World Applications

Regression algorithms power countless applications you interact with every day, often without realizing it. When you check Zillow for house prices, regression models estimate values based on square footage, location, bedrooms, and hundreds of other features. When Spotify predicts how long you'll listen, regression algorithms analyze your history to forecast streaming minutes. Amazon's delivery date estimates use regression to predict shipping time based on distance, warehouse inventory, and carrier performance. Even your weather app's temperature forecasts rely on regression models processing atmospheric data.

Beyond consumer applications, regression drives critical decisions across every industry. Healthcare uses it to predict patient recovery times and optimize drug dosages. Finance relies on regression for loan amount calculations and credit risk scoring. Insurance companies set premium prices using regression models, while manufacturers forecast production yields. Understanding regression opens doors to solving problems in virtually every domain where numerical predictions are needed.

Types of Regression

Just like choosing different tools from a toolbox depending on the job, there are several types of regression algorithms designed for different situations. Simple linear regression uses a single feature to predict with a straight line, perfect for understanding the basics. Multiple linear regression extends this to handle many features simultaneously, which is what most real-world problems require. Polynomial regression fits curved lines instead of straight ones, capturing non-linear relationships like diminishing returns or U-shaped patterns.

When models become too complex and start memorizing noise, regularization techniques come to the rescue. Ridge, Lasso, and Elastic Net add guardrails that prevent overfitting by penalizing large coefficients. For truly complex patterns that don't follow any mathematical formula, tree-based methods like Random Forest and Gradient Boosting offer powerful alternatives, though they sacrifice some interpretability. The key principle is to start simple with linear regression and only add complexity when you have clear evidence that simpler models aren't capturing the underlying patterns.

# Overview of regression types in scikit-learn
from sklearn.linear_model import LinearRegression, Ridge, Lasso, ElasticNet
from sklearn.preprocessing import PolynomialFeatures
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn.tree import DecisionTreeRegressor

# Simple/Multiple Linear Regression
linear_model = LinearRegression()

# Regularized regression (prevent overfitting)
ridge_model = Ridge(alpha=1.0)           # L2 regularization
lasso_model = Lasso(alpha=1.0)           # L1 regularization
elastic_model = ElasticNet(alpha=1.0)    # L1 + L2 combined

# For non-linear relationships
poly = PolynomialFeatures(degree=2)      # Add polynomial features
tree_model = DecisionTreeRegressor()     # Non-parametric
forest_model = RandomForestRegressor()   # Ensemble of trees
gbm_model = GradientBoostingRegressor()  # Boosted ensemble

print("Regression algorithms loaded successfully!")

Linear regression is the foundation of regression analysis and often the first algorithm every data scientist learns. Think of it as drawing the best-fit line through a scatter plot with a ruler, finding the single straight line that gets as close as possible to all data points simultaneously. Despite its simplicity, linear regression remains one of the most widely used techniques in production systems because it offers three crucial advantages: interpretability that lets you explain exactly what the model learned, computational speed that enables training in milliseconds even on large datasets, and surprising effectiveness on many real-world problems where linear relationships dominate.

The core mathematical idea is elegantly simple. You want to find the straight line that minimizes the total distance between your predictions and the actual values. This optimization problem has a clean, closed-form solution that computers can calculate directly without any iterative guessing. Understanding linear regression deeply will build your intuition for all the more complex algorithms you'll encounter later, as most advanced techniques either extend or address the limitations of this fundamental approach.

Simple Linear Regression

Simple linear regression models the relationship between a single input feature and the target variable using the familiar equation y = mx + b from algebra class. The slope m represents how much the output changes when the input increases by one unit, while the intercept b represents the baseline output when the input equals zero. For example, if predicting salary based on years of experience, a slope of 5000 means each additional year adds $5,000 to the predicted salary, and an intercept of 40000 means the starting salary with zero experience is $40,000.

The algorithm finds optimal values for these parameters by minimizing Mean Squared Error, which measures how far off each prediction is, squares those errors to make them all positive, and finds the line that makes this total as small as possible. The mathematical solution, known as the Normal Equation, provides exact values without any iterative guessing, making linear regression computationally efficient for datasets of any reasonable size.

# Simple Linear Regression from scratch
import numpy as np
import matplotlib.pyplot as plt

# Generate sample data: House size vs Price
np.random.seed(42)
X = np.random.rand(100) * 2000 + 500  # Size: 500-2500 sq ft
y = 50000 + 150 * X + np.random.randn(100) * 20000  # Price with noise

# Calculate coefficients using Normal Equation
X_mean, y_mean = np.mean(X), np.mean(y)
numerator = np.sum((X - X_mean) * (y - y_mean))
denominator = np.sum((X - X_mean) ** 2)

slope = numerator / denominator  # β₁
intercept = y_mean - slope * X_mean  # β₀

print(f"Equation: Price = {intercept:.2f} + {slope:.2f} × Size")
print(f"Interpretation: Each sq ft adds ${slope:.2f} to price")

Multiple Linear Regression

Real-world predictions rarely depend on a single factor. House prices depend on size, bedrooms, bathrooms, location, age, and dozens of other features simultaneously. Multiple linear regression extends the simple model to handle any number of input features, with each feature receiving its own coefficient that represents its independent contribution to the prediction. For example, a house price model might assign $100 per square foot, $15,000 per additional bedroom, and negative $500 for each year of age.

The key interpretive phrase is "holding all else constant." When the model assigns a $15,000 coefficient to bedrooms, it means that if two houses are identical in every other way but one has an extra bedroom, that house will be predicted to cost $15,000 more. This ability to isolate the effect of individual features while controlling for others makes multiple regression incredibly powerful for understanding what truly drives your predictions and for making actionable business recommendations.

# Multiple Linear Regression with scikit-learn
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.datasets import make_regression
import pandas as pd

# Create dataset with multiple features
X, y = make_regression(n_samples=200, n_features=4, noise=10, random_state=42)
feature_names = ['Size', 'Bedrooms', 'Age', 'Location_Score']
df = pd.DataFrame(X, columns=feature_names)
df['Price'] = y

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

# Train model
model = LinearRegression()
model.fit(X_train, y_train)

# View coefficients
print("Coefficients (feature importance):")
for name, coef in zip(feature_names, model.coef_):
    print(f"  {name}: {coef:.4f}")
print(f"\nIntercept: {model.intercept_:.4f}")
print(f"R² Score: {model.score(X_test, y_test):.4f}")

Assumptions of Linear Regression

Linear regression makes several key assumptions that should be validated for reliable results. The linearity assumption requires that the relationship between features and target actually follows a straight line pattern. Independence assumes each data point is unrelated to others, which is violated when dealing with time series data where today's value depends on yesterday's. Homoscedasticity means your prediction errors should be roughly the same size across all predictions, not larger for expensive houses and smaller for cheap ones.

The normality assumption suggests that residuals should follow a bell curve distribution, though mild violations are generally acceptable. Finally, multicollinearity should be avoided, meaning features shouldn't be highly correlated with each other like having both height in inches and height in centimeters. While linear regression is fairly robust to minor assumption violations, severe violations can make your coefficients untrustworthy and indicate you might need a different modeling approach.

Feature Scaling

When comparing features on vastly different scales, like age ranging from 0-100 versus income ranging from 0-1,000,000, the raw coefficients become impossible to compare meaningfully. Feature scaling solves this problem by transforming all features to a common scale. Standardization, the most common approach, transforms each feature to have a mean of 0 and standard deviation of 1, placing all values roughly in the -3 to +3 range.

After scaling, you can compare coefficients directly to determine which features matter most for your predictions. This becomes essential when using regularization techniques like Ridge and Lasso, which penalize large coefficients. Without scaling, regularization would unfairly penalize features that simply happen to have larger numerical values. Scikit-learn's StandardScaler handles this transformation easily, and wrapping it in a Pipeline ensures you never accidentally leak test data information into your preprocessing, a common mistake that leads to overly optimistic performance estimates.

# Feature scaling for fair coefficient comparison
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline

# Create a pipeline with scaling
pipeline = Pipeline([
    ('scaler', StandardScaler()),
    ('regressor', LinearRegression())
])

# Fit the pipeline
pipeline.fit(X_train, y_train)

# Get scaled coefficients (now comparable)
scaled_coefs = pipeline.named_steps['regressor'].coef_
print("Scaled Coefficients (comparable importance):")
for name, coef in zip(feature_names, scaled_coefs):
    print(f"  {name}: {coef:.4f}")

# Predictions work the same way
y_pred = pipeline.predict(X_test)
print(f"\nTest R² Score: {pipeline.score(X_test, y_test):.4f}")

Not all relationships in the real world follow straight lines. Consider how happiness increases rapidly when income grows from poverty to middle class, but additional millions barely register on the happiness scale — that's the curve of diminishing returns. Or consider plant growth versus sunlight, where too little sun kills the plant, the optimal amount helps it thrive, and too much sun causes burning, creating a U-shaped relationship. Fuel efficiency versus car speed follows a similar pattern, peaking around 55 mph then declining at higher speeds. When your data shows curves, bends, or U-shapes instead of straight lines, polynomial regression provides the solution.

The elegant trick behind polynomial regression is creating new features by raising existing features to various powers. By including x² and x³ alongside the original x, you enable linear regression to fit curved lines to your data. The result captures non-linear patterns that a simple straight line would completely miss, opening up a much wider range of real-world problems that can be modeled effectively.

From Linear to Polynomial

Polynomial regression is actually linear regression in disguise, with the transformation happening in the features rather than the algorithm. Starting with an original feature x, you create additional features by squaring it to get x², cubing it to get x³, and so on. These derived features join the original in a multiple linear regression, producing an equation like y = β₀ + β₁x + β₂x² that traces a parabola instead of a straight line.

Although the equation contains squared and cubed terms, it remains linear in the coefficients — we're still just multiplying and adding, which means all standard linear regression tools work perfectly. Higher polynomial degrees produce more flexible curves that can capture increasingly complex patterns, but this flexibility comes with serious risk. Too high a degree causes the model to overfit, memorizing random noise in the training data rather than learning the true underlying pattern. Finding the right balance between flexibility and generalization is the central challenge of polynomial regression.

# Polynomial Regression with scikit-learn
import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import Pipeline

# Generate non-linear data
np.random.seed(42)
X = np.linspace(0, 10, 100).reshape(-1, 1)
y = 2 + 3*X.flatten() - 0.5*X.flatten()**2 + np.random.randn(100)*3

# Compare different polynomial degrees
degrees = [1, 2, 3, 10]
plt.figure(figsize=(12, 3))

for i, degree in enumerate(degrees):
    plt.subplot(1, 4, i+1)
    
    # Create pipeline with polynomial features
    model = Pipeline([
        ('poly', PolynomialFeatures(degree=degree)),
        ('linear', LinearRegression())
    ])
    
    model.fit(X, y)
    y_pred = model.predict(X)
    
    plt.scatter(X, y, alpha=0.5, s=10)
    plt.plot(X, y_pred, 'r-', linewidth=2)
    plt.title(f'Degree {degree}')
    
plt.tight_layout()
plt.show()

Choosing the Right Degree

Selecting the optimal polynomial degree requires finding the balance between model simplicity and complexity. A degree that's too low results in underfitting, where the model is too simple to capture the true pattern in your data, like trying to describe a roller coaster's path with a straight line. A degree that's too high causes overfitting, where the model memorizes every tiny bump including random noise, fitting training data perfectly but failing miserably on new data.

The solution is cross-validation, which tests each candidate degree on data the model hasn't seen during training. You systematically try degrees 1, 2, 3, and beyond, measuring performance on held-out validation folds each time. The degree producing the lowest validation error, not training error, wins. Training error is deceptive because even a terrible overfitting model can memorize training data perfectly. Only validation performance reveals which degree will generalize well to new predictions.

# Finding optimal degree with cross-validation
from sklearn.model_selection import cross_val_score
import numpy as np

# Test degrees 1 through 10
degrees = range(1, 11)
cv_scores = []

for degree in degrees:
    model = Pipeline([
        ('poly', PolynomialFeatures(degree=degree)),
        ('linear', LinearRegression())
    ])
    
    # 5-fold cross-validation, using negative MSE (sklearn convention)
    scores = cross_val_score(model, X, y, cv=5, scoring='neg_mean_squared_error')
    cv_scores.append(-scores.mean())  # Convert to positive MSE

# Find best degree
best_degree = degrees[np.argmin(cv_scores)]
print(f"Best polynomial degree: {best_degree}")
print(f"CV MSE by degree:")
for d, score in zip(degrees, cv_scores):
    marker = " <-- Best" if d == best_degree else ""
    print(f"  Degree {d}: {score:.4f}{marker}")

Polynomial with Multiple Features

When you have multiple input features, polynomial transformation becomes significantly more complex because it creates not just powers of each feature but also interaction terms that multiply features together. For two features x₁ and x₂ with degree 2, PolynomialFeatures generates six terms: the constant 1, both original features x₁ and x₂, both squared terms x₁² and x₂², plus the interaction term x₁×x₂. These interaction terms capture situations where the effect of one feature depends on another's value, like how sunscreen effectiveness depends on both SPF level and sun exposure duration.

The number of generated features grows explosively with more input features and higher degrees. With 10 original features and degree 3, you end up with 286 polynomial features. With 20 features and degree 5, you would have over 53,000 features. This explosion makes regularization techniques like Ridge and Lasso absolutely essential for polynomial regression with multiple features, preventing the model from overfitting to the massive feature space.

# Polynomial features with multiple inputs
from sklearn.preprocessing import PolynomialFeatures
import pandas as pd

# Sample data with 2 features
X_sample = np.array([[2, 3], [4, 5]])
feature_names = ['x1', 'x2']

# Create polynomial features (degree 2)
poly = PolynomialFeatures(degree=2, include_bias=True)
X_poly = poly.fit_transform(X_sample)

# Show what features are created
poly_names = poly.get_feature_names_out(feature_names)
print("Original features:", feature_names)
print("Polynomial features:", list(poly_names))
print(f"\nFeature count: {len(feature_names)} -> {len(poly_names)}")
print("\nTransformed data:")
print(pd.DataFrame(X_poly, columns=poly_names))

When to Use Polynomial Regression

Use polynomial regression when you observe curved patterns in scatter plots, when domain knowledge suggests non-linear relationships (like diminishing returns), or when linear regression residuals show systematic patterns. However, consider alternatives: splines offer more flexibility with less overfitting risk, decision trees naturally capture non-linearity, and neural networks can learn complex patterns automatically. Polynomial regression shines when interpretability matters and the degree of non-linearity is moderate.

When models become too complex, they start memorizing every detail in the training data, including random noise that doesn't represent the true underlying pattern. This overfitting causes excellent performance on training data but miserable failure on new predictions. Regularization addresses this fundamental problem by adding a penalty term that discourages large coefficient values, effectively constraining the model's complexity and forcing it to focus on the most important patterns rather than every minor fluctuation.

The penalty works by making large coefficients expensive in terms of the optimization objective. Without regularization, a model might assign enormous weight to a single feature, claiming it's extraordinarily important. Regularization pushes back, requiring the model to justify large coefficients by demonstrably improving predictions. The three main techniques are Ridge (L2), which shrinks all coefficients toward zero while keeping all features active; Lasso (L1), which can shrink coefficients all the way to exactly zero, effectively removing unimportant features; and Elastic Net, which combines both approaches for situations where you're uncertain which strategy is best.

Ridge Regression (L2 Regularization)

Ridge regression adds a penalty proportional to the sum of squared coefficients, effectively shrinking all coefficients toward zero without ever reaching exactly zero. This approach works well when you have many features that each contribute at least something to the prediction and you don't want to eliminate any entirely. The alpha hyperparameter controls the regularization strength, with higher values producing smaller coefficients and simpler models, while alpha of zero reverts to standard linear regression.

Ridge particularly excels at handling multicollinearity, the problematic situation where features are highly correlated with each other. Standard linear regression might arbitrarily assign all the weight to one correlated feature while ignoring others, but Ridge distributes weight among correlated features more fairly. This produces more stable and reliable coefficient estimates that don't swing wildly when the training data changes slightly.

# Ridge Regression
from sklearn.linear_model import Ridge, RidgeCV
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.model_selection import train_test_split
import numpy as np

# Generate data with many features (some may be irrelevant)
np.random.seed(42)
X = np.random.randn(100, 20)  # 20 features
true_coef = np.array([1, 2, 0, 0, 0.5, 0, 0, -1, 0, 0.3] + [0]*10)
y = X @ true_coef + np.random.randn(100) * 0.5

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

# Ridge with different alpha values
alphas = [0.01, 0.1, 1, 10, 100]
print("Ridge Regression with different alpha values:")
for alpha in alphas:
    model = Pipeline([
        ('scaler', StandardScaler()),
        ('ridge', Ridge(alpha=alpha))
    ])
    model.fit(X_train, y_train)
    score = model.score(X_test, y_test)
    n_small_coef = np.sum(np.abs(model.named_steps['ridge'].coef_) < 0.1)
    print(f"  α={alpha:5}: R²={score:.4f}, Near-zero coefficients: {n_small_coef}")

Lasso Regression (L1 Regularization)

Lasso regression adds a penalty proportional to the sum of absolute coefficient values, creating a fundamentally different behavior than Ridge. Because of the mathematical properties of the absolute value function, Lasso can shrink coefficients all the way to exactly zero, effectively removing features from the model entirely. This automatic feature selection makes Lasso invaluable when you suspect many features are irrelevant noise and want the algorithm to identify which ones truly matter.

The mechanism works because the absolute value penalty creates corners in the optimization landscape where coefficients naturally settle at exactly zero. When you have 100 features but only 10 actually contribute to predictions, Lasso will identify those important features, assign them meaningful coefficients, and set the rest precisely to zero. This produces simpler, more interpretable models that focus only on the features that genuinely predict the outcome.

# Lasso Regression - Feature Selection
from sklearn.linear_model import Lasso, LassoCV

# Use same data as Ridge example
model = Pipeline([
    ('scaler', StandardScaler()),
    ('lasso', Lasso(alpha=0.1))
])
model.fit(X_train, y_train)

# Check which features were selected
coef = model.named_steps['lasso'].coef_
selected = np.where(coef != 0)[0]
eliminated = np.where(coef == 0)[0]

print(f"Lasso selected {len(selected)} features out of 20:")
print(f"  Selected features: {selected}")
print(f"  Eliminated features: {eliminated}")
print(f"\nNon-zero coefficients:")
for i in selected:
    print(f"  Feature {i}: {coef[i]:.4f}")
print(f"\nTest R²: {model.score(X_test, y_test):.4f}")

Elastic Net (L1 + L2 Combined)

Elastic Net combines both L1 and L2 penalties, blending the feature selection capability of Lasso with Ridge's ability to handle correlated features gracefully. The l1_ratio parameter controls the mix between the two approaches, with a value of 0 producing pure Ridge behavior, a value of 1 producing pure Lasso, and intermediate values creating a hybrid. This flexibility makes Elastic Net particularly valuable when you're uncertain which regularization approach suits your problem best.

Elastic Net excels when your dataset contains groups of correlated features, such as multiple measurements of the same underlying phenomenon. Pure Lasso tends to arbitrarily select just one feature from each correlated group while ignoring the others, but Elastic Net handles these groups more intelligently by either including or excluding entire groups together. Scikit-learn's ElasticNetCV class can automatically tune both the regularization strength alpha and the l1_ratio using cross-validation, removing the guesswork from hyperparameter selection.

# Elastic Net - Combining L1 and L2
from sklearn.linear_model import ElasticNet, ElasticNetCV

# Try different L1 ratios
l1_ratios = [0.1, 0.5, 0.7, 0.9, 0.95]
print("Elastic Net with different L1 ratios:")
print("(0 = Ridge, 1 = Lasso)")

for l1_ratio in l1_ratios:
    model = Pipeline([
        ('scaler', StandardScaler()),
        ('elastic', ElasticNet(alpha=0.1, l1_ratio=l1_ratio))
    ])
    model.fit(X_train, y_train)
    n_zero = np.sum(model.named_steps['elastic'].coef_ == 0)
    score = model.score(X_test, y_test)
    print(f"  l1_ratio={l1_ratio}: R²={score:.4f}, Zero coefs: {n_zero}")

# Use cross-validation to find best parameters
elastic_cv = ElasticNetCV(l1_ratio=[0.1, 0.5, 0.7, 0.9, 0.95], cv=5)
elastic_cv.fit(X_train, y_train)
print(f"\nBest l1_ratio: {elastic_cv.l1_ratio_}")
print(f"Best alpha: {elastic_cv.alpha_:.4f}")

Comparing Regularization Methods

Choosing the right regularization method depends on your specific problem. Use Ridge when you believe all features are relevant and want to reduce their impact without eliminating any. Use Lasso when you need feature selection and expect many features to be irrelevant. Use Elastic Net when you have correlated features and need both selection and stability. In practice, cross-validation is the best way to determine which method and hyperparameter values work best for your data.

Building a regression model is only half the battle — now you need to grade it! Think of evaluation metrics like a report card for your model. But here's the tricky part: unlike classification where "82% accuracy" is pretty clear, regression has multiple grades that each measure something different.

Regression evaluation relies on four fundamental metrics, each revealing different aspects of model performance. Mean Squared Error (MSE) penalizes large mistakes heavily, making it sensitive to outliers since one huge error hurts more than many small ones. Its companion metric, Root Mean Squared Error (RMSE), provides the same information but in human-readable units like dollars or degrees. Mean Absolute Error (MAE) represents the average size of mistakes while treating all errors equally, making it more robust to outliers. Finally, R² (R-squared) answers the question of how much pattern the model captured, ranging from 1.0 for perfect predictions to 0 for a model no better than guessing the mean.

Why use multiple metrics? Because each tells a different story about your model's behavior. An RMSE of $25,000 sounds concerning for a $100,000 house but excellent for a $5,000,000 mansion. R² provides a scale-independent percentage of explained variance, allowing you to compare models across vastly different prediction domains. Understanding when to prioritize each metric is essential for proper model evaluation.

Mean Squared Error (MSE) & RMSE

Mean Squared Error is the most popular regression metric, and understanding why it squares errors is key. Squaring serves two purposes: it makes negative errors positive (otherwise +10 and -10 would cancel to zero!) and acts like a megaphone for big mistakes, where an error of 10 becomes 100 but an error of 100 becomes 10,000. This property makes MSE particularly valuable when large mistakes are catastrophic, such as in medical diagnosis or bridge engineering applications. However, this sensitivity cuts both ways — a single outlier can dominate your entire MSE score and make it misleading.

RMSE (Root Mean Squared Error) is simply the square root of MSE, which converts the metric back to the original units. If you're predicting house prices in dollars, RMSE is also in dollars, making interpretation intuitive. An RMSE of $25,000 means "your predictions are typically off by about $25K" — much easier to understand than an MSE of 625,000,000. Use RMSE as your primary metric when large errors are disproportionately costly to your application.

# Calculating regression metrics
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
import numpy as np

# Example predictions
y_true = np.array([100, 150, 200, 250, 300])
y_pred = np.array([110, 145, 195, 260, 280])

# Calculate metrics
mse = mean_squared_error(y_true, y_pred)
rmse = np.sqrt(mse)
mae = mean_absolute_error(y_true, y_pred)
r2 = r2_score(y_true, y_pred)

print("Regression Metrics:")
print(f"  MSE:  {mse:.2f}")
print(f"  RMSE: {rmse:.2f}")
print(f"  MAE:  {mae:.2f}")
print(f"  R²:   {r2:.4f}")

# Detailed error analysis
errors = y_true - y_pred
print(f"\nError Analysis:")
print(f"  Errors: {errors}")
print(f"  Mean Error: {np.mean(errors):.2f} (bias)")
print(f"  Max Abs Error: {np.max(np.abs(errors)):.2f}")

Mean Absolute Error (MAE)

MAE is the "what you see is what you get" metric, simply measuring how far off your predictions are on average. Unlike MSE, MAE doesn't square the errors — an error of $100 counts as exactly $100, not $10,000. This property makes MAE far more robust to outliers that would otherwise distort your evaluation. The interpretability is also superior: an MAE of $15,000 directly means "on average, I'm off by $15K" with no mental math required.

Choosing between MAE and RMSE depends on your error philosophy. If a $100 error is simply 10x worse than a $10 error, use MAE since it treats all errors proportionally. However, if a $100 error is 100x worse than a $10 error because big mistakes are catastrophic in your domain, RMSE is the appropriate choice. MAE also shines when your data contains outliers that would unfairly inflate MSE, giving you a more representative picture of typical model performance.

# Comparing MSE vs MAE sensitivity to outliers
import numpy as np
from sklearn.metrics import mean_squared_error, mean_absolute_error

# Normal predictions
y_true_normal = np.array([100, 105, 110, 115, 120])
y_pred_normal = np.array([102, 103, 112, 118, 119])

# Same but with one outlier prediction
y_true_outlier = np.array([100, 105, 110, 115, 120])
y_pred_outlier = np.array([102, 103, 112, 118, 170])  # 170 is way off!

print("Without outlier:")
print(f"  MSE:  {mean_squared_error(y_true_normal, y_pred_normal):.2f}")
print(f"  MAE:  {mean_absolute_error(y_true_normal, y_pred_normal):.2f}")

print("\nWith outlier (170 vs 120):")
print(f"  MSE:  {mean_squared_error(y_true_outlier, y_pred_outlier):.2f}")
print(f"  MAE:  {mean_absolute_error(y_true_outlier, y_pred_outlier):.2f}")

print("\n→ MSE increased 100x, MAE only 10x - MSE is more outlier-sensitive!")

R² Score (Coefficient of Determination)

R² provides a percentage-style grade for your model, answering the fundamental question of what proportion of the pattern your model captured. An R² of 1.0 represents perfect predictions, though this should raise suspicion about potential data leakage. A score of 0.85 indicates your model explains 85% of why the target variable varies, with the remaining 15% attributable to uncaptured factors or random noise. An R² of 0.0 means your model performs no better than simply predicting the average value every time, while negative R² values indicate your model is actually worse than this naive baseline — a clear sign something has gone wrong.

The primary advantage of R² is its scale-independence, enabling direct comparison across models predicting vastly different quantities. You can meaningfully compare a model predicting $500K house prices (R² = 0.9) to a model predicting 72°F temperatures (R² = 0.9) because the metric represents explained variance rather than absolute error. This makes R² invaluable for communicating model quality to stakeholders who may not understand domain-specific error units.

Cross-Validation for Regression

A fundamental problem with single train-test splits is their susceptibility to luck. Perhaps by chance, all the easy examples ended up in training while difficult ones landed in testing, or vice versa. Cross-validation addresses this by testing your model multiple times on different data partitions, then averaging the results. In k-fold cross-validation, the data is divided into k equal parts (commonly 5 or 10). The model trains on k-1 parts and tests on the remaining part, repeating this process k times so every data point serves as a test example exactly once.

The averaged cross-validation score provides a much more reliable estimate of real-world performance than any single split. Equally important is the standard deviation across folds — if scores vary wildly from fold to fold, your model may be unstable or overly sensitive to the specific training examples it receives. Always use cross-validation when comparing models or tuning hyperparameters, as a single train-test split simply doesn't provide enough evidence to make confident decisions about model selection.

# Cross-validation for reliable evaluation
from sklearn.model_selection import cross_val_score, KFold
from sklearn.linear_model import LinearRegression, Ridge, Lasso
from sklearn.datasets import make_regression
import numpy as np

# Generate dataset
X, y = make_regression(n_samples=200, n_features=10, noise=20, random_state=42)

# Compare models using cross-validation
models = {
    'Linear Regression': LinearRegression(),
    'Ridge (α=1.0)': Ridge(alpha=1.0),
    'Ridge (α=10.0)': Ridge(alpha=10.0),
    'Lasso (α=0.1)': Lasso(alpha=0.1)
}

print("5-Fold Cross-Validation Results:")
print("-" * 50)
for name, model in models.items():
    # Use negative MSE (sklearn convention) and R²
    mse_scores = -cross_val_score(model, X, y, cv=5, scoring='neg_mean_squared_error')
    r2_scores = cross_val_score(model, X, y, cv=5, scoring='r2')
    
    print(f"{name}:")
    print(f"  RMSE: {np.sqrt(mse_scores.mean()):.2f} ± {np.sqrt(mse_scores.std()):.2f}")
    print(f"  R²:   {r2_scores.mean():.4f} ± {r2_scores.std():.4f}")

Residual Analysis

Residuals (actual - predicted) reveal problems that metrics alone might miss. Plotting residuals against predicted values should show random scatter around zero. Patterns in residuals indicate model issues: a funnel shape suggests heteroscedasticity (non-constant variance), a curve suggests missing non-linear terms, and systematic bias (residuals consistently above or below zero) suggests systematic prediction errors. Always visualize residuals as part of your evaluation workflow.

# Residual Analysis
import matplotlib.pyplot as plt
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split

# Generate data with non-linear relationship
np.random.seed(42)
X = np.random.rand(100, 1) * 10
y = 2 + 3*X.flatten() + 0.5*X.flatten()**2 + np.random.randn(100)*3

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3)

# Fit linear model to non-linear data
model = LinearRegression()
model.fit(X_train, y_train)
y_pred = model.predict(X_test)

# Residual analysis
residuals = y_test - y_pred

plt.figure(figsize=(12, 4))

# Residuals vs Predicted
plt.subplot(1, 3, 1)
plt.scatter(y_pred, residuals, alpha=0.5)
plt.axhline(y=0, color='r', linestyle='--')
plt.xlabel('Predicted Values')
plt.ylabel('Residuals')
plt.title('Residuals vs Predicted')

# Residual distribution
plt.subplot(1, 3, 2)
plt.hist(residuals, bins=20, edgecolor='black')
plt.xlabel('Residual Value')
plt.ylabel('Frequency')
plt.title('Residual Distribution')

# Q-Q plot for normality
plt.subplot(1, 3, 3)
from scipy import stats
stats.probplot(residuals, dist="norm", plot=plt)
plt.title('Q-Q Plot')

plt.tight_layout()
plt.show()