錯誤發(fā)現(xiàn)率的經(jīng)驗估計和應(yīng)用*

王 婷，曾 平，黃水平，趙華碩

1)徐州醫(yī)學(xué)院公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計學(xué)教研室徐州 221002 2)南京醫(yī)科大學(xué)公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計學(xué)教研室南京 210029

#通訊作者，男，1982年7月生，碩士，助教，研究方向:高緯數(shù)據(jù)分析和貝葉斯統(tǒng)計，E-mail:zpstat@xzmc.edu.cn

錯誤發(fā)現(xiàn)率的經(jīng)驗估計和應(yīng)用*

王 婷1)，曾 平1，2)#，黃水平1)，趙華碩1)

1)徐州醫(yī)學(xué)院公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計學(xué)教研室徐州 221002 2)南京醫(yī)科大學(xué)公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計學(xué)教研室南京 210029

#通訊作者，男，1982年7月生，碩士，助教，研究方向:高緯數(shù)據(jù)分析和貝葉斯統(tǒng)計，E-mail:zpstat@xzmc.edu.cn

微陣列數(shù)據(jù);錯誤發(fā)現(xiàn)率;經(jīng)驗貝葉斯;密度估計

目的:研究大規(guī)模數(shù)據(jù)中的密度、無效分布和錯誤發(fā)現(xiàn)率的經(jīng)驗貝葉斯估計和應(yīng)用。方法:對2個微陣列數(shù)據(jù)的貝葉斯模型，采用Poisson回歸方法估計密度函數(shù)，并在此基礎(chǔ)上經(jīng)驗估計貝葉斯錯誤發(fā)現(xiàn)率和局部錯誤發(fā)現(xiàn)率。結(jié)果:基于Poisson回歸方法的密度估計為無效分布和錯誤發(fā)現(xiàn)率的經(jīng)驗貝葉斯估計提供了恰當(dāng)?shù)姆椒ㄟx擇。結(jié)論:大規(guī)模數(shù)據(jù)的平行結(jié)構(gòu)使得對錯誤發(fā)現(xiàn)率和無效分布的估計變得可能。

錯誤發(fā)現(xiàn)率(false discovery rate，F(xiàn)DR)在大規(guī)模數(shù)據(jù)分析中起著十分重要的作用［1-2］，被越來越多地應(yīng)用在微陣列和功能磁共振成像等領(lǐng)域［3-4］。FDR可以根據(jù)尾部面積和密度定義，前者以最初提出的FDR和陽性錯誤發(fā)現(xiàn)率(positives false discovery rate，pFDR)為代表［5］，后者主要指局部FDR(local false discovery rate，locfdr)［6］。除了具有現(xiàn)實的應(yīng)用價值外，F(xiàn)DR一個吸引人的地方在于，它同時具有頻率統(tǒng)計的性質(zhì)和可解釋的貝葉斯含義［5］，此外，大規(guī)模數(shù)據(jù)使得直接估計FDR變得可能，因此FDR也具有經(jīng)驗貝葉斯的意義。Benjamini和Hochberg［1］提供了一個十分有用的控制程序，保證FDR不大于一個預(yù)先設(shè)定的水準(zhǔn)，與控制過程相反，該文主要研究在貝葉斯統(tǒng)計的框架下對FDR進(jìn)行經(jīng)驗估計，包括估計基于尾部面積的FDR和locfdr。

1 方法與原理

1.1 理論無效分布 設(shè)m個檢驗的無效假設(shè)(zeroassumption，ZA)和備擇假設(shè)為:H0i=0和H1i=1，統(tǒng)計量zi可由其他統(tǒng)計量轉(zhuǎn)換而來，例如:

Φ、Tv分別表示正態(tài)和自由度為v的t變量累計分布函數(shù)。圖1給出了2個微陣列數(shù)據(jù)z值的直方圖［7-8］，前者描述了50個正常對照和52個前列腺癌患者6 033個基因的表達(dá)水平，后者描述了45例急性淋巴細(xì)胞白血病患者和27例急性髓性白血病患者7 128個基因的表達(dá)水平。目的都是希望發(fā)現(xiàn)哪些基因在2組人群中存在表達(dá)差異。顯然如果H0i=0成立，則Z～N(0，1)，見圖1的虛線，此后稱N(0，1)為理論無效分布。圖1顯示前列腺癌數(shù)據(jù)中0附近的z值和N(0，1)比較吻合，但白血病數(shù)據(jù)中0附近的z值明顯異于N(0，1)。但后文中仍然首先假設(shè)白血病數(shù)據(jù)的理論無效分布是滿足的。

