xbe.at

Machine learning / Regression

Mastering Linear Regression for Machine Learning

8 min read, 14 sections

Contents
  1. Linear Regression Fundamentals
  2. Loss Function Implementation
  3. Gradient Descent Optimization
  4. Feature Scaling and Normalization
  5. R-Squared Implementation
  6. Statistical Significance Testing
  7. Cross-Validation Implementation
  8. Regularization Techniques
  9. Multicollinearity Detection
  10. Real-world Example - Housing Price Prediction
  11. Performance Metrics Implementation
  12. Residual Analysis Implementation
  13. Real-world Example - Time Series Regression
  14. Additional Resources

Linear Regression Fundamentals

Linear regression forms the backbone of predictive modeling by establishing relationships between variables through a linear equation. The fundamental concept involves fitting a line that minimizes the distance between predicted and actual values using ordinary least squares optimization.

import numpy as np
import matplotlib.pyplot as plt

class LinearRegression:
    def __init__(self):
        self.weights = None
        self.bias = None
        
    def fit(self, X, y, learning_rate=0.01, epochs=1000):
        n_samples = X.shape[0]
        self.weights = np.zeros(X.shape[1])
        self.bias = 0
        
        # Training loop
        for _ in range(epochs):
            # Forward pass
            y_pred = np.dot(X, self.weights) + self.bias
            
            # Compute gradients
            dw = (1/n_samples) * np.dot(X.T, (y_pred - y))
            db = (1/n_samples) * np.sum(y_pred - y)
            
            # Update parameters
            self.weights -= learning_rate * dw
            self.bias -= learning_rate * db
            
    def predict(self, X):
        return np.dot(X, self.weights) + self.bias

Loss Function Implementation

The loss function quantifies prediction errors during model training. Mean Squared Error (MSE) is commonly used, calculating the average squared difference between predicted and actual values to penalize larger errors more heavily.

def compute_mse_loss(y_true, y_pred):
    """
    Compute Mean Squared Error loss
    Formula: MSE = (1/n) * Σ(y_true - y_pred)²
    """
    n_samples = len(y_true)
    mse = np.sum((y_true - y_pred) ** 2) / n_samples
    return mse

# Example usage
y_true = np.array([1, 2, 3, 4, 5])
y_pred = np.array([1.1, 2.2, 2.9, 4.1, 5.2])
loss = compute_mse_loss(y_true, y_pred)
print(f"MSE Loss: {loss:.4f}")

Gradient Descent Optimization

The gradient descent algorithm iteratively adjusts model parameters by computing the gradient of the loss function with respect to each parameter. This optimization technique guides the model toward the global minimum of the cost function.

def gradient_descent(X, y, learning_rate=0.01, epochs=1000):
    n_samples, n_features = X.shape
    weights = np.zeros(n_features)
    bias = 0
    losses = []
    
    for epoch in range(epochs):
        # Compute predictions
        y_pred = np.dot(X, weights) + bias
        
        # Compute gradients
        dw = (2/n_samples) * np.dot(X.T, (y_pred - y))
        db = (2/n_samples) * np.sum(y_pred - y)
        
        # Update parameters
        weights -= learning_rate * dw
        bias -= learning_rate * db
        
        # Track loss
        loss = compute_mse_loss(y, y_pred)
        losses.append(loss)
        
    return weights, bias, losses

Feature Scaling and Normalization

Feature scaling ensures all variables contribute equally to the model and accelerates convergence. This implementation demonstrates standardization (z-score normalization) and min-max scaling methods for preprocessing features.

class FeatureScaler:
    def __init__(self, method='standardization'):
        self.method = method
        self.params = {}
        
    def fit_transform(self, X):
        if self.method == 'standardization':
            self.params['mean'] = np.mean(X, axis=0)
            self.params['std'] = np.std(X, axis=0)
            return (X - self.params['mean']) / self.params['std']
        
        elif self.method == 'minmax':
            self.params['min'] = np.min(X, axis=0)
            self.params['max'] = np.max(X, axis=0)
            return (X - self.params['min']) / (self.params['max'] - self.params['min'])
    
    def transform(self, X):
        if self.method == 'standardization':
            return (X - self.params['mean']) / self.params['std']
        elif self.method == 'minmax':
            return (X - self.params['min']) / (self.params['max'] - self.params['min'])

R-Squared Implementation

The coefficient of determination (R²) measures the proportion of variance in the dependent variable explained by the independent variables. This implementation calculates R² and provides statistical interpretation.

def calculate_r_squared(y_true, y_pred):
    """
    Calculate R-squared score
    R² = 1 - (Sum of squared residuals / Total sum of squares)
    """
    ss_residual = np.sum((y_true - y_pred) ** 2)
    ss_total = np.sum((y_true - np.mean(y_true)) ** 2)
    r_squared = 1 - (ss_residual / ss_total)
    
    return r_squared

def calculate_adjusted_r_squared(r_squared, n_samples, n_features):
    """
    Calculate Adjusted R-squared
    Adj_R² = 1 - [(1 - R²)(n-1)/(n-k-1)]
    where n is sample size and k is number of features
    """
    adjusted_r_squared = 1 - (1 - r_squared) * ((n_samples - 1) / (n_samples - n_features - 1))
    return adjusted_r_squared

