如何優(yōu)雅地校正多重檢驗的p值？

Jessiee_ 2018-09-24

展開全文

組學(xué)派的大家庭越來越熱鬧，派派也是干勁滿滿！作為一個公益服務(wù)性的平臺，為組學(xué)戰(zhàn)友們排憂解難始終是派派小組至高無上的責(zé)任和追求！更加要慶賀的是加入派派的小伙伴越來多，我們有了總?cè)海€有了總?cè)?號（二維碼見文末），還有各個分群。我們有了新的LOGO，相信omicsPie會越來越好！

上周，群里酷愛科研的大佬們對差異性分析結(jié)果p value的校正問題掀起了熱情洋溢的民主大討論。既然問題出現(xiàn)了，那么這寶貴的科普機會派派又怎么能夠放過，今天就讓苯寶寶系統(tǒng)地為大家梳理一下有關(guān)False Positive（假陽性）的前生今世，為后入的小伙伴們點亮這盞指坑。。。啊，不，脫坑明燈。

一、什么是False positive（假陽性）？

提及假陽性，不得不從統(tǒng)計假設(shè)檢驗的兩類錯誤講起。第一類錯誤（Type I error），就是原假設(shè)H0成立（也就是檢驗的兩組之間沒有差別），但是統(tǒng)計結(jié)果卻顯示拒絕原假設(shè)（可以理解為結(jié)果p<0.05），這樣就得到了“假”的“陽性”結(jié)果。相應(yīng)的第二類錯誤（Type II error），則是原假設(shè)H0不成立（也就是檢驗的兩組之間有差別，是不同的），但是統(tǒng)計結(jié)果卻顯示接受原假設(shè)（可以理解為結(jié)果p>0.05），這樣就得到了“假”的“陰性”結(jié)果。

今天，我們主要談一談“假陽性”。關(guān)于“假陰性”和“把握度”的問題，下次咱們再接著聊。

二、假陽性是如何產(chǎn)生的？

科學(xué)研究中，研究者往往很難得到全部的研究對象的實驗數(shù)據(jù)，所以隨機抽樣、統(tǒng)計推斷的思想貫穿于實驗設(shè)計的始終。既然樣本是從總體中抽出來的，即使是隨機化做得非常好，依然難免會存在抽樣誤差，導(dǎo)致抽樣數(shù)據(jù)的隨機波動，不能完全準確地反應(yīng)總體的數(shù)據(jù)水平。這種情況下，很顯然，即使是從同一個總體中反復(fù)抽樣，也有機會得到樣本之間的差異經(jīng)統(tǒng)計學(xué)檢驗p<0.05的結(jié)果（如圖1）。因此，不難得到這樣的結(jié)論：假陽性結(jié)果會隨著抽樣次數(shù)、比較次數(shù)的增多而大量出現(xiàn)。這里派派做了個簡單的模擬實驗，從同一組正態(tài)分布總體中隨機抽取兩組進行t-test，比較1000次累積的陽性結(jié)果（p<0.05）次數(shù)如圖2所示。

另外，如果同時分析的指標增加，假陽性率又會有怎樣的變化呢？若只分析一個指標，假陽性率為0.05，那么不犯I類錯誤的概率為95%；若同時分析兩個指標，假陽性率則變?yōu)?-（0.95）²=0.0975；若同時分析三個指標，假陽性率則變?yōu)?-（0.95）³=0.0143......若同時分析十個指標，假陽性率則變?yōu)?-（0.95）¹⁰=0.4013，假陽性率已經(jīng)超過了40%！

在組學(xué)數(shù)據(jù)的研究過程中，由于統(tǒng)計分析的變量維度較高，往往會進行指數(shù)級的假設(shè)檢驗比較，假陽性的結(jié)果則更加常見。

圖1

圖2

三、假陽性的校正方法

（一）傳統(tǒng)統(tǒng)計上的多重檢驗結(jié)果校正方法，主要用于與多組間方差分析（ANOVA）相結(jié)合的兩兩比較。簡單總結(jié)為以下幾種：

1.LSD（least significant difference）最小顯著差異t檢驗校正。這種校正方法一般用于事前比較，也就是實驗設(shè)計時已經(jīng)確定進行某些組之間的比較（如A和B，C和D），而其他組間不必進行比較。假如所有組間都要比較用LSD的話會增加假陽性的概率。

2. Dunnet-t校正。這種方法適用于多個實驗組均數(shù)與對照組均數(shù)間的比較，也就是指定其中一組如A組為對照組，B、C、D組均與A組比較，但B、C、D組之間不進行比較。

3.SNK（Student-Newman-Keuls）檢驗。此方法適用于ANOVA之后的多組間兩兩比較（事后的全局比較，不指限定比較分組），常用于探索性研究。但結(jié)果只告訴有無差異，不提供精確p值。常見的結(jié)果形式為

分組	1	2
A	mean(A)
B	mean(B)	mean(B)
C		mean(C)

