xbe.at

Machine learning / Regression

Exploring Linear Regression for Predictive Modeling

12 min read, 15 sections

Contents
  1. Understanding Linear Regression Fundamentals
  2. Mathematical Foundation of Linear Regression
  3. Implementing Linear Regression from Scratch
  4. Assessing Model Performance
  5. Validating Regression Assumptions
  6. Multiple Linear Regression Implementation
  7. Feature Engineering and Polynomial Regression
  8. Cross-Validation and Model Selection
  9. Regularization Techniques
  10. Real-World Application - Housing Price Prediction
  11. Real-World Application - Stock Price Analysis
  12. Results Analysis and Visualization
  13. Advanced Regularization Techniques
  14. Diagnostic Tools and Model Validation
  15. Additional Resources

Understanding Linear Regression Fundamentals

Linear regression forms the backbone of predictive modeling by establishing relationships between variables through a linear equation. The fundamental concept involves finding the best-fitting line that minimizes the sum of squared residuals between predicted and actual values.

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression

# Generate sample data
np.random.seed(42)
X = np.random.randn(100, 1) * 10
y = 2 * X + 1 + np.random.randn(100, 1) * 2

# Create and fit the model
model = LinearRegression()
model.fit(X, y)

# Print model parameters
print(f"Slope: {model.coef_[0][0]:.2f}")
print(f"Intercept: {model.intercept_[0]:.2f}")

# Example output:
# Slope: 1.98
# Intercept: 0.97

Mathematical Foundation of Linear Regression

The mathematical foundation of linear regression relies on the method of least squares, which minimizes the sum of squared differences between observed and predicted values, forming the basis for parameter estimation.

# Mathematical representation of Linear Regression
"""
Simple Linear Regression Equation:
$$y = \beta_0 + \beta_1x + \epsilon$$

Cost Function (Mean Squared Error):
$$MSE = \frac{1}{n}\sum_{i=1}^{n}(y_i - (\beta_0 + \beta_1x_i))^2$$

Parameter Estimation:
$$\beta_1 = \frac{\sum_{i=1}^{n}(x_i - \bar{x})(y_i - \bar{y})}{\sum_{i=1}^{n}(x_i - \bar{x})^2}$$
$$\beta_0 = \bar{y} - \beta_1\bar{x}$$
"""

Implementing Linear Regression from Scratch

Understanding the inner workings of linear regression requires implementing it from scratch. This implementation demonstrates the core computational aspects without relying on external libraries for the regression calculations.

class SimpleLinearRegression:
    def fit(self, X, y):
        # Calculate means
        X_mean = np.mean(X)
        y_mean = np.mean(y)
        
        # Calculate β1 (slope)
        numerator = np.sum((X - X_mean) * (y - y_mean))
        denominator = np.sum((X - X_mean) ** 2)
        self.beta1 = numerator / denominator
        
        # Calculate β0 (intercept)
        self.beta0 = y_mean - self.beta1 * X_mean
        
    def predict(self, X):
        return self.beta0 + self.beta1 * X

# Example usage
model = SimpleLinearRegression()
model.fit(X.flatten(), y.flatten())
print(f"Slope: {model.beta1:.2f}")
print(f"Intercept: {model.beta0:.2f}")

Assessing Model Performance

Model evaluation is crucial for understanding the reliability and accuracy of linear regression predictions. We implement various metrics to quantify model performance and validate our assumptions.

def evaluate_model(y_true, y_pred):
    # Calculate R-squared
    ss_tot = np.sum((y_true - np.mean(y_true)) ** 2)
    ss_res = np.sum((y_true - y_pred) ** 2)
    r2 = 1 - (ss_res / ss_tot)
    
    # Calculate Mean Squared Error
    mse = np.mean((y_true - y_pred) ** 2)
    
    # Calculate Root Mean Squared Error
    rmse = np.sqrt(mse)
    
    # Calculate Mean Absolute Error
    mae = np.mean(np.abs(y_true - y_pred))
    
    return {
        'R2': r2,
        'MSE': mse,
        'RMSE': rmse,
        'MAE': mae
    }