Statistical Significance Testing

This implementation calculates p-values and confidence intervals for regression coefficients using the t-distribution. The statistical testing helps determine the reliability and significance of each predictor variable.

def calculate_statistical_significance(X, y, weights, y_pred):
    """
    Calculate p-values and confidence intervals for regression coefficients
    """
    n_samples, n_features = X.shape
    
    # Calculate standard errors
    mse = np.sum((y - y_pred) ** 2) / (n_samples - n_features)
    var_coef = mse * np.linalg.inv(np.dot(X.T, X)).diagonal()
    std_errors = np.sqrt(var_coef)
    
    # Calculate t-statistics
    t_stats = weights / std_errors
    
    # Calculate p-values
    from scipy import stats
    p_values = 2 * (1 - stats.t.cdf(np.abs(t_stats), n_samples - n_features))
    
    # Calculate 95% confidence intervals
    t_critical = stats.t.ppf(0.975, n_samples - n_features)
    ci_lower = weights - t_critical * std_errors
    ci_upper = weights + t_critical * std_errors
    
    return p_values, (ci_lower, ci_upper)

Cross-Validation Implementation

Cross-validation provides a robust method for assessing model performance and generalization capabilities. This implementation includes k-fold cross-validation with performance metrics calculation.

def k_fold_cross_validation(X, y, k=5):
    """
    Perform k-fold cross-validation for linear regression
    """
    n_samples = len(y)
    fold_size = n_samples // k
    indices = np.arange(n_samples)
    np.random.shuffle(indices)
    
    scores = []
    for i in range(k):
        # Create train/test splits
        test_start = i * fold_size
        test_end = (i + 1) * fold_size
        test_indices = indices[test_start:test_end]
        train_indices = np.concatenate([indices[:test_start], indices[test_end:]])
        
        # Split data
        X_train, X_test = X[train_indices], X[test_indices]
        y_train, y_test = y[train_indices], y[test_indices]
        
        # Train model
        model = LinearRegression()
        model.fit(X_train, y_train)
        
        # Evaluate
        y_pred = model.predict(X_test)
        score = calculate_r_squared(y_test, y_pred)
        scores.append(score)
    
    return np.mean(scores), np.std(scores)

Regularization Techniques

Regularization prevents overfitting by adding penalty terms to the loss function. This implementation includes both L1 (Lasso) and L2 (Ridge) regularization methods for linear regression.

class RegularizedLinearRegression:
    def __init__(self, alpha=1.0, regularization='l2'):
        self.alpha = alpha
        self.regularization = regularization
        self.weights = None
        self.bias = None
    
    def fit(self, X, y, learning_rate=0.01, epochs=1000):
        n_samples, n_features = X.shape
        self.weights = np.zeros(n_features)
        self.bias = 0
        
        for _ in range(epochs):
            y_pred = np.dot(X, self.weights) + self.bias
            
            # Compute gradients with regularization
            if self.regularization == 'l2':
                reg_term = self.alpha * self.weights
            elif self.regularization == 'l1':
                reg_term = self.alpha * np.sign(self.weights)
            
            dw = (1/n_samples) * (np.dot(X.T, (y_pred - y)) + reg_term)
            db = (1/n_samples) * np.sum(y_pred - y)
            
            self.weights -= learning_rate * dw
            self.bias -= learning_rate * db
    
    def predict(self, X):
        return np.dot(X, self.weights) + self.bias

Multicollinearity Detection

Multicollinearity occurs when predictor variables are highly correlated, affecting model stability. This implementation detects multicollinearity using Variance Inflation Factor (VIF) analysis.

def calculate_vif(X):
    """
    Calculate Variance Inflation Factor for each feature
    """
    from sklearn.linear_model import LinearRegression
    vif_dict = {}
    
    for i in range(X.shape[1]):
        # Select all columns except the current one
        X_other = np.delete(X, i, axis=1)
        
        # Regression of feature i on other features
        regressor = LinearRegression()
        regressor.fit(X_other, X[:, i])
        
        # Calculate R² and VIF
        y_pred = regressor.predict(X_other)
        r_squared = calculate_r_squared(X[:, i], y_pred)
        vif = 1 / (1 - r_squared)
        
        vif_dict[f"Feature_{i}"] = vif
    
    return vif_dict

Real-world Example - Housing Price Prediction

# Load and preprocess housing dataset
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split

# Load data
housing = fetch_california_housing()
X, y = housing.data, housing.target

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

# Scale features
scaler = FeatureScaler(method='standardization')
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

# Train model with regularization
model = RegularizedLinearRegression(alpha=0.1, regularization='l2')
model.fit(X_train_scaled, y_train)

# Evaluate model
y_pred = model.predict(X_test_scaled)
r2_score = calculate_r_squared(y_test, y_pred)
print(f"R² Score: {r2_score:.4f}")

# Calculate feature importance
feature_importance = np.abs(model.weights)
for i, importance in enumerate(feature_importance):
    print(f"Feature {housing.feature_names[i]}: {importance:.4f}")

