anova.hpp

Go to the documentation of this file.

 
#pragma once
 
#include "statcpp/continuous_distributions.hpp"
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <limits>
#include <numeric>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
 
namespace statcpp {
 
// ============================================================================
// ANOVA Result Structures
// ============================================================================
 
struct anova_row {
    std::string source;     
    double ss;              
    double df;              
    double ms;              
    double f_statistic;     
    double p_value;         
};
 
struct one_way_anova_result {
    anova_row between;      
    anova_row within;       
    double ss_total;        
    double df_total;        
    std::size_t n_groups;   
    std::size_t n_total;    
    double grand_mean;      
    std::vector<double> group_means;       
    std::vector<std::size_t> group_sizes;  
};
 
struct two_way_anova_result {
    anova_row factor_a;         
    anova_row factor_b;         
    anova_row interaction;      
    anova_row error;            
    double ss_total;            
    double df_total;            
    std::size_t levels_a;       
    std::size_t levels_b;       
    std::size_t n_total;        
    double grand_mean;          
};
 
struct posthoc_comparison {
    std::size_t group1;     
    std::size_t group2;     
    double mean_diff;       
    double se;              
    double statistic;       
    double p_value;         
    double lower;           
    double upper;           
    bool significant;       
};
 
struct posthoc_result {
    std::string method;                          
    std::vector<posthoc_comparison> comparisons; 
    double alpha;                                
    double mse;                                  
    double df_error;                             
};
 
// ============================================================================
// One-Way ANOVA
// ============================================================================
 
inline one_way_anova_result one_way_anova(const std::vector<std::vector<double>>& groups)
{
    std::size_t k = groups.size();  // Number of groups
    if (k < 2) {
        throw std::invalid_argument("statcpp::one_way_anova: need at least 2 groups");
    }
 
    // Calculate size and mean for each group
    std::vector<std::size_t> group_sizes(k);
    std::vector<double> group_means(k);
    std::size_t n_total = 0;
    double grand_sum = 0.0;
 
    for (std::size_t i = 0; i < k; ++i) {
        if (groups[i].empty()) {
            throw std::invalid_argument("statcpp::one_way_anova: empty group detected");
        }
        group_sizes[i] = groups[i].size();
        n_total += group_sizes[i];
        double group_sum = std::accumulate(groups[i].begin(), groups[i].end(), 0.0);
        group_means[i] = group_sum / static_cast<double>(group_sizes[i]);
        grand_sum += group_sum;
    }
 
    if (n_total <= k) {
        throw std::invalid_argument("statcpp::one_way_anova: need more observations than groups");
    }
 
    double grand_mean = grand_sum / static_cast<double>(n_total);
 
    // Calculate sum of squares
    double ss_between = 0.0;  // Between-group sum of squares
    double ss_within = 0.0;   // Within-group sum of squares
 
    for (std::size_t i = 0; i < k; ++i) {
        double diff = group_means[i] - grand_mean;
        ss_between += static_cast<double>(group_sizes[i]) * diff * diff;
 
        for (double x : groups[i]) {
            double diff_within = x - group_means[i];
            ss_within += diff_within * diff_within;
        }
    }
 
    double ss_total = ss_between + ss_within;
 
    // Degrees of freedom
    double df_between = static_cast<double>(k - 1);
    double df_within = static_cast<double>(n_total - k);
    double df_total = static_cast<double>(n_total - 1);
 
    // Mean squares
    double ms_between = ss_between / df_between;
    double ms_within = ss_within / df_within;
 
    // F statistic and p-value (handle degenerate case ms_within == 0)
    double f_statistic;
    double p_value;
    if (ms_within == 0.0) {
        f_statistic = (ms_between == 0.0) ? 0.0 : std::numeric_limits<double>::infinity();
        p_value = (ms_between == 0.0) ? 1.0 : 0.0;
    } else {
        f_statistic = ms_between / ms_within;
        p_value = 1.0 - f_cdf(f_statistic, df_between, df_within);
    }
 
    anova_row between{"Between Groups", ss_between, df_between, ms_between, f_statistic, p_value};
    anova_row within{"Within Groups", ss_within, df_within, ms_within, 0.0, 0.0};
 
    return {between, within, ss_total, df_total, k, n_total, grand_mean, group_means, group_sizes};
}
 
// ============================================================================
// Two-Way ANOVA
// ============================================================================
 
inline two_way_anova_result two_way_anova(
    const std::vector<std::vector<std::vector<double>>>& data)
{
    std::size_t a = data.size();  // Number of levels for factor A
    if (a < 2) {
        throw std::invalid_argument("statcpp::two_way_anova: need at least 2 levels for factor A");
    }
 
    std::size_t b = data[0].size();  // Number of levels for factor B
    if (b < 2) {
        throw std::invalid_argument("statcpp::two_way_anova: need at least 2 levels for factor B");
    }
 
    // Check if all levels have the same number of replications
    std::size_t n_rep = data[0][0].size();  // Number of replications
    if (n_rep < 2) {
        throw std::invalid_argument("statcpp::two_way_anova: at least 2 replications per cell are required (n_rep >= 2)");
    }
    std::size_t n_total = 0;
    double grand_sum = 0.0;
 
    for (std::size_t i = 0; i < a; ++i) {
        if (data[i].size() != b) {
            throw std::invalid_argument("statcpp::two_way_anova: inconsistent number of levels for factor B");
        }
        for (std::size_t j = 0; j < b; ++j) {
            if (data[i][j].size() != n_rep) {
                throw std::invalid_argument("statcpp::two_way_anova: unequal cell sizes not supported");
            }
            n_total += n_rep;
            for (double x : data[i][j]) {
                grand_sum += x;
            }
        }
    }
 
    double grand_mean = grand_sum / static_cast<double>(n_total);
    double n_rep_d = static_cast<double>(n_rep);
    double a_d = static_cast<double>(a);
    double b_d = static_cast<double>(b);
 
    // Calculate means for each level
    std::vector<double> mean_a(a, 0.0);  // Mean for each level of factor A
    std::vector<double> mean_b(b, 0.0);  // Mean for each level of factor B
    std::vector<std::vector<double>> mean_ab(a, std::vector<double>(b, 0.0));  // Cell means
 
    for (std::size_t i = 0; i < a; ++i) {
        for (std::size_t j = 0; j < b; ++j) {
            double cell_sum = std::accumulate(data[i][j].begin(), data[i][j].end(), 0.0);
            mean_ab[i][j] = cell_sum / n_rep_d;
            mean_a[i] += cell_sum;
            mean_b[j] += cell_sum;
        }
        mean_a[i] /= (b_d * n_rep_d);
    }
    for (std::size_t j = 0; j < b; ++j) {
        mean_b[j] /= (a_d * n_rep_d);
    }
 
    // Calculate sum of squares
    double ss_a = 0.0;      // Sum of squares for factor A
    double ss_b = 0.0;      // Sum of squares for factor B
    double ss_ab = 0.0;     // Sum of squares for interaction
    double ss_error = 0.0;  // Error sum of squares
 
    for (std::size_t i = 0; i < a; ++i) {
        double diff_a = mean_a[i] - grand_mean;
        ss_a += diff_a * diff_a;
    }
    ss_a *= b_d * n_rep_d;
 
    for (std::size_t j = 0; j < b; ++j) {
        double diff_b = mean_b[j] - grand_mean;
        ss_b += diff_b * diff_b;
    }
    ss_b *= a_d * n_rep_d;
 
    for (std::size_t i = 0; i < a; ++i) {
        for (std::size_t j = 0; j < b; ++j) {
            double interaction = mean_ab[i][j] - mean_a[i] - mean_b[j] + grand_mean;
            ss_ab += interaction * interaction;
 
            for (double x : data[i][j]) {
                double error = x - mean_ab[i][j];
                ss_error += error * error;
            }
        }
    }
    ss_ab *= n_rep_d;
 
    double ss_total = ss_a + ss_b + ss_ab + ss_error;
 
    // Degrees of freedom
    double df_a = a_d - 1.0;
    double df_b = b_d - 1.0;
    double df_ab = df_a * df_b;
    double df_error = static_cast<double>(n_total) - a_d * b_d;
    double df_total = static_cast<double>(n_total) - 1.0;
 
    // Mean squares
    double ms_a = ss_a / df_a;
    double ms_b = ss_b / df_b;
    double ms_ab = ss_ab / df_ab;
    double ms_error = ss_error / df_error;
 
    // F statistics and p-values (handle degenerate case ms_error == 0)
    double f_a, f_b, f_ab;
    double p_a, p_b, p_ab;
    if (ms_error == 0.0) {
        f_a = (ms_a == 0.0) ? 0.0 : std::numeric_limits<double>::infinity();
        f_b = (ms_b == 0.0) ? 0.0 : std::numeric_limits<double>::infinity();
        f_ab = (ms_ab == 0.0) ? 0.0 : std::numeric_limits<double>::infinity();
        p_a = (ms_a == 0.0) ? 1.0 : 0.0;
        p_b = (ms_b == 0.0) ? 1.0 : 0.0;
        p_ab = (ms_ab == 0.0) ? 1.0 : 0.0;
    } else {
        f_a = ms_a / ms_error;
        f_b = ms_b / ms_error;
        f_ab = ms_ab / ms_error;
        p_a = 1.0 - f_cdf(f_a, df_a, df_error);
        p_b = 1.0 - f_cdf(f_b, df_b, df_error);
        p_ab = 1.0 - f_cdf(f_ab, df_ab, df_error);
    }
 
    anova_row factor_a{"Factor A", ss_a, df_a, ms_a, f_a, p_a};
    anova_row factor_b{"Factor B", ss_b, df_b, ms_b, f_b, p_b};
    anova_row interaction{"A x B", ss_ab, df_ab, ms_ab, f_ab, p_ab};
    anova_row error{"Error", ss_error, df_error, ms_error, 0.0, 0.0};
 
    return {factor_a, factor_b, interaction, error,
            ss_total, df_total, a, b, n_total, grand_mean};
}
 
// ============================================================================
// Post-hoc Comparisons
// ============================================================================
 
inline posthoc_result tukey_hsd(const one_way_anova_result& anova_result,
                                 const std::vector<std::vector<double>>& groups,
                                 double alpha = 0.05)
{
    (void)groups;  // Currently unused; all statistics derived from anova_result
 
    if (alpha <= 0.0 || alpha >= 1.0) {
        throw std::invalid_argument("statcpp::tukey_hsd: alpha must be in (0, 1)");
    }
 
    std::size_t k = anova_result.n_groups;
    double mse = anova_result.within.ms;
    double df_error = anova_result.within.df;
    double k_d = static_cast<double>(k);
 
    double q_crit = studentized_range_quantile(1.0 - alpha, k_d, df_error);
 
    std::vector<posthoc_comparison> comparisons;
 
    // Compare all pairs
    for (std::size_t i = 0; i < k; ++i) {
        for (std::size_t j = i + 1; j < k; ++j) {
            double mean_diff = anova_result.group_means[i] - anova_result.group_means[j];
            double n_i = static_cast<double>(anova_result.group_sizes[i]);
            double n_j = static_cast<double>(anova_result.group_sizes[j]);
 
            // Tukey-Kramer standard error
            double se = std::sqrt(mse * 0.5 * (1.0 / n_i + 1.0 / n_j));
 
            double q_stat, p_value, lower, upper;
            bool significant;
 
            if (se == 0.0) {
                // Degenerate case: within-group variance is zero
                if (mean_diff == 0.0) {
                    q_stat = 0.0; p_value = 1.0;
                    lower = 0.0; upper = 0.0; significant = false;
                } else {
                    q_stat = std::numeric_limits<double>::infinity();
                    p_value = 0.0;
                    lower = mean_diff; upper = mean_diff; significant = true;
                }
            } else {
                // q statistic (studentized range statistic)
                q_stat = std::abs(mean_diff) / se;
 
                // p-value from studentized range distribution
                p_value = 1.0 - studentized_range_cdf(q_stat, k_d, df_error);
                p_value = std::max(0.0, std::min(1.0, p_value));
 
                double margin = q_crit * se;
                lower = mean_diff - margin;
                upper = mean_diff + margin;
                significant = (q_stat > q_crit);
            }
 
            comparisons.push_back({i, j, mean_diff, se, q_stat, p_value, lower, upper, significant});
        }
    }
 
    return {"Tukey HSD", comparisons, alpha, mse, df_error};
}
 
inline posthoc_result bonferroni_posthoc(const one_way_anova_result& anova_result,
                                          double alpha = 0.05)
{
    if (alpha <= 0.0 || alpha >= 1.0) {
        throw std::invalid_argument("statcpp::bonferroni_posthoc: alpha must be in (0, 1)");
    }
 
    std::size_t k = anova_result.n_groups;
    double mse = anova_result.within.ms;
    double df_error = anova_result.within.df;
 
    double n_comparisons = static_cast<double>(k * (k - 1) / 2);
    double alpha_adj = alpha / n_comparisons;
    double t_crit = t_quantile(1.0 - alpha_adj / 2.0, df_error);
 
    std::vector<posthoc_comparison> comparisons;
 
    for (std::size_t i = 0; i < k; ++i) {
        for (std::size_t j = i + 1; j < k; ++j) {
            double mean_diff = anova_result.group_means[i] - anova_result.group_means[j];
            double n_i = static_cast<double>(anova_result.group_sizes[i]);
            double n_j = static_cast<double>(anova_result.group_sizes[j]);
 
            double se = std::sqrt(mse * (1.0 / n_i + 1.0 / n_j));
 
            double t_stat, p_value, lower, upper;
            bool significant;
 
            if (se == 0.0) {
                // Degenerate case: within-group variance is zero
                if (mean_diff == 0.0) {
                    t_stat = 0.0; p_value = 1.0;
                    lower = 0.0; upper = 0.0; significant = false;
                } else {
                    t_stat = std::copysign(std::numeric_limits<double>::infinity(), mean_diff);
                    p_value = 0.0;
                    lower = mean_diff; upper = mean_diff; significant = true;
                }
            } else {
                t_stat = mean_diff / se;
 
                p_value = 2.0 * (1.0 - t_cdf(std::abs(t_stat), df_error));
                p_value = std::min(1.0, p_value * n_comparisons);
 
                double margin = t_crit * se;
                lower = mean_diff - margin;
                upper = mean_diff + margin;
                significant = (std::abs(t_stat) > t_crit);
            }
 
            comparisons.push_back({i, j, mean_diff, se, t_stat, p_value, lower, upper, significant});
        }
    }
 
    return {"Bonferroni", comparisons, alpha, mse, df_error};
}
 
inline posthoc_result dunnett_posthoc(const one_way_anova_result& anova_result,
                                       std::size_t control_group = 0,
                                       double alpha = 0.05)
{
    if (alpha <= 0.0 || alpha >= 1.0) {
        throw std::invalid_argument("statcpp::dunnett_posthoc: alpha must be in (0, 1)");
    }
    if (control_group >= anova_result.n_groups) {
        throw std::invalid_argument("statcpp::dunnett_posthoc: invalid control group index");
    }
 
    std::size_t k = anova_result.n_groups;
    double mse = anova_result.within.ms;
    double df_error = anova_result.within.df;
 
    double n_comparisons = static_cast<double>(k - 1);
    double alpha_adj = alpha / n_comparisons;  // Bonferroni approximation
    double t_crit = t_quantile(1.0 - alpha_adj / 2.0, df_error);
 
    std::vector<posthoc_comparison> comparisons;
 
    for (std::size_t i = 0; i < k; ++i) {
        if (i == control_group) continue;
 
        double mean_diff = anova_result.group_means[i] - anova_result.group_means[control_group];
        double n_i = static_cast<double>(anova_result.group_sizes[i]);
        double n_c = static_cast<double>(anova_result.group_sizes[control_group]);
 
        double se = std::sqrt(mse * (1.0 / n_i + 1.0 / n_c));
 
        double t_stat, p_value, lower, upper;
        bool significant;
 
        if (se == 0.0) {
            // Degenerate case: within-group variance is zero
            if (mean_diff == 0.0) {
                t_stat = 0.0; p_value = 1.0;
                lower = 0.0; upper = 0.0; significant = false;
            } else {
                t_stat = std::copysign(std::numeric_limits<double>::infinity(), mean_diff);
                p_value = 0.0;
                lower = mean_diff; upper = mean_diff; significant = true;
            }
        } else {
            t_stat = mean_diff / se;
 
            p_value = 2.0 * (1.0 - t_cdf(std::abs(t_stat), df_error));
            p_value = std::min(1.0, p_value * n_comparisons);
 
            double margin = t_crit * se;
            lower = mean_diff - margin;
            upper = mean_diff + margin;
            significant = (std::abs(t_stat) > t_crit);
        }
 
        comparisons.push_back({i, control_group, mean_diff, se, t_stat, p_value, lower, upper, significant});
    }
 
    return {"Dunnett (Bonferroni approximation)", comparisons, alpha, mse, df_error};
}
 
inline posthoc_result scheffe_posthoc(const one_way_anova_result& anova_result,
                                       double alpha = 0.05)
{
    if (alpha <= 0.0 || alpha >= 1.0) {
        throw std::invalid_argument("statcpp::scheffe_posthoc: alpha must be in (0, 1)");
    }
 
    std::size_t k = anova_result.n_groups;
    double mse = anova_result.within.ms;
    double df_between = anova_result.between.df;
    double df_error = anova_result.within.df;
 
    // Scheffe critical value
    double f_crit = f_quantile(1.0 - alpha, df_between, df_error);
    double scheffe_crit = std::sqrt(df_between * f_crit);
 
    std::vector<posthoc_comparison> comparisons;
 
    for (std::size_t i = 0; i < k; ++i) {
        for (std::size_t j = i + 1; j < k; ++j) {
            double mean_diff = anova_result.group_means[i] - anova_result.group_means[j];
            double n_i = static_cast<double>(anova_result.group_sizes[i]);
            double n_j = static_cast<double>(anova_result.group_sizes[j]);
 
            double se = std::sqrt(mse * (1.0 / n_i + 1.0 / n_j));
 
            double t_stat, p_value, lower, upper;
            bool significant;
 
            if (se == 0.0) {
                // Degenerate case: within-group variance is zero
                if (mean_diff == 0.0) {
                    t_stat = 0.0; p_value = 1.0;
                    lower = 0.0; upper = 0.0; significant = false;
                } else {
                    t_stat = std::copysign(std::numeric_limits<double>::infinity(), mean_diff);
                    p_value = 0.0;
                    lower = mean_diff; upper = mean_diff; significant = true;
                }
            } else {
                t_stat = mean_diff / se;
 
                // Scheffe's F statistic
                double f_stat = (t_stat * t_stat) / df_between;
                p_value = 1.0 - f_cdf(f_stat, df_between, df_error);
 
                double margin = scheffe_crit * se;
                lower = mean_diff - margin;
                upper = mean_diff + margin;
                significant = (std::abs(t_stat) > scheffe_crit);
            }
 
            comparisons.push_back({i, j, mean_diff, se, t_stat, p_value, lower, upper, significant});
        }
    }
 
    return {"Scheffe", comparisons, alpha, mse, df_error};
}
 
// ============================================================================
// ANCOVA (Analysis of Covariance)
// ============================================================================
 
struct ancova_result {
    double ss_covariate;        
    double ss_treatment;        
    double ss_error;            
    double df_covariate;        
    double df_treatment;        
    double df_error;            
    double ms_covariate;        
    double ms_treatment;        
    double ms_error;            
    double f_covariate;         
    double f_treatment;         
    double p_covariate;         
    double p_treatment;         
    std::vector<double> adjusted_means;  
};
 
inline ancova_result one_way_ancova(
    const std::vector<std::vector<std::pair<double, double>>>& groups)
{
    std::size_t k = groups.size();
    if (k < 2) {
        throw std::invalid_argument("statcpp::one_way_ancova: need at least 2 groups");
    }
 
    // Organize data
    std::size_t n_total = 0;
    for (const auto& g : groups) {
        if (g.empty()) {
            throw std::invalid_argument("statcpp::one_way_ancova: empty group detected");
        }
        n_total += g.size();
    }
 
    if (n_total <= k + 1) {
        throw std::invalid_argument("statcpp::one_way_ancova: insufficient observations");
    }
 
    // Calculate overall means
    double sum_y = 0.0, sum_x = 0.0;
    for (const auto& g : groups) {
        for (const auto& pair : g) {
            sum_y += pair.first;
            sum_x += pair.second;
        }
    }
    double grand_mean_y = sum_y / static_cast<double>(n_total);
    double grand_mean_x = sum_x / static_cast<double>(n_total);
 
    // Calculate means for each group
    std::vector<double> group_mean_y(k);
    std::vector<double> group_mean_x(k);
    std::vector<std::size_t> group_sizes(k);
 
    for (std::size_t i = 0; i < k; ++i) {
        group_sizes[i] = groups[i].size();
        double sy = 0.0, sx = 0.0;
        for (const auto& pair : groups[i]) {
            sy += pair.first;
            sx += pair.second;
        }
        group_mean_y[i] = sy / static_cast<double>(group_sizes[i]);
        group_mean_x[i] = sx / static_cast<double>(group_sizes[i]);
    }
 
    // Calculate total sums of squares and cross-products
    double sst_y = 0.0;     // Total sum of squares (y)
    double sst_x = 0.0;     // Total sum of squares (x)
    double spt = 0.0;       // Total sum of cross-products
    double ssw_y = 0.0;     // Within-group sum of squares (y)
    double ssw_x = 0.0;     // Within-group sum of squares (x)
    double spw = 0.0;       // Within-group sum of cross-products
 
    for (std::size_t i = 0; i < k; ++i) {
        for (const auto& pair : groups[i]) {
            double dy_t = pair.first - grand_mean_y;
            double dx_t = pair.second - grand_mean_x;
            sst_y += dy_t * dy_t;
            sst_x += dx_t * dx_t;
            spt += dy_t * dx_t;
 
            double dy_w = pair.first - group_mean_y[i];
            double dx_w = pair.second - group_mean_x[i];
            ssw_y += dy_w * dy_w;
            ssw_x += dx_w * dx_w;
            spw += dy_w * dx_w;
        }
    }
 
    // Common regression coefficient (within groups)
    double b_within = (ssw_x > 0.0) ? spw / ssw_x : 0.0;
 
    // Calculate sum of squares
    double ss_error = ssw_y - b_within * spw;  // Adjusted within-group sum of squares
    double ss_covariate = b_within * spw;       // Variance explained by covariate
 
    // Overall regression coefficient
    double b_total = (sst_x > 0.0) ? spt / sst_x : 0.0;
    double ss_total_adj = sst_y - b_total * spt;
 
    // Sum of squares for treatment effect
    double ss_treatment = ss_total_adj - ss_error;
 
    // Degrees of freedom
    double df_covariate = 1.0;
    double df_treatment = static_cast<double>(k - 1);
    double df_error = static_cast<double>(n_total - k - 1);
 
    // Mean squares
    double ms_covariate = ss_covariate / df_covariate;
    double ms_treatment = ss_treatment / df_treatment;
    double ms_error = ss_error / df_error;
 
    // F statistics and p-values (handle degenerate case ms_error == 0)
    double f_covariate, f_treatment;
    double p_covariate, p_treatment;
    if (ms_error == 0.0) {
        f_covariate = (ms_covariate == 0.0) ? 0.0 : std::numeric_limits<double>::infinity();
        f_treatment = (ms_treatment == 0.0) ? 0.0 : std::numeric_limits<double>::infinity();
        p_covariate = (ms_covariate == 0.0) ? 1.0 : 0.0;
        p_treatment = (ms_treatment == 0.0) ? 1.0 : 0.0;
    } else {
        f_covariate = ms_covariate / ms_error;
        f_treatment = ms_treatment / ms_error;
        p_covariate = 1.0 - f_cdf(f_covariate, df_covariate, df_error);
        p_treatment = 1.0 - f_cdf(f_treatment, df_treatment, df_error);
    }
 
    // Adjusted means
    std::vector<double> adjusted_means(k);
    for (std::size_t i = 0; i < k; ++i) {
        adjusted_means[i] = group_mean_y[i] - b_within * (group_mean_x[i] - grand_mean_x);
    }
 
    return {
        ss_covariate, ss_treatment, ss_error,
        df_covariate, df_treatment, df_error,
        ms_covariate, ms_treatment, ms_error,
        f_covariate, f_treatment,
        p_covariate, p_treatment,
        adjusted_means
    };
}
 
// ============================================================================
// Effect Size for ANOVA
// ============================================================================
 
inline double eta_squared(const one_way_anova_result& result)
{
    if (result.ss_total == 0.0) {
        return 0.0;
    }
    return result.between.ss / result.ss_total;
}
 
inline double partial_eta_squared_a(const two_way_anova_result& result)
{
    double denom = result.factor_a.ss + result.error.ss;
    if (denom == 0.0) {
        return 0.0;
    }
    return result.factor_a.ss / denom;
}
 
inline double partial_eta_squared_b(const two_way_anova_result& result)
{
    double denom = result.factor_b.ss + result.error.ss;
    if (denom == 0.0) {
        return 0.0;
    }
    return result.factor_b.ss / denom;
}
 
inline double partial_eta_squared_interaction(const two_way_anova_result& result)
{
    double denom = result.interaction.ss + result.error.ss;
    if (denom == 0.0) {
        return 0.0;
    }
    return result.interaction.ss / denom;
}
 
inline double omega_squared(const one_way_anova_result& result)
{
    double ss_between = result.between.ss;
    double df_between = result.between.df;
    double ms_within = result.within.ms;
    double ss_total = result.ss_total;
 
    double denom = ss_total + ms_within;
    if (denom == 0.0) {
        return 0.0;
    }
    return (ss_between - df_between * ms_within) / denom;
}
 
inline double cohens_f(const one_way_anova_result& result)
{
    double eta_sq = eta_squared(result);
    if (eta_sq >= 1.0) {
        return std::numeric_limits<double>::infinity();
    }
    return std::sqrt(eta_sq / (1.0 - eta_sq));
}
 
} // namespace statcpp