Micorarray 分析：P-value 的 Multiple test correction and False Discovery Rate

p-value 在統計上通常代表虛無假設（null hypothesis）為真的機率，用以顯示判斷的正確性，有人會翻譯成信心水準。用個例子來說明， 假如我們假設箱子中有十粒球 全部是白色的(hypothesis)，那麼箱子裡面至少有一顆其它顏色的球就是虛無假設(null hypothesis)。然後我們開始做實驗，一次抽一顆，記下顏色後放回去，連抽了十次，全是白色的，所以你可以這樣說，我抽了十次，都是白球，所以箱子裡都是白球。但這麼說也可能是錯的，也許你只是剛好沒抽到那顆黑球。我們可以由 null hypothesis 來推論一下錯誤的機率，如果箱子中有一粒黑球，連抽十次都沒有抽到的機率是(9/10)的10次方(Ｐ約為 0.35)，所以剛剛那個箱子裡全是白球的推論有0.35的機會錯誤,換句話說信心水準為0.35。

這是一個極度簡化的實驗，在現實世界中遠比這個例子複雜，除了取樣數目的問題外，還有各種環境及個體差異的因子存在，像是白球裡可能會滲了一粒被弄黑的白球等等。不過，這些問題可以透過統計的檢定，用不同的機率的分布來推測null hypothesis 的機率，科學界也普遍的使用p value 來描述結果的可信度，你常會看到說癌細胞A gene 表現量比正常細胞高 (p <0.001)，這就表示你可能測個上千個樣品，才可能會有一個例外。

如果是單一實驗來看 p-value 是沒問題的，但你如果一次要檢查一千個箱子事情可能更複雜。如果這一千個箱子中只有三個真得全是白球。而你一如往常的用在＂p-Value < 0.05 的信心水準下＂,你相信箱子實驗的結果。這時就會有很嚴重的問題了。在這樣的標準下，就算每個箱子都有黑球，仍會有五十個箱子(1000 x 0.05)被判定全是白球（偽真 False positive)，這遠比真正全是白球的箱子（true positive)還多....


##ReadMore##
因為這個原因，所有在大量平行檢定統計時（例如：Microarray 實驗)，對p-value 的要求較單一實驗又更嚴苛了， 我們會做額外的 p-Value 調整，以減少 false positive 的可能。這類的方法叫 Multiple test correction 或False Discovery Rate. 校正的方法大致可以分成三類:：

1.　Multiple test correction (MTC)

 MTC 算是比較嚴謹的方法，它的精神在於控制 Experiment-wise error Type I rate，也就是說上述例子中，MTC 的修正p-value一個代表讓整個實驗（所有一千個箱子的實驗）總錯誤機率(ex. 1000次的 total error rate <0.05)，常用的方法有： Bonferroni, Dunn-Sidak, Bonferroni Step-down (Holm)，大概說一下原理，如果想知道更詳細請參閱 <Bonferroni and Dunn Sidak by R>


首先,假定Microarry上測試基因數目共有k 個。整個實驗的信心強度是 (Experiment-wise significant level) αe，單一個基因的測試信心強度是（comparison-wise significance level）αc

Bonferroni correction
Bonferroni correction 是最直覺的校正法，他的做法很簡單，就是把檢定的 p-Value (t-test, ANOVA....)標準乘以測試的次數。例如用一個有1000個基因的Microrray做基因表現變異分析中，我們用t-test來找出兩樣品群有差異的基因，p-value 小於 0.05判定為有顯著差異，但因為我們一次做了1000次，這個p-value 的門檻可能會有很多false positive的基因，就整個實驗而言可能出現的錯誤 0.05 x 1000 = 50, 超過1,為了把整個j實驗的錯誤率修正回 0.05, 最簡單作法就是把p-value 的門檻除以 1000, (p-value<0.00005)，可能出現的錯誤 0.00005 x 1000 = 0.05。這個方法很直覺，但會發生的問題也可以從這個例子中看的出來，因為大部份Microarray 實驗能符合p-value<0.00005　的基因，用手指數的出來，剩下來的99.99% 都被丟棄了。Bonferroni correction太過嚴謹，所以不是很普遍，特別是在歧異度高的臨床研究。

Bonferroni correction
K independent test, so all α is equal. so overall number of type I errors = K x α
αe = αc/K => αc = αe x K

Correct p = p x K