1.2 貝葉斯模型 假設(shè)所有基因只屬于差別表達(dá)和無差別表達(dá)兩種情況［9］，設(shè)H0=0的先驗概率為π0，H1=1的先驗概率為1－π0，H0=0時z值的密度和分布函數(shù)分別為f0、F0，H1=1時z值的密度和分布函數(shù)分別為f1、F1。則z值具有混合邊際密度f(z)=π0f0(z)+(1－π0)f1(z)和混合分布函數(shù)F(z)=π0F0(z)+(1－π0)F1(z)。如以Γ=(Z≤z)作為拒絕域，根據(jù)貝葉斯定理FDR(z)=π0F0(z)/F(z)和locfdr(z)=Pr(H0|Z=z)=π0f0(z)/f (z)［6］。FDR(z)可看作是后驗尾部面積，locfdr(z)可看作是后驗概率，因此從貝葉斯角度看locfdr(z)更具有可解釋性。兩者之間的關(guān)系為E(locfdr)= FDR。

圖1 前列腺癌(上)和白血病(下)數(shù)

1.3 經(jīng)驗貝葉斯 按照前面的假設(shè)，如H0=0成立，則f0(z)=N(0，1)，F(xiàn)0(z)=Φ，在微陣列數(shù)據(jù)中由于稀疏的原因，π0多大于0.90，此時即使取π0= 1也不會對FDR的估計產(chǎn)生很大的影響，因此FDR的完全貝葉斯分析只需要給f(z)或F(z)指定先驗，但是微陣列的大規(guī)模平行結(jié)構(gòu)能夠直接應(yīng)用數(shù)據(jù)估計f(z)或F(z)，由此得到FDR的經(jīng)驗貝葉斯模型。累積分布函數(shù)F最直接的估計是經(jīng)驗分布:ˉF(z)= #{Z≤z}/m。FDR的非參數(shù)經(jīng)驗貝葉斯估計值為:

FDR(z)=π0F0(z)/ˉF(z)

例如對前列腺癌數(shù)據(jù)，如果選擇Γ={Z≤－3}，取π0=1，有F0(－3)=1.35×10－3，#{Z∈Γ}=49， FDR(－3)=mF0/49=0.166，這意味著在49個差別基因中大約有1/6屬于錯誤識別。

1.5 經(jīng)驗無效分布 假設(shè)f0=N(0，1)，從圖1可見理論無效分布對前列腺癌數(shù)據(jù)是合理的，但對白血病微陣列那樣的高維數(shù)據(jù)，在這種情況有必要重新選擇更加合適的無效分布。在傳統(tǒng)的單個假設(shè)檢驗中，應(yīng)用者只能被動接受理論無效分布，但在大規(guī)模數(shù)據(jù)分析中能夠利用數(shù)據(jù)估計無效分布，稱之為經(jīng)驗無效分布［10］。如圖1，雖然不能明確哪些基因來源于H0=0或是H1=1，但是幾乎可以肯定的是，絕大多數(shù)存在于0附近的基因應(yīng)該來自H0=0，這些基因的z值可用來估計無效分布，上述的假設(shè)稱為零假設(shè)(zero assumption，ZA): f1(z)=0，z∈［－a，a］。a為一個固定的截點(diǎn)，比如a=0.5。仍然假設(shè)經(jīng)驗無效分布為正態(tài)分布，但具有不一樣的參數(shù):f0ZA=N(z|μ，σ2)。ZA暗示，如果f0ZA為正態(tài)分布，那么0附近的z值應(yīng)該和f=π0f0ZA一致，根據(jù)這個原理可采用中心匹配的方法估計參數(shù)N0、μ、σ2［11］。

2 結(jié)果

見表1。選擇a=0.5時的參數(shù)估計值，此時前列腺癌和白血病數(shù)據(jù)中#{z∈(－0.5，0.5)}分別為2 112和1 571。圖2給出了不同z值的locfdr，此時前列腺癌和白血病數(shù)據(jù)中l(wèi)ocfdr(z)≤0.2的基因數(shù)分別為42和202個，可見經(jīng)驗無效分布對白血病數(shù)據(jù)結(jié)果的影響之大，而前列腺癌數(shù)據(jù)結(jié)果的影響主要來自于π0。

表1 不同a值對應(yīng)的經(jīng)驗無效分布參數(shù)估計值

圖2 前列腺癌(上)和白血病(下)數(shù)據(jù)的locfdr

3 討論

FDR在大規(guī)模數(shù)據(jù)分析中具有十分現(xiàn)實的應(yīng)用性，雖然是在頻率統(tǒng)計下發(fā)展起來的，但同時也具有可解釋的貝葉斯和經(jīng)驗貝葉斯含義。頻率統(tǒng)計意義下的FDR和經(jīng)典的假設(shè)檢驗中基于尾部面積的統(tǒng)計決策思維是一致的，而locfdr則具有貝葉斯后驗概率的意義。Benjamini和Hochberg關(guān)于FDR的控制程序及貝葉斯FDR都是相對一個拒絕域而言的，但它們并沒有對單個 z值給予任何陳述，而locfdr基于單個統(tǒng)計量，能夠量化zi＞zj時存在的FDR差別，而這種差別可能正是研究者所關(guān)心的。