# Example usage
y_pred = model.predict(X.flatten())
metrics = evaluate_model(y.flatten(), y_pred)
for metric, value in metrics.items():
    print(f"{metric}: {value:.4f}")

Validating Regression Assumptions

Statistical validation of linear regression assumptions is essential for ensuring model reliability. This implementation provides methods to check linearity, normality, and homoscedasticity assumptions.

def check_assumptions(X, y, y_pred):
    import scipy.stats as stats
    
    # Calculate residuals
    residuals = y - y_pred
    
    # Normality test (Shapiro-Wilk)
    _, normality_p = stats.shapiro(residuals)
    
    # Homoscedasticity (Breusch-Pagan test)
    resi_sq = residuals ** 2
    _, homo_p = stats.pearsonr(X.flatten(), resi_sq.flatten())
    
    # Linearity test (Ramsey RESET test)
    y_pred_sq = y_pred ** 2
    X_extended = np.column_stack([X, y_pred_sq])
    
    results = {
        'Normality p-value': normality_p,
        'Homoscedasticity p-value': homo_p,
        'Mean residual': np.mean(residuals)
    }
    
    return results

# Example usage
assumption_tests = check_assumptions(X, y, y_pred)
for test, value in assumption_tests.items():
    print(f"{test}: {value:.4f}")

Multiple Linear Regression Implementation

Multiple linear regression extends the simple linear model by incorporating multiple independent variables. This implementation demonstrates how to handle multiple predictors while maintaining computational efficiency.

class MultipleLinearRegression:
    def fit(self, X, y):
        # Add column of ones for intercept
        X_b = np.column_stack([np.ones(len(X)), X])
        
        # Calculate parameters using normal equation
        # β = (X^T X)^(-1) X^T y
        self.theta = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)
        
        self.intercept_ = self.theta[0]
        self.coef_ = self.theta[1:]
    
    def predict(self, X):
        X_b = np.column_stack([np.ones(len(X)), X])
        return X_b.dot(self.theta)

# Example usage with multiple features
X_multi = np.random.randn(100, 3)
y_multi = 2 * X_multi[:, 0] + 1.5 * X_multi[:, 1] - 0.5 * X_multi[:, 2] + 2 + np.random.randn(100) * 0.1

model_multi = MultipleLinearRegression()
model_multi.fit(X_multi, y_multi)
print(f"Intercept: {model_multi.intercept_:.2f}")
print(f"Coefficients: {model_multi.coef_}")

Feature Engineering and Polynomial Regression

Polynomial regression extends linear regression by incorporating higher-order terms, enabling the modeling of non-linear relationships while maintaining the linear regression framework.

def create_polynomial_features(X, degree=2):
    n_samples, n_features = X.shape
    def combinations_with_replacement(n, r):
        from itertools import combinations_with_replacement
        return list(combinations_with_replacement(range(n), r))
    
    combinations = []
    for d in range(1, degree + 1):
        combinations.extend(combinations_with_replacement(n_features, d))
    
    n_output_features = len(combinations)
    X_poly = np.empty((n_samples, n_output_features))
    
    for i, combination in enumerate(combinations):
        X_poly[:, i] = np.prod(X[:, combination], axis=1)
    
    return X_poly

# Example usage
X_poly = create_polynomial_features(X_multi, degree=2)
model_poly = MultipleLinearRegression()
model_poly.fit(X_poly, y_multi)
print(f"Number of polynomial features: {X_poly.shape[1]}")
print(f"First few coefficients: {model_poly.coef_[:3]}")

Cross-Validation and Model Selection

Cross-validation provides a robust method for assessing model performance and selecting optimal hyperparameters. This implementation demonstrates k-fold cross-validation for linear regression.

