sts_correction.cpp

Go to the documentation of this file.

//=============================================================================================================
 
//=============================================================================================================
// INCLUDES
//=============================================================================================================
 
#include "sts_correction.h"
 
#include <algorithm>
#include <numeric>
#include <vector>
#include <cmath>
 
//=============================================================================================================
// USED NAMESPACES
//=============================================================================================================
 
using namespace STSLIB;
using namespace Eigen;
 
//=============================================================================================================
// DEFINE METHODS
//=============================================================================================================
 
MatrixXd StatsMcCorrection::bonferroni(const MatrixXd& pValues)
{
    const int n = static_cast<int>(pValues.size());
    MatrixXd corrected = pValues * static_cast<double>(n);
    corrected = corrected.cwiseMin(1.0);
    return corrected;
}
 
//=============================================================================================================
 
MatrixXd StatsMcCorrection::holmBonferroni(const MatrixXd& pValues)
{
    const int n = static_cast<int>(pValues.size());
 
    // Flatten to a vector with indices
    std::vector<std::pair<double, int>> indexed(n);
    for (int i = 0; i < n; ++i) {
        indexed[i] = {pValues.data()[i], i};
    }
 
    // Sort ascending by p-value
    std::sort(indexed.begin(), indexed.end());
 
    // Compute adjusted p-values: adjusted_p[i] = p * (n - rank_i + 1), rank is 1-based
    std::vector<double> adjusted(n);
    for (int i = 0; i < n; ++i) {
        adjusted[i] = indexed[i].first * static_cast<double>(n - i);
    }
 
    // Enforce monotonicity: from first to last, corrected[i] = max(corrected[i], corrected[i-1])
    for (int i = 1; i < n; ++i) {
        adjusted[i] = std::max(adjusted[i], adjusted[i - 1]);
    }
 
    // Cap at 1.0 and place back in original order
    MatrixXd corrected(pValues.rows(), pValues.cols());
    for (int i = 0; i < n; ++i) {
        corrected.data()[indexed[i].second] = std::min(adjusted[i], 1.0);
    }
 
    return corrected;
}
 
//=============================================================================================================
 
MatrixXd StatsMcCorrection::fdr(const MatrixXd& pValues, double alpha)
{
    Q_UNUSED(alpha);
 
    const int n = static_cast<int>(pValues.size());
 
    // Flatten to a vector with indices
    std::vector<std::pair<double, int>> indexed(n);
    for (int i = 0; i < n; ++i) {
        indexed[i] = {pValues.data()[i], i};
    }
 
    // Sort ascending by p-value
    std::sort(indexed.begin(), indexed.end());
 
    // Compute adjusted p-values: adjusted_p = p * n / rank (rank is 1-based)
    std::vector<double> adjusted(n);
    for (int i = 0; i < n; ++i) {
        int rank = i + 1;
        adjusted[i] = indexed[i].first * static_cast<double>(n) / static_cast<double>(rank);
    }
 
    // Enforce monotonicity: from last to first, corrected[i] = min(corrected[i], corrected[i+1])
    for (int i = n - 2; i >= 0; --i) {
        adjusted[i] = std::min(adjusted[i], adjusted[i + 1]);
    }
 
    // Cap at 1.0 and place back in original order
    MatrixXd corrected(pValues.rows(), pValues.cols());
    for (int i = 0; i < n; ++i) {
        corrected.data()[indexed[i].second] = std::min(adjusted[i], 1.0);
    }
 
    return corrected;
}