但是估計locfdr要比估計基于尾部面積的FDR更加困難，前者涉及密度估計，后者只需要估計經(jīng)驗分布即可。Poisson回歸方法為估計密度函數(shù)提供了足夠高的精度和準(zhǔn)確度，即使是在z值具有相關(guān)性時仍然能夠得到滿意的結(jié)果(微陣列數(shù)據(jù)常常存在相關(guān))。除 Poisson回歸方法外，其他方法如Grenander密度估計也可用來估計邊際密度［12］，前者的優(yōu)勢在于將密度估計轉(zhuǎn)化為了更加熟悉的回歸理論，并且能夠用來進(jìn)一步對FDR進(jìn)行光滑估計，能在常用的軟件中執(zhí)行，如R軟件的glm函數(shù)。但Poisson回歸方法中引入了額外的參數(shù)如組段和多項式(或樣條函數(shù))的自由度，對某些異常數(shù)據(jù)可能需要更加細(xì)致地選擇和調(diào)整。Efron等［6］采用了基于置換檢驗的logistic回歸，用f0(z)/f(z)的比值間接估計locfdr。

高維數(shù)據(jù)分析中另一個十分重要的問題是，理論無效分布常常被違背，如白血病數(shù)據(jù)，在這種情況FDR的估計明顯錯誤，但大規(guī)模平行數(shù)據(jù)結(jié)構(gòu)能夠?qū)o效分布進(jìn)行經(jīng)驗估計。作者選擇ZA條件下的匹配估計方法發(fā)現(xiàn)，無效分布的參數(shù)依賴選擇的固定常數(shù)a，但是在一定范圍內(nèi)a的影響有限。a越小，對ZA條件的信心越高，此時偏倚減小，但同時用于估計參數(shù)的數(shù)據(jù)量變少，因此導(dǎo)致方差增加，a越大導(dǎo)致的結(jié)果相反，也即是存在一個偏倚和方差的平衡選擇;此外，ZA條件也對混合分布模型的可識別性起了重要作用。作者選擇了a=0.5這個相對保險的截點(diǎn)，此時的數(shù)據(jù)量也足以精確估計參數(shù)，但是關(guān)于截點(diǎn)a的自適應(yīng)選擇仍是有必要的。

［1］Benjamini Y，Hochberg Y.Controlling the false discovery rate:a practical and powerful approach to multiple testing［J］.J Royal Statist Soc:Series B，1995，57(1):289

［2］Benjamini Y.Discovering the false discovery rate［J］.J Royal Statist Soc:Series B，2010，72(4):405

［3］Dudoit S，Shaffer JP，Boldrick JC.Multiple hypothesis testing in microarray experiments［J］.Statist Sci，2003，18(1):71

［4］Lazar N.The statistical analysis of functional MRI data［M］.New York:Springer，2008.

［5］Storey JD.The positive false discovery rate:a Bayesian interpretation and the q-value［J］.Ann Statist，2003，31 (6):2013

［6］Efron B，Tibshirani R，Storey JD，et al.Empirical Bayes analysis of a microarray experiment［J］.J Am Statist Ass，2001，96(456):1151

［7］Singh D，F(xiàn)ebbo PG，Ross K，et al.Gene expression correlates of clinical prostate cancer behavior［J］.Cancer Cell，2002，1(2):203

［8］Golub TR，Slonim DK，Tamayo P，et al.Molecular classification of cancer:class discovery and class prediction by gene expression monitoring［J］.Science，1999，286(5439):531

［9］Efron B.Microarrays，empirical Bayes，and the twogroups model［J］.Statistical Science，2008，23(11):1

［10］Efron B.Large-scale stimultaneous hypothesis testing:the choice of a hull hypothesis［J］.J Am Statist Ass，2004，99(1):96

［11］Efron B.Doing thousands of hypothesis tests at the same time［J］.Metron Int J Statist，2007，65(1):3

［12］Strimmer K.A unified approach to false discovery rate estimation［J］.BMC Bioinformatics，2008，9:303

Empirical estimation and application of false discovery rate

WANG Ting1)，ZENG Ping1，2)，HUANG Shuiping1)，ZHAO Huashuo1)1)Department of Epidemiology and Health Statistics，School of Public Health，Xuzhou Medical College，Xuzhou 221002 2)Department of Epidemiology and Health Statistics，School of Public Health，Nanjing Medical University，Nanjing 210029

microarray data;false discovery rate;empirical Bayes;density estimation

Aim:To investigate the empirical Bayesian estimation and application for density，null distribution and false discovery rate in large scale data.Methods:A Bayesian two-group model was constructed for two microarray data，density function was estimated using the method of Poisson regression，and then the empirical Bayes was applied to estimate false discovery rate and local false discovery rate based on the estimated density.Results:The method of Poisson regression for density estimation was an appropriate alternative for the empirical Bayesian estimation for null distribution and false discovery rate.Conclusion:It is feasible to estimate false discovery rate and null distribution empirically by employing the parallel data structure in large scale data.

R195.1

10.3969/j.issn.1671-6825.2012.05.014

*江蘇省教育廳高校哲學(xué)社會科學(xué)研究基金資助項目2010SJB790037;徐州醫(yī)學(xué)院公共衛(wèi)生學(xué)院科研課題資助項目201107，201115

(2011-09-07收稿 責(zé)任編輯李沛寰)