Python SciPy stats.false_discovery_control用法及代碼示例

本文簡要介紹 python 語言中 scipy.stats.false_discovery_control 的用法。

用法:

scipy.stats.false_discovery_control(ps, *, axis=0, method='bh')#

調整 p 值以控製錯誤發現率。

錯誤發現率 (FDR) 是被拒絕的原假設實際上為真的預期比例。如果當調整後的 p 值低於指定水平時拒絕原假設，則錯誤發現率將控製在該水平。

參數 ：：

ps： 一維數組

要調整的 p 值。元素必須是 0 到 1 之間的實數。

axis： int

執行調整所沿的軸。調整是沿著每個axis-slice獨立執行的。如果axis為None，則在執行調整之前先清除ps。

method： {‘bh’, ‘by’}

要應用的錯誤發現速率控製過程：'bh'用於Benjamini-Hochberg[1](等式1)，'by'用於Benjaminini-Yekutieli[2](定理1.3)。後者更為保守，但即使 p 值不是來自獨立測試，也能保證控製 FDR。

返回 ：：

ps_adusted： array_like

調整後的 p 值。如果零假設在低於指定水平的情況下被拒絕，則錯誤發現率將被控製在該水平。

注意：

在多重假設檢驗中，錯誤發現控製程序往往比係列錯誤率控製程序(例如 Bonferroni 校正 [1])提供更高的功效。

如果 p 值與獨立測試相對應(或使用 “positive regression dependencies” [2] 進行測試)，則拒絕與低於 \(q\) 的 Benjamini-Hochberg-adjusted p 值相對應的零假設可控製錯誤發現率，其水平小於或等於\(q m_0 / m\) ，其中 \(m_0\) 是真零假設的數量，\(m\) 是測試的零假設的總數。當根據更保守的 Benjaminini-Yekutieli 過程調整 p 值時，即使對於相關檢驗也是如此。

該函數生成的調整後的 p 值與 R 函數生成的 p 值相當p.adjust和 statsmodels 函數statsmodels.stats.multitest.multipletests。請考慮後者以獲得更高級的多重比較校正方法。

參考：

[1] (1,2,3,4,5)

本傑明尼、約阿夫和約瑟夫·霍赫伯格。 “控製錯誤發現率：一種實用而強大的多重測試方法。” 《皇家統計學會雜誌》：B 係列(方法論)57.1 (1995)：289-300。

[2] (1,2)

本傑明尼、約夫和丹尼爾·耶庫蒂利。 “依賴關係下多重測試中錯誤發現率的控製。”統計年鑒(2001)：1165-1188。

[3]

瓷磚統計。羅斯福 - Benjamini-Hochberg 解釋 - YouTube。 https://www.youtube.com/watch?v=rZKa4tW2NKs。

[4]

諾伊豪斯，Karl-Ludwig，等人。 “front-loaded 給予阿替普酶可改善急性心肌梗塞的溶栓：rt-PA-APSAC 通暢研究 (TAPS) 的結果。”美國心髒病學會雜誌 19.5 (1992)：885-891。

例子：

我們遵循[1]中的示例。

Thrombolysis with recombinant tissue-type plasminogen activator (rt-PA) and anisoylated plasminogen streptokinase activator (APSAC) in myocardial infarction has been proved to reduce mortality. [4] investigated the effects of a new front-loaded administration of rt-PA versus those obtained with a standard regimen of APSAC, in a randomized multicentre trial in 421 patients with acute myocardial infarction.

該研究測試了四個假設係列，最後一個是“血栓治療開始後的心髒和其他事件”。在這一係列假設中可能需要 FDR 控製，因為如果 front-loaded 治療僅相當於先前的治療，則不適合得出 front-loaded 治療更好的結論。

與該族中 15 個假設相對應的 p 值為

>>> ps = [0.0001, 0.0004, 0.0019, 0.0095, 0.0201, 0.0278, 0.0298, 0.0344,
...       0.0459, 0.3240, 0.4262, 0.5719, 0.6528, 0.7590, 1.000]

如果所選的顯著性水平為 0.05，我們可能會傾向於拒絕與前 9 個 p 值對應的檢驗的原假設，因為前 9 個 p 值低於所選的顯著性水平。然而，這會忽略“multiplicity”的問題：如果我們未能糾正正在執行多重比較的事實，我們更有可能錯誤地拒絕真零假設。

解決多重性問題的一種方法是控製 family-wise 錯誤率 (FWER)，即原假設在實際為真時被拒絕的比率。這種常見的過程是 Bonferroni 校正 [1]。我們首先將 p 值乘以測試假設的數量。

>>> import numpy as np
>>> np.array(ps) * len(ps)
array([1.5000e-03, 6.0000e-03, 2.8500e-02, 1.4250e-01, 3.0150e-01,
       4.1700e-01, 4.4700e-01, 5.1600e-01, 6.8850e-01, 4.8600e+00,
       6.3930e+00, 8.5785e+00, 9.7920e+00, 1.1385e+01, 1.5000e+01])

為了將 FWER 控製在 5%，我們僅拒絕與調整後的 p 值小於 0.05 相對應的假設。在這種情況下，隻能拒絕與前三個 p 值相對應的假設。根據[1]，這三個假設涉及“allergic reaction”和“出血的兩個不同方麵”。

另一種方法是控製錯誤發現率：被拒絕的原假設實際上為真的預期比例。這種方法的優點在於，它通常提供更大的功效：當原假設確實為假時，拒絕原假設的概率更高。為了將錯誤發現率控製在 5%，我們應用了 Benjamini-Hochberg p 值調整。

>>> from scipy import stats
>>> stats.false_discovery_control(ps)
array([0.0015    , 0.003     , 0.0095    , 0.035625  , 0.0603    ,
       0.06385714, 0.06385714, 0.0645    , 0.0765    , 0.486     ,
       0.58118182, 0.714875  , 0.75323077, 0.81321429, 1.        ])

現在，前四個調整後的 p 值低於 0.05，因此我們將拒絕與這四個 p 值相對應的原假設。拒絕第四個零假設對於最初的研究尤其重要，因為它得出的結論是新療法“大大降低了 in-hospital 死亡率”。