A組和B組在類別上都被定義為1類，說明兩組mean值之間的差異經(jīng)校正后沒有統(tǒng)計學(xué)意義；同理B、C兩組間的差異也沒有統(tǒng)計學(xué)意義。但A、C兩組分別被定義為1類和2類，說明A、C兩組間mean值的差異具有統(tǒng)計學(xué)意義。

以上三種校正方法，使用的前提是各組均滿足正態(tài)性和方差齊性，即僅適用于與ANOVA、多重t檢驗結(jié)合使用。而另外兩種適用性更強的萬能校正方法，也備受推崇：

4.Bonferroni 校正。這種校正方法比較簡單：校正檢驗水準α'=α/m（這里α通常為0.05，m為檢驗次數(shù)），也相當(dāng)于如果保持α=0.05不變，校正p值為p'=m*p。這種校正方式最為保守嚴格，也就是得到的陽性結(jié)果把握度會很高，但靈敏度比較低，往往不適用于潛在研究目標的篩選。

5.Sidak法校正。既可用于事前比較，也可用于事后比較。校正檢驗水準α'=1-（1-α）^1/m，m為比較次數(shù)，更適用于比較次數(shù)較多時使用。

當(dāng)然，還有例Tukey法、Duncan法等諸多方法，小編研究得不多就暫不在此枚舉了。

（二）更加適用于高維數(shù)據(jù)的校正方法Benjamini-Hochberg false discovery rate (FDR)，是基于對假陽性發(fā)現(xiàn)率的控制來決定p值的閾值。

相對Bonferroni來說，F(xiàn)DR校正更加溫和。其目標是在假陽性和假陰性間達到平衡，將假/真陽性的比例控制在一定范圍之內(nèi)。例如，如果檢驗100次，我們設(shè)定的FDR閾值為0.05（5%），那么無論我們得到多少個差異特征變量，這些差異特征變量中出現(xiàn)假陽性的概率將保持在5%之內(nèi)，這就是控制FDR＜5%。這種校正方法意義明確，更加適用于新的差異物質(zhì)的發(fā)現(xiàn)，例如用FDR<5%的標準篩選到100個差異基因、蛋白或代謝物，那么可以理解為這100個物質(zhì)里面大約有95個是在兩組間存在真實差異。

FDR的計算公式也十分簡單：q=p *m/rank，這里的m為p值的總個數(shù)，rank為p值從大到小的排序。簡單來說，若進行了100次差異分析得到了100個p值（m=100），則只有最小的p值rank=100，校正之后保持不變，而其他較大的p值排序靠前<100，都會由于檢驗次數(shù)的增加而受到懲罰。

當(dāng)然，F(xiàn)DR還有很多其他的估計方法和計算公式如：SAM法、經(jīng)驗貝葉斯法等，感興趣的童鞋們可以自行查閱相關(guān)文獻，如哈醫(yī)大衛(wèi)生統(tǒng)計教研室李康教授的一篇中文綜述《多重假設(shè)檢驗中FDR的控制與估計方法》中就有更多細節(jié)的介紹。

四、FDR校正的R語言實現(xiàn)

能做FDR校正的R packages還是有很多的，派派給大家介紹一個自己常用的“fdrtool”：

1.計算p value，合并為一個向量。

如果有100個p value,則可記錄為data = c(p1:p100)

2. 載入程序包

library('fdrtool')

3.計算FDR校正后的p value

FDR <- fdrtool(data,statistic='pvalue',plot=F)$qval

4.計算結(jié)果中FDR的順序與data中p value的順序一一對應(yīng)。

五、p value校正小建議

千萬不要迷信p value！千萬不要迷信p value！千萬不要迷信p value！

無論是p value還是校正后的p value，都只是代表了統(tǒng)計學(xué)上的概率，而概率這個東西在一定程度上并不意味著最終結(jié)果的真實與否。所以從統(tǒng)計學(xué)的角度，派派更傾向于在研究設(shè)計之前就確認好研究的主要目標，圍繞主要目標進行的合理嚴謹?shù)目蒲性O(shè)計才是提高結(jié)果可靠度的最佳途徑，而在研究過程中意外發(fā)現(xiàn)的陽性結(jié)果，只有經(jīng)過更多的重復(fù)試驗，或者更有針對性的實驗設(shè)計重復(fù)研究，才可以認為是靠譜的結(jié)果。

同樣的，即使發(fā)現(xiàn)的差異物質(zhì)p value不夠顯著，但呈現(xiàn)了預(yù)期的變化趨勢，在考慮到樣本量可能不足，實驗設(shè)計的混雜因素控制得不夠理想的情況下，建議不要輕易放棄，最好進行下實驗研究的把握度評價。若本次實驗的把握度不夠理想，而從專業(yè)角度上該物質(zhì)發(fā)揮效用的可能性甚高，那么也值得重新設(shè)計或擴大實驗樣本量，進行深入的驗證和探索。