def cross_validate(X, y, k_folds=5):
    from sklearn.model_selection import KFold
    kf = KFold(n_splits=k_folds, shuffle=True, random_state=42)
    
    scores = {
        'r2_scores': [],
        'mse_scores': []
    }
    
    for train_idx, test_idx in kf.split(X):
        X_train, X_test = X[train_idx], X[test_idx]
        y_train, y_test = y[train_idx], y[test_idx]
        
        model = MultipleLinearRegression()
        model.fit(X_train, y_train)
        y_pred = model.predict(X_test)
        
        # Calculate metrics
        metrics = evaluate_model(y_test, y_pred)
        scores['r2_scores'].append(metrics['R2'])
        scores['mse_scores'].append(metrics['MSE'])
    
    return {
        'mean_r2': np.mean(scores['r2_scores']),
        'std_r2': np.std(scores['r2_scores']),
        'mean_mse': np.mean(scores['mse_scores']),
        'std_mse': np.std(scores['mse_scores'])
    }

# Example usage
cv_results = cross_validate(X_multi, y_multi)
for metric, value in cv_results.items():
    print(f"{metric}: {value:.4f}")

Regularization Techniques

Regularization helps prevent overfitting by adding penalty terms to the cost function. This implementation includes Ridge (L2) and Lasso (L1) regularization options.

class RegularizedLinearRegression:
    def __init__(self, alpha=1.0, regularization='ridge'):
        self.alpha = alpha
        self.regularization = regularization
    
    def fit(self, X, y):
        X_b = np.column_stack([np.ones(len(X)), X])
        n_features = X_b.shape[1]
        
        if self.regularization == 'ridge':
            # Ridge regression (L2)
            I = np.eye(n_features)
            I[0, 0] = 0  # Don't regularize intercept
            self.theta = np.linalg.inv(X_b.T.dot(X_b) + 
                                     self.alpha * I).dot(X_b.T).dot(y)
        
        elif self.regularization == 'lasso':
            # Lasso regression (L1) using coordinate descent
            self.theta = np.zeros(n_features)
            max_iter = 1000
            for _ in range(max_iter):
                theta_old = self.theta.copy()
                for j in range(n_features):
                    if j == 0:  # Don't regularize intercept
                        self.theta[j] = np.mean(y - X_b[:, 1:].dot(self.theta[1:]))
                    else:
                        r = y - X_b.dot(self.theta) + X_b[:, j] * self.theta[j]
                        self.theta[j] = np.sign(X_b[:, j].dot(r)) * \
                                      max(0, abs(X_b[:, j].dot(r)) - self.alpha) / \
                                      (X_b[:, j].dot(X_b[:, j]))
                if np.allclose(self.theta, theta_old):
                    break
        
        self.intercept_ = self.theta[0]
        self.coef_ = self.theta[1:]
    
    def predict(self, X):
        X_b = np.column_stack([np.ones(len(X)), X])
        return X_b.dot(self.theta)

# Example usage
ridge_model = RegularizedLinearRegression(alpha=1.0, regularization='ridge')
ridge_model.fit(X_multi, y_multi)
print(f"Ridge coefficients: {ridge_model.coef_}")

lasso_model = RegularizedLinearRegression(alpha=0.1, regularization='lasso')
lasso_model.fit(X_multi, y_multi)
print(f"Lasso coefficients: {lasso_model.coef_}")

Real-World Application - Housing Price Prediction

This implementation demonstrates a complete workflow for predicting housing prices using multiple features, including data preprocessing, feature selection, and model evaluation.

import pandas as pd
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split

# Generate synthetic housing data
np.random.seed(42)
n_samples = 1000

data = {
    'size_sqft': np.random.normal(2000, 500, n_samples),
    'bedrooms': np.random.randint(1, 6, n_samples),
    'bathrooms': np.random.randint(1, 4, n_samples),
    'age_years': np.random.normal(20, 10, n_samples),
    'distance_downtown': np.random.normal(5, 2, n_samples),
    'price': None
}

