山西醫(yī)科大學公共衛(wèi)生學院衛(wèi)生統(tǒng)計教研室(030001)
馬金沙 董曉強 高 倩 陶 然 許樹紅 李艷艷 王 彤△
【提 要】 目的 對高維數(shù)據(jù)進行變量篩選并構建預測模型是組學數(shù)據(jù)分析的研究熱點之一。本研究旨在為結局為二分類變量的高維組學數(shù)據(jù)篩選自變量并構建預測結局的稀疏統(tǒng)計模型。方法 本研究通過模擬研究和實例分析闡釋基于non-local先驗的貝葉斯變量選擇方法——乘積逆矩先驗(product inverse moment,piMOM)相較于懲罰類方法ISIS-光滑平切絕對偏差(iterative sure independence screening-smoothly clipped absolute deviation,ISIS-SCAD)和ISIS-最小最大凹懲罰(iterative sure independence screening-minimax concave penalty,ISIS-MCP)在高維數(shù)據(jù)中變量篩選及其預測效果的性能優(yōu)劣。結果 模擬研究發(fā)現(xiàn):在高維的情況下,經piMOM、ISIS-SCAD和ISIS-MCP方法篩選所得變量的平均真陽性數(shù)和受試者工作特征曲線下面積(AUC,area under curve)基本相等,ISIS-SCAD、ISIS-MCP的平均假陽性數(shù)、回歸系數(shù)均方誤差以及預測均方誤差明顯高于基于non-local先驗的貝葉斯變量方法所獲得的對應值。piMOM方法分析彌漫大B細胞淋巴瘤實例數(shù)據(jù)共識別5個有意義的基因,AUC為0.996;ISIS-SCAD識別7個基因,AUC為0.975;ISIS-MCP識別7個基因,AUC為0.968。結論 在模型選擇相合性和預測準確性方面,piMOM方法與ISIS-SCAD和ISIS-MCP相比,具有優(yōu)勢,在一定意義上可有效控制假陽性率。
測序技術的發(fā)展催生大量組學數(shù)據(jù),此類數(shù)據(jù)中,變量個數(shù)p遠大于觀測個體數(shù)n(p=O(n))且變量間往往呈現(xiàn)高度相關性,稱為高維數(shù)據(jù)[1]。目前,解決高維數(shù)據(jù)變量篩選問題有三種基本思路:最優(yōu)子集選擇、懲罰以及貝葉斯方法。最優(yōu)子集選擇方法在低維數(shù)據(jù)中表現(xiàn)良好,但在高維數(shù)據(jù)中,估計結果不穩(wěn)定,加之其較大的計算量,不再適用。
懲罰類方法以懲罰項為依據(jù),分為光滑平切絕對偏差[2](smoothly clipped absolute deviation,SCAD)、自適應Lasso[3](adaptive lasso,aLasso)、最小最大凹懲罰[4](minimax concave penalty,MCP)等多種類型。在高維數(shù)據(jù)中,因計算量較大使懲罰類方法的運用受到限制?;谠搯栴},F(xiàn)an和Lv[1]于2008年提出具有sure screening 性質的(I)SIS(sure independence screening)方法,并通過模擬研究闡明了其應用于高維數(shù)據(jù)分析時的良好效果。2016年,王彤團隊[5]將ISIS+Lasso、ISIS+SCAD、ISIS+MCP三種方法分別應用到模擬研究和實例分析中,探討了其應用于高維數(shù)據(jù)二分類結局資料中實現(xiàn)變量選擇目的時的表現(xiàn)。
貝葉斯方法使用馬爾可夫鏈蒙特卡洛算法(Markov chain Monte Carlo algorithm,MCMC),從后驗分布中抽取待估計的參數(shù)值。2010年,基于貝葉斯假設檢驗理論,Johnson和Rossell等人[6]將貝葉斯先驗分為兩類:local 與non-local 先驗。2010年,Johnson和Rossel等[6]引入乘積逆矩(product inverse moment,piMOM)和乘積矩(product moment,pMOM)。2013年Rossell等[7]提出乘積指數(shù)矩(product exponential moment,peMOM),同時探究各貝葉斯方法用于高維數(shù)據(jù)變量選擇的性能。上述三種先驗均為non-local先驗,見圖1。該圖提示在高維組學數(shù)據(jù)變量選擇方面,基于non-local先驗的貝葉斯方法優(yōu)于local先驗。
本文首先闡釋基于non-local先驗的貝葉斯變量篩選方法——piMOM,繼而通過模擬研究和實例分析探究其在高維數(shù)據(jù)中的應用,從而比較piMOM、ISIS-SCAD、ISIS-MCP方法在高維數(shù)據(jù)中進行變量篩選的性能優(yōu)劣。
圖1 備擇假設下模型參數(shù)先驗分布的比較左:local先驗 右:non-local先驗
1.基于ISIS的懲罰類方法
該方法實現(xiàn)分為兩個階段:第一階段,從高維數(shù)據(jù)的多個自變量中識別出與結局變量相關性相對較強的變量。第二階段,基于懲罰類思想,從與結局變量強相關的變量中篩選出真正具有意義的自變量,本研究選用SCAD、MCP兩種懲罰類方法。其中SCAD[2]的懲罰項為:
(1)
MCP[4]的懲罰項為:
(2)
α,λ為調整參數(shù)。
2.基于non-local先驗貝葉斯變量選擇方法:piMOM
應用該方法進行變量篩選時,需對模型參數(shù)先驗分布、包括模型空間先驗進行合理設定,并使用恰當方法以先驗和似然獲得的模型為依據(jù)估計后驗概率,這兩個過程是保證貝葉斯方法能夠合理進行變量選擇的關鍵。
(1)模型參數(shù)先驗分布及其超參數(shù)設置
本文采用的non-local先驗為乘積逆矩先驗,表達式為:
(3)
其中,τ,γ>0為piMOM的兩個超參數(shù)。τ為尺度參數(shù),γ為形狀參數(shù)。上述兩個超參數(shù)分別決定先驗函數(shù)0附近和兩端尾部的分布情況。某種意義上,所構建模型中參數(shù)的最小值由尺度參數(shù)τ決定。針對“如何對τ值進行合理選擇”這一問題,Nikooienejad[8]于2016年給出相關建議:數(shù)據(jù)經標準化后,能使原假設下和備擇假設下概率密度函數(shù)交叉面積低于一定閾值(p-α)的最大τ值,即為合理τ值。合理選取該值能在有效控制模型的假陽性率(兩者密度函數(shù)交叉部分)的同時,保證模型具有較高的靈敏度,見圖2。
圖2 基于non-local先驗貝葉斯變量選擇方法中超參數(shù)選擇的模擬圖
(2)模型空間先驗
模型k的beta-binomial先驗[9]表達式為:
(4)
(2)模型的后驗概率
基于貝葉斯理論知識,后驗概率公式可記為:
(5)
通過(5)式可獲得最高后驗概率模型(highest posterior probability model,HPPM)。鑒于參數(shù)空間的高維屬性和后驗概率密度函數(shù)的復雜屬性,本研究運用MCMC算法實現(xiàn)HPPM模型的最大化。同時,對各邊際概率密度函數(shù)進行普拉斯(Laplace)近似處理以達到降維目的。其中,對應于模型k的邊際概率密度函數(shù)如下所示:
(6)
(7)
最大后驗估計值可借助函數(shù)nlminb()完成計算。通過上述所示赫氏矩陣表達式可求得目標函數(shù)f(yn,βk)的局部極小值,在此基礎上利用該局部極小值得到邊際概率密度函數(shù)最大值(mk(yn)),最終獲得后驗概率。Johnson等人[10-11]于1998年改進收斂性評價方法,提出校正成對診斷(modified coupling diagnostic),本文即選用該方法對馬爾可夫鏈蒙特卡洛算法的收斂性進行評價。
3.軟件實現(xiàn)
本文運用R 3.3.2軟件進行統(tǒng)計分析。采用Fan[12]2015年提供的SIS包實現(xiàn)兩類懲罰類方法ISIS-SCAD和ISIS-MCP,并采用十折交叉驗證方法進行調整參數(shù)選擇[13],根據(jù)Fan等的建議[2],對SCAD方法,令a=3.7,對MCP方法,令a=3。采用2012年Johnson等人[14]提供的mombf包以實現(xiàn)基于non-local先驗(piMOM)的貝葉斯變量篩選方法。同時調用snowfall并行運算包以達到提高貝葉斯方法運算效率的目的。
1.模擬研究設計
本研究樣本量設置5個水平,分別為50,100,200,400和600。自變量個數(shù)設置兩個水平,分別為1000和3000,該設置滿足Li和Fan提出的高維數(shù)據(jù)條件[1]。以Minsuk Shin[15](2015)模擬情境設置作為參考,假定模型k包含5個有效(即真正具有意義)自變量,并設自變量x服從均值0,協(xié)方差為∑的多元正態(tài)分布,即xi~Np(0,∑),其中∑可分為如下不同情況:
方案一:令各自變量間互不相關,即∑jj′=0,j≠j′且∑jj=1,1≤j,j′≤p。
方案二:令各自變量間呈復合對稱相關,即∑jj′=0.5,j≠j′且∑jj=1,1≤j,j′≤p。
方案三:令各自變量間呈自回歸相關,即∑jj′=0.5|j-j′|,j≠j′且∑jj=1,1≤j,j′≤p。
2.模型評價指標
所構建模型優(yōu)劣的評價從兩個方面進行:模型選擇相合性和模型預測準確性。
(1)模型選擇相合性
模型選擇相合性的定義為三種方法所構建模型與真實模型相一致的程度。可使用假陽性數(shù)(false positive,FP)和真陽性數(shù)(true positive,TP)進行評價。本研究中平均真陽性數(shù)越接近模擬情境中設置的真正有意義的變量數(shù)5,平均假陽性數(shù)越接近0,則模型選擇相合度越高。
(2)模型預測準確性
模型預測準確性定義為三種方法所構建模型預測的結局與真實結局的一致程度??墒褂没貧w系數(shù)均方誤差(mean squared error of regression coefficient,RMSE)、預測均方誤差(mean squared error of prediction,PMSE)、以及受試者工作特征曲線下面積(area under curve,AUC)進行評價。PMSE和RMSE越小,則所構建模型越接近真實模型,預測越準確。AUC值越大,提示模型分辨能力越高,預測越準確。
3.模擬結果
(1)超參數(shù)選擇的模擬結果
piMOM方法在方案一中的超參數(shù)選擇結果如表1,方案二和方案三超參數(shù)選擇結果表2,表3。
表1 方案一piMOM方法的超參數(shù)選擇
表2 方案二piMOM方法的超參數(shù)選擇
表3 方案三piMOM方法的超參數(shù)選擇
由以上超參數(shù)選擇結果可知:超參數(shù)τ和γ隨著數(shù)據(jù)維數(shù)的增高而增大;反之,即數(shù)據(jù)維數(shù)降低時,超參數(shù)τ和γ隨之逐漸減小。
(2)自變量選擇的模擬結果
方案一模擬結果見圖3(方案二和方案三見圖4,5)。
對于piMOM方法,協(xié)方差矩陣結構和樣本量一定時,模型TP、FP、PMSE、RMSE、AUC以及所獲得模型總變量數(shù)隨自變量維數(shù)變化較小。協(xié)方差結構和自變量維數(shù)一定時,樣本量增加,TP增加明顯;FP始終處于較低水平,相對穩(wěn)定;PMSE在樣本量達200之前下降較快,其后相對穩(wěn)定;RMSE相對穩(wěn)定,始終處于較低水平;AUC相對穩(wěn)定,均>0.8;所獲模型總變量數(shù)隨樣本量增加逐漸增加,直至獲得所有的真陽性變量。自變量維數(shù)和樣本量一定時,協(xié)方差結構改變,預測模型性能指標變化較小。
圖3 方案一的模擬結果
圖4 方案二的模擬結果
圖5 方案三的模擬結果
對懲罰類方法ISIS-SCAD和ISIS-MCP,協(xié)方差矩陣結構和樣本量一定時,模型TP、PMSE、RMSE、AUC以及所獲得模型總變量數(shù)隨自變量維數(shù)變化較小。協(xié)方差結構和自變量維數(shù)一定時,樣本量增加,TP明顯增加;FP呈現(xiàn)先上升后下降的趨勢,且均在篩得所有真陽性變量前取得最大值;PMSE和RMSE在樣本量達100之前下降較快,其后相對穩(wěn)定;AUC變化不大,始終>0.8;所獲得模型總變量數(shù)逐漸增加,直至獲得所有的真陽性變量。自變量維數(shù)和樣本量一定時,協(xié)方差結構改變,預測模型效性能指標變化較小。
給定某一特定高維情況,在模型選擇相合性方面,piMOM方法所獲得模型的TP與ISIS-SCAD、ISIS-MCP大致相等,F(xiàn)P小于兩類懲罰類方法。在模型選擇準確性方面,piMOM方法所獲得模型的PMSE、RMSE以及AUC小于兩類懲罰類方法。在模型選擇穩(wěn)定性方面,piMOM方法各評價指標相較于兩類懲罰類方法變化幅度較小。
1.數(shù)據(jù)來源
本研究所使用數(shù)據(jù)集可在高通量基因表達數(shù)據(jù)庫(gene expression omnibus,GEO)中獲得,GEO數(shù)據(jù)庫鏈接為https://www.ncbi.nlm.nih.gov/geo/,數(shù)據(jù)集相應編號為GSE10846。該數(shù)據(jù)集包含性別、治療方案、國際預后指數(shù)(international prognostic index,IPI)、中位生存期以及基因表達信息,共420個樣本,54675個基因。因現(xiàn)有認知中,DLBCL共分為3型:生發(fā)中心B細胞樣型(germinal center B cell-like,GCB),預后較好;表達活化的外周血B細胞和漿細胞特征基因的活化外周血B細胞樣型(activated B cell-like,ABC),預后較差;無明確特征的異源性類型,即第3型,預后較差[16-17]。因第3型不明確,為排除性定義,且該型樣本量僅為總樣本量的1/6(70/420),故本研究考慮將第3型排除,選擇DLBCL的兩個分子學亞型:ABC(167)和GCB(183)進行探究。
2.研究方法
首先使用R軟件包Bioconductor中的limma包對原始數(shù)據(jù)進行經驗貝葉斯差異表達分析[18-19],繼而借助四分位數(shù)間距[20]對數(shù)據(jù)進行降維處理,從而獲得包括350個樣本(ABC:167;GCB:183),3237個基因進行最終變量選擇的樣本集。本研究采用十折交叉驗證方法識別有效基因,以保證結果穩(wěn)定性。經10次篩選后,將頻數(shù)大于5的變量作為有效基因選入最終模型。
3.結果
piMOM方法篩選結果顯示:共有5個有意義的基因(TNFRSF13B,BATF,LMO2,SYTL4,PALD1),AUC為0.996;ISIS-SCAD結果顯示:共有7個有意義的基因(MYBL1,CYB5R2,MAML3,S1PR2,ASB13,ZBTB46,BATF),AUC為0.975;ISIS-MCP結果顯示:共有7個有意義的基因(CYB5R2,MAML3,S1PR2,ASB13,ACVR2A,TNFRSF13B,BATF),AUC為0.968。詳見表4
表4 三種方法基因選擇結果比較
*:“_”為已有文獻報道與DLBCL分型有關的基因。
本文將piMOM方法應用到高維數(shù)據(jù)中構建二分類logistic回歸模型并實現(xiàn)變量篩選,繼而以模擬研究和實例分析兩種途徑探討其相對于兩種懲罰類方法ISIS-SCAD、ISIS-MCP的優(yōu)劣。
模擬研究表明:相比于ISIS-SCAD和ISIS-MCP,piMOM方法在模型選擇相合性方面能準確篩選有效自變量,具有較小的假陽性率;模型預測準確性方面其PMSE和RMSE均較小。這種結果是由貝葉斯方法和傳統(tǒng)懲罰類方法在較小非0模型參數(shù)附近不同的概率密度函數(shù)所導致:從貝葉斯角度看,兩類懲罰類方法SCAD、MCP方法傾向于使得較小非0模型參數(shù)附近先驗概率更大,相當于在模型參數(shù)上施加local先驗;而piMOM方法則不同,其通過BIC模型賦予較小非0的模型參數(shù)以較小的先驗概率,原理如下:BIC模型選擇的目標函數(shù)為:
(8)
(9)
d>0,通過對(8)、(9)式的比較可知:以BIC為基礎,增加一個懲罰項后即可獲得piMOM目標函數(shù),這使得piMOM方法給予0附近(包含0)模型參數(shù)以無窮大懲罰,從而保證所建立終模型的假陽性率較小。
在模擬研究中,經ISIS-SCAD和ISIS-MCP方法篩選所獲得模型,變量的假陽性數(shù)呈現(xiàn)先上升后下降的趨勢,推測原因如下:在自變量個數(shù)p不變的情況下,隨著樣本量的增加,數(shù)據(jù)提供的信息會越多,所能篩得的總的自變量個數(shù)相對增加,起初真陽性數(shù)和假陽性數(shù)都隨之增加,至獲得所有真陽性自變量前,模型中所獲得的真陽性變量漸能代表樣本信息,假陽性數(shù)繼而出現(xiàn)下降趨勢,直至獲得所有真陽性數(shù)時,假陽性數(shù)降至最低。而對于貝葉斯而言,則可以通過先驗分布中超參數(shù)值的選擇控制假陽性數(shù)始終處于較低水平。
實例研究表明:基于non-local先驗的貝葉斯變量篩選方法所構建的模型能夠以較少的變量數(shù),獲得與ISIS-SCAD、ISIS-MCP近乎相同的AUC。若將貝葉斯變量篩選的結果應用于實踐,可縮小后期進行實驗驗證的探索范圍,降低經濟成本,為更好地進行疾病機制研究提供便利。從這一角度上講,基于non-local先驗的貝葉斯方法相較于懲罰類方法更可取。
本探究尚有不足之處:(1)貝葉斯變量篩選方法運算量較大,將其應用于高維組學數(shù)據(jù)時,該問題帶來的計算困難尤為突出,如何提高貝葉斯方法的計算效率,為該領域亟待解決的問題;(2)本文未將基于local先驗的貝葉斯變量選擇方法與non-local先驗同時應用到二分類結局變量資料的高維數(shù)據(jù)中,以完成兩種先驗之間的對比,關于貝葉斯方法變量選擇性能尚需作進一步探究。