Dunn-Sidak
Dunn-Sidak 的目標與 Bonferroni 一樣，希望控制整個實驗的錯誤率，只是算法有點不同，
Dunn-Sidak 假定單一實驗的信心強度是 αc時， 那麼這個實驗正確的機率為 1-αc，那麼整個實驗(k 個平行分析）的錯誤率（即信心強度 αe） 就等於一減去所有單次數驗都是正確的機率 (1-αc)的K次方.

αe = 1-(1-αc)^K

交換一下左右象，也就是說

αc = 1-(1-αe)^1/k

所以校正的 p 值等於 1-(1-p)^1/k

Dunn-Sidak 的作法是假設完全不能有false positive的，這對一個次執行上萬平行基因檢測的實驗太過嚴苛。所以跟Bonferroni 一樣，有著曲高合寡的問題。


Bonferroni Step-down (Holm)
為了解決太過嚴謹的問題，Bonferroni Step-down方法作了一點讓步來。Bonferroni Step-down 的作法是先以 p-value 由小到大排序，然後逐個校正如下：

Smallest

1st. corrected p-value = p-value x n(gene number)

2nd. corrected p-value = p-value x n-1

3rd. corrected p-value = p-value x n-2

....

nth. corrected p-value = p-value x 1

Biggest

2.False positive discovery (FDR)
跟MTC不同，FDR不是校正p-value 的方法，FDR是用來描述以一個p-value 為標準所判定的結果，有多少是False positive的機會，也就是所謂Family-wise error rate。例如一千個箱子以p-value<0.05為信心準來判定箱內全是白球，FDR 算出來p-value <0.05 的 q-value = 0.2，這就表示，所有判定全為白球的箱中，存在20%的箱子可能是錯的。常用的方法有 qValue (Storey)，False Discovery Rate，False Discovery Rate 又可以分Step up(step up, Benjamini & Hochberg, 1995) , Step down　(Benjamini & Liu, 1999),計算方法如下。

Step up False Discovery Rate
(參考資料來源）
假設我們同時做了 m 次實驗（eg. 例如Microarray 上的 m 個點），就會算出 m 個 p-value接下來我們按照 p 值昇冪排序 ( i = 1~m)，第 i 個 p-value 調整為：

簡單的說明這個公式，裡面那個 min( ) 是為了確保 P 不會大於１，外面的min{} 是從第 i 個算到最後一個 m (k=i~m, 所以叫 step up), 取最小值．
p-Value sort ascending
1st. Corrected p = p x n/(n-0)
2nd. Corrected p = p x n/(n-1)
3rd. Corrected p = p x n/(n-2)
......
nth. Corrected p = p x n /1


Step down False Discovery Rate
(參考資料來源）
Step down 跟前項類似，只是 p 值調整的方法不同，公式如下列：

這個公式 P值算法不同，外面的min{} 是反過來從第 i個算到第一個 i (k=1~i, 所以叫 step Down), 取最小值．

1st. Corrected p = p x n/(n)
2nd. Corrected p = p x n/(n+1-2)
3rd. Corrected p = p x n/(n+1-3)
......
nth. Corrected p = p x n /1

(2013/5/15 ：感謝 AlvinLu 的指正，修正了 step-up, step-down 的計算方式，我不是統計領域的人，閱讀這些符號對我來說也是很痛苦事情，只能儘可能的就我理解範圍整理收集到的資訊，如果有錯誤，歡迎各位專家指正，我虛心接受並隨時修改）

q-Value
qValue (positive false discovery rate) 的計算方法比較複雜，他是利用 bayesian 分布去推估的，實際上的算法，我研究了很久，但是最後還是看不懂，隔行如隔山呀！不過，q-Value 在array 分析還蠻常用的，很多軟體都有支援，不過我比較常用Ｒ來算。



3.Boot strap
Boot strap (Westfall and Young Permutation)是很不一樣的方法,他撇開數值方法，直接去作實驗來驗証，例如我們用t-test測試兩群族的A gene 表現有差異，p-Value = 0.01，為了確定這個差異真的出現在兩個樣品群中，我們就隨機的把樣品分成兩個組別，如果用隨機分組的t-test都顯著程度皆無法超越真正的樣品分組，那我們就確認 這個 p-value 可能是真的，不是偶發的 false positive, 這種隨機分組方法有時候在缺乏資料模型時很好用，不過缺點是非常的慢，而其 false positive rate 相近放前Bonferroni Step-down方法。


信心水準要設多少才可信?

介紹完三種p-Value的校正問題,接下來最常被問到的問題是:那麼信心水準要設多少才可信?老實說,並沒有標準答案,我們當然門檻越高越好,但如果花了幾十萬作實驗,到頭來一個基因也選不到,相信你也會不甘心吧!務實一點,我通常是以數量來選基因,如果做Biomark我會選前一百吧!功能註解富集分析,我大概會選一千個以內,至於pValue則是用來提醒自已,以這些資料為基礎的分析結果,可信度有多少,如果是薄弱的,在後續就需要的更多的實驗證明。P-Value畢竟是個形容性的數字,不必拘泥一個固定的數字。