# Create price with realistic relationships
data['price'] = (
    200000 + 
    150 * data['size_sqft'] + 
    25000 * data['bedrooms'] +
    35000 * data['bathrooms'] -
    2000 * data['age_years'] -
    15000 * data['distance_downtown'] +
    np.random.normal(0, 25000, n_samples)
)

# Create DataFrame
df = pd.DataFrame(data)

# Preprocess data
X = df.drop('price', axis=1)
y = df['price']

# Scale features
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

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

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

# Evaluate model
y_pred = model.predict(X_test)
metrics = evaluate_model(y_test, y_pred)
print("\nModel Performance Metrics:")
for metric, value in metrics.items():
    print(f"{metric}: {value:.4f}")

# Feature importance
feature_importance = pd.DataFrame({
    'Feature': X.columns,
    'Coefficient': model.coef_
}).sort_values('Coefficient', key=abs, ascending=False)

print("\nFeature Importance:")
print(feature_importance)

Real-World Application - Stock Price Analysis

Implementation of linear regression for analyzing stock price trends and making predictions based on historical data and technical indicators.

def create_technical_indicators(prices, window=14):
    df = pd.DataFrame(prices, columns=['close'])
    
    # Calculate Moving Averages
    df['MA20'] = df['close'].rolling(window=20).mean()
    df['MA50'] = df['close'].rolling(window=50).mean()
    
    # Calculate RSI
    delta = df['close'].diff()
    gain = (delta.where(delta > 0, 0)).rolling(window=window).mean()
    loss = (-delta.where(delta < 0, 0)).rolling(window=window).mean()
    rs = gain / loss
    df['RSI'] = 100 - (100 / (1 + rs))
    
    # Calculate volatility
    df['volatility'] = df['close'].rolling(window=window).std()
    
    # Create target variable (next day's return)
    df['target'] = df['close'].shift(-1) / df['close'] - 1
    
    return df.dropna()

# Generate synthetic stock data
np.random.seed(42)
dates = pd.date_range(start='2020-01-01', end='2023-12-31', freq='D')
prices = np.random.normal(0.001, 0.02, len(dates)).cumsum()
prices = np.exp(prices) * 100  # Start price at 100

# Prepare data
df = create_technical_indicators(prices)
features = ['MA20', 'MA50', 'RSI', 'volatility']
X = df[features]
y = df['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
)

scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

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

# Evaluate predictions
y_pred = model.predict(X_test_scaled)
metrics = evaluate_model(y_test, y_pred)

print("\nStock Price Prediction Metrics:")
for metric, value in metrics.items():
    print(f"{metric}: {value:.4f}")

# Feature importance
feature_importance = pd.DataFrame({
    'Feature': features,
    'Coefficient': model.coef_
}).sort_values('Coefficient', key=abs, ascending=False)

print("\nFeature Importance in Stock Prediction:")
print(feature_importance)

Results Analysis and Visualization

Comprehensive visualization and analysis of regression results, including residual plots, Q-Q plots, and feature importance visualization.

def visualize_regression_diagnostics(y_true, y_pred, X, feature_names):
    fig, axes = plt.subplots(2, 2, figsize=(15, 12))
    
    # Residual plot
    residuals = y_true - y_pred
    axes[0, 0].scatter(y_pred, residuals, alpha=0.5)
    axes[0, 0].axhline(y=0, color='r', linestyle='--')
    axes[0, 0].set_xlabel('Predicted values')
    axes[0, 0].set_ylabel('Residuals')
    axes[0, 0].set_title('Residual Plot')
    
    # Q-Q plot
    from scipy import stats
    stats.probplot(residuals, dist="norm", plot=axes[0, 1])
    axes[0, 1].set_title('Q-Q Plot')
    
    # Feature importance plot
    model = RegularizedLinearRegression(alpha=0.1)
    model.fit(X, y_true)
    importance = np.abs(model.coef_)
    importance = importance / np.sum(importance)
    
    axes[1, 0].barh(feature_names, importance)
    axes[1, 0].set_title('Feature Importance')
    
    # Prediction vs Actual
    axes[1, 1].scatter(y_true, y_pred, alpha=0.5)
    axes[1, 1].plot([y_true.min(), y_true.max()], 
                    [y_true.min(), y_true.max()], 
                    'r--', lw=2)
    axes[1, 1].set_xlabel('Actual values')
    axes[1, 1].set_ylabel('Predicted values')
    axes[1, 1].set_title('Prediction vs Actual')
    
    plt.tight_layout()
    return fig

