FDR（False Discovery Rate）矫正

FDR（False Discovery Rate）矫正是一种用于多重假设检验的统计方法，旨在控制错误发现率，即错误地将假设测试判定为显著的比例。与Bonferroni矫正相比，FDR矫正通常更加强大，因为它尝试平衡发现真实效应的能力和控制错误发现的风险。

FDR矫正的概念

• 错误发现率（FDR）：在多次假设检验中，FDR是错误拒绝的零假设数（即假阳性）与拒绝的总假设数（即所有判定为显著的结果数）的比例的期望值。
• q值：每项测试的FDR相对应的p值，可以理解为在给定的假设检验中，该特定发现被错误地认定为显著的最大预期比例。

FDR矫正的方法

最常用的FDR控制方法是Benjamini-Hochberg（B-H）步骤上升程序。这个过程如下：

​​​​​​​1..将一系列的p值按照从大到小排序，然后利用下述公式计算每个p值所对应的FDR值。
公式：p * (n/i)， p是这一次检验的pvalue，n是检验的次数，i是排序后的位置ID（如最大的P值的i值肯定为n，第二大则是n-1，依次至最小为1）。
2.将计算出来的FDR值赋予给排序后的p值，如果某一个p值所对应的FDR值大于前一位p值（排序的前一位）所对应的FDR值，则放弃公式计算出来的FDR值，选用与它前一位相同的值。因此会产生连续相同FDR值的现象；反之则保留计算的FDR值。
3. 将FDR值按照最初始的p值的顺序进行重新排序，返回结果。
最后我们就可以使用校正后的P值进行后续的分析了。

FDR矫正的优势

• 平衡性：相对于更保守的Bonferroni矫正，FDR矫正在控制第一类错误和保持统计检验力之间提供了更好的平衡。
• 适应性：FDR矫正可以适应不同的科研场景，特别是当测试数量很大时。
• 灵活性：研究者可以根据自己的需要调整FDR阈值，以便在探索性数据分析中发现可能的效应。

应用

FDR矫正在基因组学、转录组学以及其他涉及大规模并行测试的生物信息学研究中特别有用。在这些领域，研究者常常需要在成千上万的变量中寻找显著的信号，而FDR矫正提供了一种控制错误发现的同时不过度减少真实发现的方法。

import numpy as np

# 给定的原始p值数组，包含空值
p_values = np.array([1.17836813e-01, np.nan, 1.45997443e-01, 1.37996794e-01,
                     3.38666687e-01, 2.64082940e-02, 7.47437478e-01, 4.84377465e-01,
                     8.72211809e-01, 7.38586652e-01, 7.16772005e-01, 2.05507199e-01,
                     4.49812094e-01, 3.99140973e-01, 1.20402197e-01, np.nan,
                     9.16996311e-01, 4.45464679e-01, 2.26163837e-01, 2.47342036e-01,
                     5.38874611e-01, 6.45708342e-01, 2.38058962e-01, 4.47739171e-01,
                     8.43079359e-01, 6.22000000e-05, 6.47445900e-03, 5.81733596e-01,
                     9.40436037e-01, 2.09028721e-01, 3.34373287e-01, np.nan,
                     8.02464418e-01, 6.02054444e-01, 8.26793611e-01, 9.42701410e-01,
                     3.35592785e-01, 6.56212501e-01, 5.23566380e-01, 4.42976900e-03,
                     1.08784900e-03, 4.77160630e-02, 3.04027928e-01, 2.25743353e-01,
                     9.45493857e-01, 1.04490171e-01, 8.11697045e-01, 4.53856410e-02,
                     2.09128709e-01, 8.22943091e-01, 1.94786109e-01, 2.16325134e-01,
                     1.54625825e-01, 6.31615146e-01, 2.01431377e-01, 9.62574050e-02,
                     8.65825612e-01, 1.65415505e-01, 9.90914069e-01, 8.81338110e-02,
                     2.66931120e-02, 1.88121800e-01, 5.49984100e-03, 5.12905425e-01,
                     5.50800253e-01, np.nan, 1.22319634e-01, 9.64692255e-01])

def benjamini_hochberg_correction_with_nan(p_values, alpha=0.05):
    """
    对可能含有空值的p值数组进行Benjamini-Hochberg FDR校正。
    """
    n = len(p_values)
    non_nan_p_values = p_values[~np.isnan(p_values)]
    sorted_indices = np.argsort(non_nan_p_values)
    sorted_p_values = non_nan_p_values[sorted_indices]
    
    m = len(non_nan_p_values)  # 实际参与计算的p值数量
    corrected_p_thresholds = alpha * (np.arange(1, m+1) / m)
    
    corrected_p_values = np.empty(n)
    corrected_p_values[:] = np.nan  # 初始化全部为nan
    corrected_p_values[~np.isnan(p_values)] = np.minimum.accumulate(sorted_p_values[::-1] / m * (np.arange(m, 0, -1)))[::-1][np.argsort(sorted_indices)]
    
    significant = np.zeros(n, dtype=bool)
    significant[~np.isnan(p_values)] = corrected_p_values[~np.isnan(p_values)] < alpha
    
    return corrected_p_values, significant

corrected_p_values, significant_mask = benjamini_hochberg_correction_with_nan(p_values)

print("校正后的p值: ", corrected_p_values)
print("显著性测试结果: ", significant_mask)

解释： 

 关于讲解FDR矫正的视频的一个链接如下，如果不太理解再看看：

【fMRI小知识:什么是多重比较校正？怎么计算FDR,FWE?一条死鱼的脑激活】 https://www.bilibili.com/video/BV1Yh411Z75c/?share_source=copy_web&vd_source=e155425fabb4ded4ffb3c3800d48043a