Performance Metrics Implementation

This comprehensive metrics implementation provides a complete evaluation toolkit for regression models, including MSE, RMSE, MAE, and MAPE, enabling thorough model assessment across different error measures.

class RegressionMetrics:
    def __init__(self, y_true, y_pred):
        self.y_true = y_true
        self.y_pred = y_pred
        self.n_samples = len(y_true)
    
    def mean_squared_error(self):
        """Calculate Mean Squared Error"""
        return np.mean((self.y_true - self.y_pred) ** 2)
    
    def root_mean_squared_error(self):
        """Calculate Root Mean Squared Error"""
        return np.sqrt(self.mean_squared_error())
    
    def mean_absolute_error(self):
        """Calculate Mean Absolute Error"""
        return np.mean(np.abs(self.y_true - self.y_pred))
    
    def mean_absolute_percentage_error(self):
        """Calculate Mean Absolute Percentage Error"""
        return np.mean(np.abs((self.y_true - self.y_pred) / self.y_true)) * 100
    
    def get_all_metrics(self):
        return {
            'MSE': self.mean_squared_error(),
            'RMSE': self.root_mean_squared_error(),
            'MAE': self.mean_absolute_error(),
            'MAPE': self.mean_absolute_percentage_error()
        }

Residual Analysis Implementation

Residual analysis helps verify regression assumptions and identify potential model issues. This implementation provides tools for analyzing residual patterns and detecting heteroscedasticity.

class ResidualAnalyzer:
    def __init__(self, y_true, y_pred):
        self.residuals = y_true - y_pred
        self.standardized_residuals = self.residuals / np.std(self.residuals)
        
    def plot_residuals(self):
        plt.figure(figsize=(10, 6))
        plt.scatter(y_pred, self.standardized_residuals)
        plt.axhline(y=0, color='r', linestyle='--')
        plt.xlabel('Predicted Values')
        plt.ylabel('Standardized Residuals')
        plt.title('Residual Plot')
        plt.show()
        
    def test_normality(self):
        """Test residuals for normality using Shapiro-Wilk test"""
        from scipy import stats
        statistic, p_value = stats.shapiro(self.residuals)
        return {'statistic': statistic, 'p_value': p_value}
    
    def test_heteroscedasticity(self):
        """Breusch-Pagan test for heteroscedasticity"""
        squared_residuals = self.residuals ** 2
        model = LinearRegression()
        model.fit(y_pred.reshape(-1, 1), squared_residuals)
        
        n = len(self.residuals)
        r_squared = calculate_r_squared(squared_residuals, model.predict(y_pred.reshape(-1, 1)))
        lm_stat = n * r_squared
        p_value = 1 - stats.chi2.cdf(lm_stat, df=1)
        
        return {'statistic': lm_stat, 'p_value': p_value}

Real-world Example - Time Series Regression

import pandas as pd
from datetime import datetime, timedelta

# Generate synthetic time series data
dates = pd.date_range(start='2023-01-01', end='2023-12-31', freq='D')
n_samples = len(dates)

# Create features
X = np.column_stack([
    np.sin(2 * np.pi * np.arange(n_samples) / 365),  # Yearly seasonality
    np.sin(2 * np.pi * np.arange(n_samples) / 7),    # Weekly seasonality
    np.arange(n_samples)                             # Trend
])

# Generate target variable with noise
y = 10 + 2 * X[:, 0] + 1.5 * X[:, 1] + 0.01 * X[:, 2] + np.random.normal(0, 0.5, n_samples)

# Train model
model = LinearRegression()
model.fit(X[:-30], y[:-30])  # Train on all but last 30 days

# Make predictions
y_pred = model.predict(X[-30:])  # Predict last 30 days

# Evaluate predictions
metrics = RegressionMetrics(y[-30:], y_pred)
results = metrics.get_all_metrics()
print("Forecast Evaluation Metrics:")
for metric, value in results.items():
    print(f"{metric}: {value:.4f}")

Additional Resources

More in Regression

  1. Poisson Regression vs. Linear Regression in Python11 min

    Regression analysis is a statistical method used to model the relationship between variables. We'll explore two popular types: Linear Regression and Poisson...

  2. Logistic Regression Fundamentals From Scratch in Python14 min

    Logistic Regression is a versatile statistical method used for predicting binary outcomes. It's commonly applied in various fields:

  3. Predicting Binary Outcomes with Logistic Regression10 min

    Logistic regression is a powerful statistical method used for predicting binary outcomes. It's particularly useful when we need to classify data into two...

  4. Linear Regression with Multiple Variables in Python3 min

    Linear regression with multiple variables is an extension of simple linear regression, where we predict a target variable (y) based on multiple input...

  5. Predictive Modeling with AutoML in Python6 min

    H2O AutoML is a powerful automated machine learning library that automates the process of building and comparing multiple machine learning models. It handles...

  6. Linear Regression Implementation from Scratch10 min

    Linear regression forms the foundation of predictive modeling by establishing relationships between variables through a linear equation. This implementation...