linear_regression.hpp

Go to the documentation of this file.

 
#pragma once
 
#include "statcpp/basic_statistics.hpp"
#include "statcpp/continuous_distributions.hpp"
#include "statcpp/correlation_covariance.hpp"
 
#include <cmath>
#include <cstddef>
#include <limits>
#include <stdexcept>
#include <string>
#include <vector>
 
namespace statcpp {
 
// ============================================================================
// Linear Regression Result Structures
// ============================================================================
 
struct simple_regression_result {
    double intercept;           
    double slope;               
    double intercept_se;        
    double slope_se;            
    double intercept_t;         
    double slope_t;             
    double intercept_p;         
    double slope_p;             
    double r_squared;           
    double adj_r_squared;       
    double residual_se;         
    double f_statistic;         
    double f_p_value;           
    double df_regression;       
    double df_residual;         
    double ss_total;            
    double ss_regression;       
    double ss_residual;         
};
 
struct multiple_regression_result {
    std::vector<double> coefficients;       
    std::vector<double> coefficient_se;     
    std::vector<double> t_statistics;       
    std::vector<double> p_values;           
    double r_squared;                       
    double adj_r_squared;                   
    double residual_se;                     
    double f_statistic;                     
    double f_p_value;                       
    double df_regression;                   
    double df_residual;                     
    double ss_total;                        
    double ss_regression;                   
    double ss_residual;                     
};
 
struct prediction_interval {
    double prediction;      
    double lower;           
    double upper;           
    double se_prediction;   
};
 
struct residual_diagnostics {
    std::vector<double> residuals;              
    std::vector<double> standardized_residuals; 
    std::vector<double> studentized_residuals;  
    std::vector<double> hat_values;             
    std::vector<double> cooks_distance;         
    double durbin_watson;                       
};
 
// ============================================================================
// Simple Linear Regression
// ============================================================================
 
template <typename IteratorX, typename IteratorY>
simple_regression_result simple_linear_regression(IteratorX x_first, IteratorX x_last,
                                                   IteratorY y_first, IteratorY y_last)
{
    auto n_x = statcpp::count(x_first, x_last);
    auto n_y = statcpp::count(y_first, y_last);
 
    if (n_x != n_y) {
        throw std::invalid_argument("statcpp::simple_linear_regression: x and y must have same length");
    }
    if (n_x < 3) {
        throw std::invalid_argument("statcpp::simple_linear_regression: need at least 3 observations");
    }
 
    std::size_t n = n_x;
    double n_d = static_cast<double>(n);
 
    // Calculate means
    double mean_x = statcpp::mean(x_first, x_last);
    double mean_y = statcpp::mean(y_first, y_last);
 
    // Calculate Sxx, Syy, Sxy
    double sxx = 0.0, syy = 0.0, sxy = 0.0;
    auto it_x = x_first;
    auto it_y = y_first;
    for (; it_x != x_last; ++it_x, ++it_y) {
        double dx = static_cast<double>(*it_x) - mean_x;
        double dy = static_cast<double>(*it_y) - mean_y;
        sxx += dx * dx;
        syy += dy * dy;
        sxy += dx * dy;
    }
 
    if (sxx == 0.0) {
        throw std::invalid_argument("statcpp::simple_linear_regression: zero variance in x");
    }
 
    // Regression coefficients
    double slope = sxy / sxx;
    double intercept = mean_y - slope * mean_x;
 
    // Sum of squares
    double ss_total = syy;
    double ss_regression = slope * sxy;
    double ss_residual = std::max(0.0, ss_total - ss_regression);
 
    if (ss_total == 0.0) {
        throw std::invalid_argument("statcpp::simple_linear_regression: zero variance in y (constant response)");
    }
 
    // Residual variance and standard error
    double df_reg = 1.0;
    double df_res = n_d - 2.0;
    double mse = ss_residual / df_res;
    double residual_se = std::sqrt(mse);
 
    // Standard errors of coefficients
    double slope_se = residual_se / std::sqrt(sxx);
    double intercept_se = residual_se * std::sqrt(1.0 / n_d + mean_x * mean_x / sxx);
 
    // t-statistics and p-values
    double slope_t = slope / slope_se;
    double intercept_t = intercept / intercept_se;
    double slope_p = 2.0 * (1.0 - t_cdf(std::abs(slope_t), df_res));
    double intercept_p = 2.0 * (1.0 - t_cdf(std::abs(intercept_t), df_res));
 
    // Coefficient of determination
    double r_squared = ss_regression / ss_total;
    double adj_r_squared = 1.0 - (1.0 - r_squared) * (n_d - 1.0) / df_res;
 
    // F-statistic
    double f_statistic = (ss_regression / df_reg) / mse;
    double f_p_value = 1.0 - f_cdf(f_statistic, df_reg, df_res);
 
    return {
        intercept, slope,
        intercept_se, slope_se,
        intercept_t, slope_t,
        intercept_p, slope_p,
        r_squared, adj_r_squared,
        residual_se,
        f_statistic, f_p_value,
        df_reg, df_res,
        ss_total, ss_regression, ss_residual
    };
}
 
// ============================================================================
// Multiple Linear Regression
// ============================================================================
 
namespace detail {
 
inline void validate_matrix_structure(const std::vector<std::vector<double>>& data,
                                      const char* func_name)
{
    if (data.empty()) {
        std::string msg = "statcpp::";
        msg += func_name;
        msg += ": empty data";
        throw std::invalid_argument(msg);
    }
 
    std::size_t p = data[0].size();
    if (p == 0) {
        std::string msg = "statcpp::";
        msg += func_name;
        msg += ": first row is empty (0 columns)";
        throw std::invalid_argument(msg);
    }
 
    // Verify all rows have the same number of columns
    for (std::size_t i = 1; i < data.size(); ++i) {
        if (data[i].size() != p) {
            std::string msg = "statcpp::";
            msg += func_name;
            msg += ": inconsistent row dimensions (row 0 has ";
            msg += std::to_string(p);
            msg += " columns, but row ";
            msg += std::to_string(i);
            msg += " has ";
            msg += std::to_string(data[i].size());
            msg += " columns)";
            throw std::invalid_argument(msg);
        }
    }
}
 
inline void validate_no_intercept_column(const std::vector<std::vector<double>>& X, const char* func_name)
{
    if (X.empty()) {
        return;
    }
 
    std::size_t n = X.size();
    std::size_t p = X[0].size();
 
    // Check if first column is all 1.0
    if (p > 0) {
        bool all_ones = true;
        for (std::size_t i = 0; i < n; ++i) {
            if (X[i].size() == 0 || std::abs(X[i][0] - 1.0) > 1e-10) {
                all_ones = false;
                break;
            }
        }
 
        if (all_ones) {
            std::string msg = "statcpp::";
            msg += func_name;
            msg += ": X should not contain intercept column (all 1s in first column detected). ";
            msg += "The intercept is added automatically.";
            throw std::invalid_argument(msg);
        }
    }
}
 
// Simple matrix operations (for small matrices)
 
inline std::vector<std::vector<double>> transpose(const std::vector<std::vector<double>>& A)
{
    if (A.empty()) return {};
 
    // Validate matrix structure (all rows have same size)
    std::size_t rows = A.size();
    std::size_t cols = A[0].size();
    for (std::size_t i = 1; i < rows; ++i) {
        if (A[i].size() != cols) {
            throw std::invalid_argument("statcpp::detail::transpose: inconsistent row dimensions");
        }
    }
 
    std::vector<std::vector<double>> result(cols, std::vector<double>(rows));
    for (std::size_t i = 0; i < rows; ++i) {
        for (std::size_t j = 0; j < cols; ++j) {
            result[j][i] = A[i][j];
        }
    }
    return result;
}
 
inline std::vector<std::vector<double>> matrix_multiply(
    const std::vector<std::vector<double>>& A,
    const std::vector<std::vector<double>>& B)
{
    if (A.empty() || B.empty()) return {};
    std::size_t m = A.size();
    std::size_t n = A[0].size();
    std::size_t p = B[0].size();
 
    if (n != B.size()) {
        throw std::invalid_argument("statcpp::detail::matrix_multiply: incompatible dimensions");
    }
 
    // Check that each row has consistent size
    for (std::size_t i = 0; i < m; ++i) {
        if (A[i].size() != n) {
            throw std::invalid_argument("statcpp::detail::matrix_multiply: matrix A has inconsistent row dimensions");
        }
    }
    for (std::size_t k = 0; k < n; ++k) {
        if (B[k].size() != p) {
            throw std::invalid_argument("statcpp::detail::matrix_multiply: matrix B has inconsistent row dimensions");
        }
    }
 
    std::vector<std::vector<double>> result(m, std::vector<double>(p, 0.0));
    for (std::size_t i = 0; i < m; ++i) {
        for (std::size_t j = 0; j < p; ++j) {
            for (std::size_t k = 0; k < n; ++k) {
                result[i][j] += A[i][k] * B[k][j];
            }
        }
    }
    return result;
}
 
inline std::vector<double> matrix_vector_multiply(
    const std::vector<std::vector<double>>& A,
    const std::vector<double>& v)
{
    if (A.empty()) return {};
    std::size_t m = A.size();
    std::size_t n = A[0].size();
 
    if (n != v.size()) {
        throw std::invalid_argument("statcpp::detail::matrix_vector_multiply: incompatible dimensions");
    }
 
    std::vector<double> result(m, 0.0);
    for (std::size_t i = 0; i < m; ++i) {
        for (std::size_t j = 0; j < n; ++j) {
            result[i] += A[i][j] * v[j];
        }
    }
    return result;
}
 
inline std::vector<std::vector<double>> cholesky(const std::vector<std::vector<double>>& A)
{
    std::size_t n = A.size();
    std::vector<std::vector<double>> L(n, std::vector<double>(n, 0.0));
 
    for (std::size_t i = 0; i < n; ++i) {
        for (std::size_t j = 0; j <= i; ++j) {
            double sum = 0.0;
            for (std::size_t k = 0; k < j; ++k) {
                sum += L[i][k] * L[j][k];
            }
            if (i == j) {
                double val = A[i][i] - sum;
                if (val <= 0.0) {
                    throw std::runtime_error("statcpp::detail::cholesky: matrix is not positive definite");
                }
                L[i][j] = std::sqrt(val);
            } else {
                L[i][j] = (A[i][j] - sum) / L[j][j];
            }
        }
    }
    return L;
}
 
inline std::vector<double> solve_cholesky(
    const std::vector<std::vector<double>>& L,
    const std::vector<double>& b)
{
    std::size_t n = L.size();
 
    // Forward substitution: L * y = b
    std::vector<double> y(n);
    for (std::size_t i = 0; i < n; ++i) {
        double sum = 0.0;
        for (std::size_t j = 0; j < i; ++j) {
            sum += L[i][j] * y[j];
        }
        y[i] = (b[i] - sum) / L[i][i];
    }
 
    // Back substitution: L^T * x = y
    std::vector<double> x(n);
    for (std::size_t i = n; i > 0; --i) {
        std::size_t idx = i - 1;
        double sum = 0.0;
        for (std::size_t j = idx + 1; j < n; ++j) {
            sum += L[j][idx] * x[j];
        }
        x[idx] = (y[idx] - sum) / L[idx][idx];
    }
    return x;
}
 
inline std::vector<std::vector<double>> inverse_cholesky(
    const std::vector<std::vector<double>>& L)
{
    std::size_t n = L.size();
    std::vector<std::vector<double>> inv(n, std::vector<double>(n, 0.0));
 
    // Solve each column
    for (std::size_t j = 0; j < n; ++j) {
        std::vector<double> e(n, 0.0);
        e[j] = 1.0;
        std::vector<double> col = solve_cholesky(L, e);
        for (std::size_t i = 0; i < n; ++i) {
            inv[i][j] = col[i];
        }
    }
    return inv;
}
 
} // namespace detail
 
inline multiple_regression_result multiple_linear_regression(
    const std::vector<std::vector<double>>& X,
    const std::vector<double>& y)
{
    // Verify that X does not contain an intercept column
    detail::validate_no_intercept_column(X, "multiple_linear_regression");
 
    std::size_t n = X.size();
    if (n == 0) {
        throw std::invalid_argument("statcpp::multiple_linear_regression: empty data");
    }
    if (n != y.size()) {
        throw std::invalid_argument("statcpp::multiple_linear_regression: X and y must have same number of observations");
    }
 
    std::size_t p = X[0].size();  // Number of predictors (excluding intercept)
    for (const auto& row : X) {
        if (row.size() != p) {
            throw std::invalid_argument("statcpp::multiple_linear_regression: inconsistent number of predictors");
        }
    }
 
    std::size_t p_full = p + 1;  // Number of coefficients including intercept
    if (n <= p_full) {
        throw std::invalid_argument("statcpp::multiple_linear_regression: need more observations than predictors");
    }
 
    double n_d = static_cast<double>(n);
 
    // Create design matrix with intercept term
    std::vector<std::vector<double>> X_design(n, std::vector<double>(p_full));
    for (std::size_t i = 0; i < n; ++i) {
        X_design[i][0] = 1.0;  // Intercept
        for (std::size_t j = 0; j < p; ++j) {
            X_design[i][j + 1] = X[i][j];
        }
    }
 
    // Calculate X^T * X
    auto Xt = detail::transpose(X_design);
    auto XtX = detail::matrix_multiply(Xt, X_design);
 
    // Calculate X^T * y
    auto Xty = detail::matrix_vector_multiply(Xt, y);
 
    // Solve for coefficients using Cholesky decomposition
    auto L = detail::cholesky(XtX);
    auto coefficients = detail::solve_cholesky(L, Xty);
 
    // Calculate (X^T X)^{-1}
    auto XtX_inv = detail::inverse_cholesky(L);
 
    // Calculate predicted values and residuals
    std::vector<double> y_hat(n);
    std::vector<double> residuals(n);
    double mean_y = statcpp::mean(y.begin(), y.end());
 
    double ss_total = 0.0;
    double ss_residual = 0.0;
 
    for (std::size_t i = 0; i < n; ++i) {
        double pred = 0.0;
        for (std::size_t j = 0; j < p_full; ++j) {
            pred += X_design[i][j] * coefficients[j];
        }
        y_hat[i] = pred;
        residuals[i] = y[i] - pred;
 
        ss_total += (y[i] - mean_y) * (y[i] - mean_y);
        ss_residual += residuals[i] * residuals[i];
    }
 
    double ss_regression = ss_total - ss_residual;
 
    if (ss_total == 0.0) {
        throw std::invalid_argument("statcpp::multiple_linear_regression: zero variance in y (constant response)");
    }
 
    // Degrees of freedom
    double df_reg = static_cast<double>(p);
    double df_res = n_d - static_cast<double>(p_full);
 
    // Residual variance
    double mse = ss_residual / df_res;
    double residual_se = std::sqrt(mse);
 
    // Standard errors, t-statistics, and p-values of coefficients
    std::vector<double> coefficient_se(p_full);
    std::vector<double> t_statistics(p_full);
    std::vector<double> p_values(p_full);
 
    for (std::size_t j = 0; j < p_full; ++j) {
        coefficient_se[j] = std::sqrt(mse * XtX_inv[j][j]);
        t_statistics[j] = coefficients[j] / coefficient_se[j];
        p_values[j] = 2.0 * (1.0 - t_cdf(std::abs(t_statistics[j]), df_res));
    }
 
    // Coefficient of determination
    double r_squared = ss_regression / ss_total;
    double adj_r_squared = 1.0 - (1.0 - r_squared) * (n_d - 1.0) / df_res;
 
    // F-statistic
    double f_statistic = (ss_regression / df_reg) / mse;
    double f_p_value = 1.0 - f_cdf(f_statistic, df_reg, df_res);
 
    return {
        coefficients, coefficient_se, t_statistics, p_values,
        r_squared, adj_r_squared,
        residual_se,
        f_statistic, f_p_value,
        df_reg, df_res,
        ss_total, ss_regression, ss_residual
    };
}
 
// ============================================================================
// Prediction
// ============================================================================
 
inline double predict(const simple_regression_result& model, double x)
{
    return model.intercept + model.slope * x;
}
 
inline double predict(const multiple_regression_result& model, const std::vector<double>& x)
{
    if (x.size() + 1 != model.coefficients.size()) {
        throw std::invalid_argument("statcpp::predict: x dimension mismatch");
    }
 
    double pred = model.coefficients[0];  // Intercept
    for (std::size_t i = 0; i < x.size(); ++i) {
        pred += model.coefficients[i + 1] * x[i];
    }
    return pred;
}
 
// ============================================================================
// Prediction Interval
// ============================================================================
 
template <typename IteratorX>
prediction_interval prediction_interval_simple(
    const simple_regression_result& model,
    IteratorX x_first, IteratorX x_last,
    double x_new,
    double confidence = 0.95)
{
    if (confidence <= 0.0 || confidence >= 1.0) {
        throw std::invalid_argument("statcpp::prediction_interval_simple: confidence must be in (0, 1)");
    }
 
    auto n = statcpp::count(x_first, x_last);
    double n_d = static_cast<double>(n);
    double mean_x = statcpp::mean(x_first, x_last);
 
    // Calculate Sxx
    double sxx = 0.0;
    for (auto it = x_first; it != x_last; ++it) {
        double dx = static_cast<double>(*it) - mean_x;
        sxx += dx * dx;
    }
 
    double y_hat = predict(model, x_new);
 
    // Standard error of prediction (for new observation)
    double dx_new = x_new - mean_x;
    double se_pred = model.residual_se * std::sqrt(1.0 + 1.0 / n_d + dx_new * dx_new / sxx);
 
    double t_crit = t_quantile(1.0 - (1.0 - confidence) / 2.0, model.df_residual);
    double margin = t_crit * se_pred;
 
    return {y_hat, y_hat - margin, y_hat + margin, se_pred};
}
 
template <typename IteratorX>
prediction_interval confidence_interval_mean(
    const simple_regression_result& model,
    IteratorX x_first, IteratorX x_last,
    double x_new,
    double confidence = 0.95)
{
    if (confidence <= 0.0 || confidence >= 1.0) {
        throw std::invalid_argument("statcpp::confidence_interval_mean: confidence must be in (0, 1)");
    }
 
    auto n = statcpp::count(x_first, x_last);
    double n_d = static_cast<double>(n);
    double mean_x = statcpp::mean(x_first, x_last);
 
    // Calculate Sxx
    double sxx = 0.0;
    for (auto it = x_first; it != x_last; ++it) {
        double dx = static_cast<double>(*it) - mean_x;
        sxx += dx * dx;
    }
 
    double y_hat = predict(model, x_new);
 
    // Standard error of mean prediction
    double dx_new = x_new - mean_x;
    double se_mean = model.residual_se * std::sqrt(1.0 / n_d + dx_new * dx_new / sxx);
 
    double t_crit = t_quantile(1.0 - (1.0 - confidence) / 2.0, model.df_residual);
    double margin = t_crit * se_mean;
 
    return {y_hat, y_hat - margin, y_hat + margin, se_mean};
}
 
// ============================================================================
// Residual Diagnostics
// ============================================================================
 
template <typename IteratorX, typename IteratorY>
residual_diagnostics compute_residual_diagnostics(
    const simple_regression_result& model,
    IteratorX x_first, IteratorX x_last,
    IteratorY y_first, IteratorY y_last)
{
    auto n = statcpp::count(x_first, x_last);
    if (n != statcpp::count(y_first, y_last)) {
        throw std::invalid_argument("statcpp::compute_residual_diagnostics: x and y must have same length");
    }
 
    double n_d = static_cast<double>(n);
    double mean_x = statcpp::mean(x_first, x_last);
 
    // Calculate Sxx
    double sxx = 0.0;
    for (auto it = x_first; it != x_last; ++it) {
        double dx = static_cast<double>(*it) - mean_x;
        sxx += dx * dx;
    }
 
    std::vector<double> residuals(n);
    std::vector<double> hat_values(n);
    std::vector<double> standardized_residuals(n);
 
    auto it_x = x_first;
    auto it_y = y_first;
    for (std::size_t i = 0; it_x != x_last; ++it_x, ++it_y, ++i) {
        double x_i = static_cast<double>(*it_x);
        double y_i = static_cast<double>(*it_y);
        double y_hat = predict(model, x_i);
        residuals[i] = y_i - y_hat;
 
        // Leverage h_ii = 1/n + (x_i - mean_x)^2 / Sxx
        double dx = x_i - mean_x;
        hat_values[i] = 1.0 / n_d + dx * dx / sxx;
 
        // Standardized residuals
        standardized_residuals[i] = residuals[i] / model.residual_se;
    }
 
    // Studentized residuals and Cook's distance
    std::vector<double> studentized_residuals(n);
    std::vector<double> cooks_distance(n);
    double p = 2.0;  // Number of coefficients (intercept + slope)
 
    for (std::size_t i = 0; i < n; ++i) {
        double h_i = hat_values[i];
        double se_i = model.residual_se * std::sqrt(1.0 - h_i);
        studentized_residuals[i] = (se_i > 0.0) ? residuals[i] / se_i : 0.0;
        cooks_distance[i] = (standardized_residuals[i] * standardized_residuals[i] / p)
                          * (h_i / ((1.0 - h_i) * (1.0 - h_i)));
    }
 
    // Durbin-Watson statistic
    double dw_num = 0.0;
    double dw_den = 0.0;
    for (std::size_t i = 0; i < n; ++i) {
        dw_den += residuals[i] * residuals[i];
        if (i > 0) {
            double diff = residuals[i] - residuals[i - 1];
            dw_num += diff * diff;
        }
    }
    double durbin_watson = (dw_den > 0.0) ? dw_num / dw_den : 0.0;
 
    return {residuals, standardized_residuals, studentized_residuals,
            hat_values, cooks_distance, durbin_watson};
}
 
inline residual_diagnostics compute_residual_diagnostics(
    const multiple_regression_result& model,
    const std::vector<std::vector<double>>& X,
    const std::vector<double>& y)
{
    std::size_t n = X.size();
    if (n != y.size()) {
        throw std::invalid_argument("statcpp::compute_residual_diagnostics: X and y must have same length");
    }
 
    std::size_t p = X[0].size();
    std::size_t p_full = p + 1;
 
    // Create design matrix
    std::vector<std::vector<double>> X_design(n, std::vector<double>(p_full));
    for (std::size_t i = 0; i < n; ++i) {
        X_design[i][0] = 1.0;
        for (std::size_t j = 0; j < p; ++j) {
            X_design[i][j + 1] = X[i][j];
        }
    }
 
    // Calculate (X^T X)^{-1}
    auto Xt = detail::transpose(X_design);
    auto XtX = detail::matrix_multiply(Xt, X_design);
    auto L = detail::cholesky(XtX);
    auto XtX_inv = detail::inverse_cholesky(L);
 
    // H = X(X^T X)^{-1}X^T leverage values (diagonal elements only)
    std::vector<double> hat_values(n);
    for (std::size_t i = 0; i < n; ++i) {
        double h_ii = 0.0;
        for (std::size_t j = 0; j < p_full; ++j) {
            for (std::size_t k = 0; k < p_full; ++k) {
                h_ii += X_design[i][j] * XtX_inv[j][k] * X_design[i][k];
            }
        }
        hat_values[i] = h_ii;
    }
 
    // Calculate residuals
    std::vector<double> residuals(n);
    for (std::size_t i = 0; i < n; ++i) {
        double pred = predict(model, X[i]);
        residuals[i] = y[i] - pred;
    }
 
    // Standardized residuals, studentized residuals, Cook's distance
    std::vector<double> standardized_residuals(n);
    std::vector<double> studentized_residuals(n);
    std::vector<double> cooks_distance(n);
 
    double p_d = static_cast<double>(p_full);
 
    for (std::size_t i = 0; i < n; ++i) {
        standardized_residuals[i] = residuals[i] / model.residual_se;
 
        double h_i = hat_values[i];
        double se_i = model.residual_se * std::sqrt(1.0 - h_i);
        studentized_residuals[i] = (se_i > 0.0) ? residuals[i] / se_i : 0.0;
 
        cooks_distance[i] = (standardized_residuals[i] * standardized_residuals[i] / p_d)
                          * (h_i / ((1.0 - h_i) * (1.0 - h_i)));
    }
 
    // Durbin-Watson statistic
    double dw_num = 0.0;
    double dw_den = 0.0;
    for (std::size_t i = 0; i < n; ++i) {
        dw_den += residuals[i] * residuals[i];
        if (i > 0) {
            double diff = residuals[i] - residuals[i - 1];
            dw_num += diff * diff;
        }
    }
    double durbin_watson = (dw_den > 0.0) ? dw_num / dw_den : 0.0;
 
    return {residuals, standardized_residuals, studentized_residuals,
            hat_values, cooks_distance, durbin_watson};
}
 
// ============================================================================
// VIF (Variance Inflation Factor)
// ============================================================================
 
inline std::vector<double> compute_vif(const std::vector<std::vector<double>>& X)
{
    std::size_t n = X.size();
    if (n < 3) {
        throw std::invalid_argument("statcpp::compute_vif: need at least 3 observations");
    }
 
    std::size_t p = X[0].size();
    if (p < 2) {
        throw std::invalid_argument("statcpp::compute_vif: need at least 2 predictors");
    }
 
    std::vector<double> vif(p);
 
    for (std::size_t j = 0; j < p; ++j) {
        // Regress j-th variable as response on other variables as predictors
        std::vector<double> y_j(n);
        std::vector<std::vector<double>> X_others(n, std::vector<double>(p - 1));
 
        for (std::size_t i = 0; i < n; ++i) {
            y_j[i] = X[i][j];
            std::size_t col = 0;
            for (std::size_t k = 0; k < p; ++k) {
                if (k != j) {
                    X_others[i][col++] = X[i][k];
                }
            }
        }
 
        auto result = multiple_linear_regression(X_others, y_j);
        double r_sq = result.r_squared;
 
        // VIF = 1 / (1 - R^2)
        if (r_sq >= 1.0) {
            vif[j] = std::numeric_limits<double>::infinity();
        } else {
            vif[j] = 1.0 / (1.0 - r_sq);
        }
    }
 
    return vif;
}
 
// ============================================================================
// Multicollinearity Diagnostics (Extended)
// ============================================================================
 
inline double correlation_matrix_determinant(const std::vector<std::vector<double>>& X)
{
    std::size_t n = X.size();
    if (n < 2) {
        throw std::invalid_argument("statcpp::correlation_matrix_determinant: need at least 2 observations");
    }
 
    std::size_t p = X[0].size();
    if (p < 2) {
        throw std::invalid_argument("statcpp::correlation_matrix_determinant: need at least 2 predictors");
    }
 
    // Calculate correlation matrix
    std::vector<std::vector<double>> corr_matrix(p, std::vector<double>(p));
    for (std::size_t i = 0; i < p; ++i) {
        for (std::size_t j = 0; j < p; ++j) {
            if (i == j) {
                corr_matrix[i][j] = 1.0;
            } else if (j > i) {
                // Extract each column
                std::vector<double> col_i(n), col_j(n);
                for (std::size_t k = 0; k < n; ++k) {
                    col_i[k] = X[k][i];
                    col_j[k] = X[k][j];
                }
                double corr = pearson_correlation(col_i.begin(), col_i.end(),
                                                   col_j.begin(), col_j.end());
                corr_matrix[i][j] = corr;
                corr_matrix[j][i] = corr;
            }
        }
    }
 
    // Calculate determinant (only supports small matrices)
    if (p == 2) {
        // 2x2 matrix: det = a11*a22 - a12*a21
        return corr_matrix[0][0] * corr_matrix[1][1] - corr_matrix[0][1] * corr_matrix[1][0];
    } else if (p == 3) {
        // 3x3 matrix: Sarrus' rule
        double a = corr_matrix[0][0] * corr_matrix[1][1] * corr_matrix[2][2];
        double b = corr_matrix[0][1] * corr_matrix[1][2] * corr_matrix[2][0];
        double c = corr_matrix[0][2] * corr_matrix[1][0] * corr_matrix[2][1];
        double d = corr_matrix[0][2] * corr_matrix[1][1] * corr_matrix[2][0];
        double e = corr_matrix[0][0] * corr_matrix[1][2] * corr_matrix[2][1];
        double f = corr_matrix[0][1] * corr_matrix[1][0] * corr_matrix[2][2];
        return a + b + c - d - e - f;
    } else {
        // Larger matrices would require LU decomposition, not supported here
        throw std::invalid_argument("statcpp::correlation_matrix_determinant: only 2 or 3 predictors supported");
    }
}
 
inline double multicollinearity_score(const std::vector<std::vector<double>>& X)
{
    double det = correlation_matrix_determinant(X);
    // Determinant can be negative, but what matters is how close to 0
    return 1.0 - std::abs(det);
}
 
// ============================================================================
// R-squared and Related Measures
// ============================================================================
 
template <typename IteratorY, typename IteratorPred>
double r_squared(IteratorY y_first, IteratorY y_last,
                  IteratorPred pred_first, IteratorPred pred_last)
{
    auto n_y = statcpp::count(y_first, y_last);
    auto n_pred = statcpp::count(pred_first, pred_last);
 
    if (n_y != n_pred) {
        throw std::invalid_argument("statcpp::r_squared: y and predictions must have same length");
    }
    if (n_y < 2) {
        throw std::invalid_argument("statcpp::r_squared: need at least 2 observations");
    }
 
    double mean_y = statcpp::mean(y_first, y_last);
 
    double ss_total = 0.0;
    double ss_residual = 0.0;
 
    auto it_y = y_first;
    auto it_pred = pred_first;
    for (; it_y != y_last; ++it_y, ++it_pred) {
        double y_i = static_cast<double>(*it_y);
        double pred_i = static_cast<double>(*it_pred);
        ss_total += (y_i - mean_y) * (y_i - mean_y);
        ss_residual += (y_i - pred_i) * (y_i - pred_i);
    }
 
    if (ss_total == 0.0) {
        return 1.0;  // Completely constant case
    }
 
    return 1.0 - ss_residual / ss_total;
}
 
template <typename IteratorY, typename IteratorPred>
double adjusted_r_squared(IteratorY y_first, IteratorY y_last,
                           IteratorPred pred_first, IteratorPred pred_last,
                           std::size_t num_predictors)
{
    auto n = statcpp::count(y_first, y_last);
    double n_d = static_cast<double>(n);
    double p = static_cast<double>(num_predictors);
 
    if (n_d <= p + 1.0) {
        throw std::invalid_argument("statcpp::adjusted_r_squared: need more observations than predictors");
    }
 
    double r_sq = r_squared(y_first, y_last, pred_first, pred_last);
    return 1.0 - (1.0 - r_sq) * (n_d - 1.0) / (n_d - p - 1.0);
}
 
} // namespace statcpp