categorical.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <stdexcept>
#include <vector>
 
namespace statcpp {
 
// ============================================================================
// Contingency Table (Cross-tabulation)
// ============================================================================
 
struct contingency_table_result {
    std::vector<std::vector<std::size_t>> table;  
    std::vector<std::size_t> row_totals;          
    std::vector<std::size_t> col_totals;          
    std::size_t total;                            
    std::size_t n_rows;                           
    std::size_t n_cols;                           
};
 
inline contingency_table_result contingency_table(
    const std::vector<std::size_t>& row_data,
    const std::vector<std::size_t>& col_data)
{
    if (row_data.size() != col_data.size()) {
        throw std::invalid_argument("statcpp::contingency_table: data lengths must match");
    }
    if (row_data.empty()) {
        throw std::invalid_argument("statcpp::contingency_table: empty data");
    }
 
    // Determine number of categories
    std::size_t max_row = *std::max_element(row_data.begin(), row_data.end());
    std::size_t max_col = *std::max_element(col_data.begin(), col_data.end());
    std::size_t n_rows = max_row + 1;
    std::size_t n_cols = max_col + 1;
 
    // Create contingency table
    std::vector<std::vector<std::size_t>> table(n_rows, std::vector<std::size_t>(n_cols, 0));
 
    for (std::size_t i = 0; i < row_data.size(); ++i) {
        table[row_data[i]][col_data[i]]++;
    }
 
    // Calculate marginal frequencies
    std::vector<std::size_t> row_totals(n_rows, 0);
    std::vector<std::size_t> col_totals(n_cols, 0);
    std::size_t total = 0;
 
    for (std::size_t i = 0; i < n_rows; ++i) {
        for (std::size_t j = 0; j < n_cols; ++j) {
            row_totals[i] += table[i][j];
            col_totals[j] += table[i][j];
            total += table[i][j];
        }
    }
 
    return {table, row_totals, col_totals, total, n_rows, n_cols};
}
 
// ============================================================================
// Odds Ratio
// ============================================================================
 
struct odds_ratio_result {
    double odds_ratio;            
    double log_odds_ratio;        
    double se_log_odds_ratio;     
    double ci_lower;              
    double ci_upper;              
};
 
inline odds_ratio_result odds_ratio(const std::vector<std::vector<std::size_t>>& table)
{
    if (table.size() != 2 || table[0].size() != 2 || table[1].size() != 2) {
        throw std::invalid_argument("statcpp::odds_ratio: table must be 2x2");
    }
 
    double a = static_cast<double>(table[0][0]);
    double b = static_cast<double>(table[0][1]);
    double c = static_cast<double>(table[1][0]);
    double d = static_cast<double>(table[1][1]);
 
    if (a == 0 || b == 0 || c == 0 || d == 0) {
        throw std::invalid_argument("statcpp::odds_ratio: zero cell count not allowed");
    }
 
    double or_val = (a * d) / (b * c);
    double log_or = std::log(or_val);
    double se_log_or = std::sqrt(1.0/a + 1.0/b + 1.0/c + 1.0/d);
 
    // 95% confidence interval
    double z = 1.96;
    double ci_lower = std::exp(log_or - z * se_log_or);
    double ci_upper = std::exp(log_or + z * se_log_or);
 
    return {or_val, log_or, se_log_or, ci_lower, ci_upper};
}
 
inline odds_ratio_result odds_ratio(std::size_t a, std::size_t b, std::size_t c, std::size_t d)
{
    return odds_ratio({{a, b}, {c, d}});
}
 
// ============================================================================
// Relative Risk (Risk Ratio)
// ============================================================================
 
struct relative_risk_result {
    double relative_risk;         
    double log_relative_risk;     
    double se_log_relative_risk;  
    double ci_lower;              
    double ci_upper;              
};
 
inline relative_risk_result relative_risk(const std::vector<std::vector<std::size_t>>& table)
{
    if (table.size() != 2 || table[0].size() != 2 || table[1].size() != 2) {
        throw std::invalid_argument("statcpp::relative_risk: table must be 2x2");
    }
 
    double a = static_cast<double>(table[0][0]);
    double b = static_cast<double>(table[0][1]);
    double c = static_cast<double>(table[1][0]);
    double d = static_cast<double>(table[1][1]);
 
    double n1 = a + b;  // Size of exposed group
    double n0 = c + d;  // Size of unexposed group
 
    if (n1 == 0 || n0 == 0) {
        throw std::invalid_argument("statcpp::relative_risk: zero row total");
    }
    if (a == 0 || c == 0) {
        throw std::invalid_argument("statcpp::relative_risk: zero risk in a group");
    }
 
    double risk1 = a / n1;  // Risk in exposed group
    double risk0 = c / n0;  // Risk in unexposed group
 
    double rr = risk1 / risk0;
    double log_rr = std::log(rr);
 
    // Standard error (Greenland-Robins formula)
    double se_log_rr = std::sqrt((1.0 - risk1)/(a) + (1.0 - risk0)/(c));
 
    // 95% confidence interval
    double z = 1.96;
    double ci_lower = std::exp(log_rr - z * se_log_rr);
    double ci_upper = std::exp(log_rr + z * se_log_rr);
 
    return {rr, log_rr, se_log_rr, ci_lower, ci_upper};
}
 
inline relative_risk_result relative_risk(std::size_t a, std::size_t b, std::size_t c, std::size_t d)
{
    return relative_risk({{a, b}, {c, d}});
}
 
// ============================================================================
// Risk Difference (Attributable Risk)
// ============================================================================
 
struct risk_difference_result {
    double risk_difference;       
    double se;                    
    double ci_lower;              
    double ci_upper;              
};
 
inline risk_difference_result risk_difference(const std::vector<std::vector<std::size_t>>& table)
{
    if (table.size() != 2 || table[0].size() != 2 || table[1].size() != 2) {
        throw std::invalid_argument("statcpp::risk_difference: table must be 2x2");
    }
 
    double a = static_cast<double>(table[0][0]);
    double b = static_cast<double>(table[0][1]);
    double c = static_cast<double>(table[1][0]);
    double d = static_cast<double>(table[1][1]);
 
    double n1 = a + b;
    double n0 = c + d;
 
    if (n1 == 0 || n0 == 0) {
        throw std::invalid_argument("statcpp::risk_difference: zero row total");
    }
 
    double risk1 = a / n1;
    double risk0 = c / n0;
 
    double rd = risk1 - risk0;
 
    // Standard error
    double se = std::sqrt(risk1 * (1.0 - risk1) / n1 + risk0 * (1.0 - risk0) / n0);
 
    // 95% confidence interval
    double z = 1.96;
    double ci_lower = rd - z * se;
    double ci_upper = rd + z * se;
 
    return {rd, se, ci_lower, ci_upper};
}
 
inline risk_difference_result risk_difference(std::size_t a, std::size_t b, std::size_t c, std::size_t d)
{
    return risk_difference({{a, b}, {c, d}});
}
 
// ============================================================================
// Number Needed to Treat (NNT)
// ============================================================================
 
inline double number_needed_to_treat(const std::vector<std::vector<std::size_t>>& table)
{
    auto rd = risk_difference(table);
    if (rd.risk_difference == 0.0) {
        throw std::invalid_argument("statcpp::number_needed_to_treat: risk difference is zero");
    }
    return 1.0 / std::abs(rd.risk_difference);
}
 
} // namespace statcpp