# Example usage
fig = visualize_regression_diagnostics(
    y_test, y_pred, X_test_scaled, features
)
plt.show()

Advanced Regularization Techniques

Advanced regularization methods incorporating Elastic Net and adaptive penalties, providing flexible control over feature selection and model complexity while maintaining predictive performance.

class ElasticNetRegression:
    def __init__(self, alpha=1.0, l1_ratio=0.5, max_iter=1000):
        self.alpha = alpha
        self.l1_ratio = l1_ratio
        self.max_iter = max_iter
    
    def fit(self, X, y):
        X_b = np.column_stack([np.ones(len(X)), X])
        n_samples, n_features = X_b.shape
        self.theta = np.zeros(n_features)
        
        for _ in range(self.max_iter):
            theta_old = self.theta.copy()
            
            for j in range(n_features):
                if j == 0:  # Don't regularize intercept
                    self.theta[j] = np.mean(y - X_b[:, 1:].dot(self.theta[1:]))
                else:
                    # Elastic Net update
                    r = y - X_b.dot(self.theta) + X_b[:, j] * self.theta[j]
                    rho = X_b[:, j].dot(X_b[:, j])
                    
                    # Soft thresholding with combined L1 and L2 penalties
                    if self.l1_ratio > 0:
                        soft_threshold = np.sign(X_b[:, j].dot(r)) * \
                            max(0, abs(X_b[:, j].dot(r)) - self.alpha * self.l1_ratio) / \
                            (rho + self.alpha * (1 - self.l1_ratio))
                        self.theta[j] = soft_threshold
                    else:
                        self.theta[j] = X_b[:, j].dot(r) / (rho + self.alpha)
            
            if np.allclose(self.theta, theta_old, rtol=1e-4):
                break
        
        self.intercept_ = self.theta[0]
        self.coef_ = self.theta[1:]
    
    def predict(self, X):
        X_b = np.column_stack([np.ones(len(X)), X])
        return X_b.dot(self.theta)

# Example usage with comparison
np.random.seed(42)
n_samples = 200
n_features = 10

# Generate sparse coefficients
true_coef = np.zeros(n_features)
true_coef[0:3] = [1.5, -2, 3]

# Generate data with noise
X = np.random.randn(n_samples, n_features)
y = X.dot(true_coef) + np.random.randn(n_samples) * 0.1

# Train and compare models
models = {
    'Ridge': RegularizedLinearRegression(alpha=1.0, regularization='ridge'),
    'Lasso': RegularizedLinearRegression(alpha=0.1, regularization='lasso'),
    'ElasticNet': ElasticNetRegression(alpha=0.1, l1_ratio=0.5)
}

results = {}
for name, model in models.items():
    model.fit(X, y)
    y_pred = model.predict(X)
    results[name] = {
        'coefficients': model.coef_,
        'metrics': evaluate_model(y, y_pred)
    }

# Print comparison results
for name, result in results.items():
    print(f"\n{name} Results:")
    print("Coefficients:", result['coefficients'])
    print("Metrics:", result['metrics'])

Diagnostic Tools and Model Validation

Advanced diagnostic tools for comprehensive model validation, including influence measures, outlier detection, and statistical tests for assumption verification.

