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

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

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

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

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

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

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

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

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

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

1 方法與原理

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

Φ、Tv分別表示正態(tài)和自由度為v的t變量累計(jì)分布函數(shù)。圖1給出了2個(gè)微陣列數(shù)據(jù)z值的直方圖［7-8］，前者描述了50個(gè)正常對(duì)照和52個(gè)前列腺癌患者6 033個(gè)基因的表達(dá)水平，后者描述了45例急性淋巴細(xì)胞白血病患者和27例急性髓性白血病患者7 128個(gè)基因的表達(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的先驗(yàn)概率為π0，H1=1的先驗(yàn)概率為1－π0，H0=0時(shí)z值的密度和分布函數(shù)分別為f0、F0，H1=1時(shí)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)可看作是后驗(yàn)尾部面積，locfdr(z)可看作是后驗(yàn)概率，因此從貝葉斯角度看locfdr(z)更具有可解釋性。兩者之間的關(guān)系為E(locfdr)= FDR。

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

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

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

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

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

2 結(jié)果

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

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

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

3 討論

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

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

高維數(shù)據(jù)分析中另一個(gè)十分重要的問題是，理論無效分布常常被違背，如白血病數(shù)據(jù)，在這種情況FDR的估計(jì)明顯錯(cuò)誤，但大規(guī)模平行數(shù)據(jù)結(jié)構(gòu)能夠?qū)o效分布進(jìn)行經(jīng)驗(yàn)估計(jì)。作者選擇ZA條件下的匹配估計(jì)方法發(fā)現(xiàn)，無效分布的參數(shù)依賴選擇的固定常數(shù)a，但是在一定范圍內(nèi)a的影響有限。a越小，對(duì)ZA條件的信心越高，此時(shí)偏倚減小，但同時(shí)用于估計(jì)參數(shù)的數(shù)據(jù)量變少，因此導(dǎo)致方差增加，a越大導(dǎo)致的結(jié)果相反，也即是存在一個(gè)偏倚和方差的平衡選擇;此外，ZA條件也對(duì)混合分布模型的可識(shí)別性起了重要作用。作者選擇了a=0.5這個(gè)相對(duì)保險(xiǎn)的截點(diǎn)，此時(shí)的數(shù)據(jù)量也足以精確估計(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é)社會(huì)科學(xué)研究基金資助項(xiàng)目2010SJB790037;徐州醫(yī)學(xué)院公共衛(wèi)生學(xué)院科研課題資助項(xiàng)目201107，201115

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