class RegressionDiagnostics:
    def __init__(self, model, X, y):
        self.model = model
        self.X = X
        self.y = y
        self.n_samples = len(y)
        self.n_features = X.shape[1]
        
        # Fit model and get predictions
        self.y_pred = model.predict(X)
        self.residuals = y - self.y_pred
        
    def calculate_leverage(self):
        X_b = np.column_stack([np.ones(self.n_samples), self.X])
        H = X_b.dot(np.linalg.inv(X_b.T.dot(X_b))).dot(X_b.T)
        return np.diag(H)
    
    def calculate_cooks_distance(self):
        leverage = self.calculate_leverage()
        standardized_residuals = self.residuals / np.sqrt(
            np.sum(self.residuals**2) / (self.n_samples - self.n_features - 1)
        )
        return (standardized_residuals**2 * leverage) / \
               (self.n_features * (1 - leverage))
    
    def test_assumptions(self):
        from scipy import stats
        
        # Normality test
        _, normality_p = stats.shapiro(self.residuals)
        
        # Homoscedasticity test (Breusch-Pagan)
        _, hetero_p = stats.levene(self.y_pred, self.residuals)
        
        # Durbin-Watson test for autocorrelation
        dw_stat = np.sum(np.diff(self.residuals)**2) / np.sum(self.residuals**2)
        
        return {
            'normality_p_value': normality_p,
            'heteroscedasticity_p_value': hetero_p,
            'durbin_watson': dw_stat
        }
    
    def get_influential_points(self, threshold=0.1):
        cooks_d = self.calculate_cooks_distance()
        return np.where(cooks_d > threshold)[0]

# Example usage
diagnostics = RegressionDiagnostics(models['Ridge'], X, y)

# Print diagnostic results
print("\nModel Diagnostics:")
assumptions = diagnostics.test_assumptions()
for test, value in assumptions.items():
    print(f"{test}: {value:.4f}")

influential_points = diagnostics.get_influential_points()
print(f"\nNumber of influential points: {len(influential_points)}")
print(f"Influential point indices: {influential_points}")

Additional Resources

  • Scientific Papers:
  • "A Unified Approach to Interpreting Model Predictions" - https://arxiv.org/abs/1705.07874
  • "Understanding Regularized Linear Models" - https://arxiv.org/abs/2003.12081
  • "Robust Regression and Outlier Detection" - https://arxiv.org/abs/2007.15975
  • "Modern Linear Regression Techniques" - https://arxiv.org/abs/1509.09169
  • Recommended Search Topics:
  • Advanced regularization techniques in linear regression
  • Robust regression methods for outlier detection
  • Modern approaches to feature selection in linear models
  • Statistical learning theory and model validation
  • Online Resources:
  • Elements of Statistical Learning (Stanford)
  • Journal of Machine Learning Research
  • Scikit-learn Documentation for Linear Models

More in Regression

  1. Exploring Multiple Linear Regression and Its Applications10 min

    Multiple Linear Regression (MLR) is a statistical technique used to model the relationship between multiple independent variables and a single dependent...

  2. Comprehensive Linear Regression11 min

    Linear regression is a fundamental machine learning technique used to model the relationship between variables. Imagine you're trying to predict the price of...

  3. Fundamental Assumptions of Linear Regression in Python12 min

    Linear regression is a fundamental statistical technique used for modeling the relationship between a dependent variable and one or more independent...

  4. Catastrophic Forgetting in Linear Regression Models10 min

    Catastrophic Forgetting in Linear Regression

  5. Introducing Simple Linear Regression12 min

    Simple Linear Regression predicts a dependent variable using one independent variable. It's ideal for straightforward relationships that can be represented...

  6. Building a Logistic Regression Algorithm from Scratch in Python6 min

    Logistic regression is a fundamental machine learning algorithm used for binary classification problems. It's a statistical method that estimates the...