基于non-local先驗(yàn)的貝葉斯變量選擇方法及其在高維數(shù)據(jù)分析中的應(yīng)用*

山西醫(yī)科大學(xué)公共衛(wèi)生學(xué)院衛(wèi)生統(tǒng)計(jì)教研室(030001)

馬金沙 董曉強(qiáng) 高 倩 陶 然 許樹紅 李艷艷 王 彤△

【提 要】 目的 對(duì)高維數(shù)據(jù)進(jìn)行變量篩選并構(gòu)建預(yù)測(cè)模型是組學(xué)數(shù)據(jù)分析的研究熱點(diǎn)之一。本研究旨在為結(jié)局為二分類變量的高維組學(xué)數(shù)據(jù)篩選自變量并構(gòu)建預(yù)測(cè)結(jié)局的稀疏統(tǒng)計(jì)模型。方法 本研究通過模擬研究和實(shí)例分析闡釋基于non-local先驗(yàn)的貝葉斯變量選擇方法——乘積逆矩先驗(yàn)(product inverse moment，piMOM)相較于懲罰類方法ISIS-光滑平切絕對(duì)偏差(iterative sure independence screening-smoothly clipped absolute deviation,ISIS-SCAD)和ISIS-最小最大凹懲罰(iterative sure independence screening-minimax concave penalty,ISIS-MCP)在高維數(shù)據(jù)中變量篩選及其預(yù)測(cè)效果的性能優(yōu)劣。結(jié)果 模擬研究發(fā)現(xiàn)：在高維的情況下，經(jīng)piMOM、ISIS-SCAD和ISIS-MCP方法篩選所得變量的平均真陽性數(shù)和受試者工作特征曲線下面積(AUC，area under curve)基本相等，ISIS-SCAD、ISIS-MCP的平均假陽性數(shù)、回歸系數(shù)均方誤差以及預(yù)測(cè)均方誤差明顯高于基于non-local先驗(yàn)的貝葉斯變量方法所獲得的對(duì)應(yīng)值。piMOM方法分析彌漫大B細(xì)胞淋巴瘤實(shí)例數(shù)據(jù)共識(shí)別5個(gè)有意義的基因，AUC為0.996；ISIS-SCAD識(shí)別7個(gè)基因，AUC為0.975；ISIS-MCP識(shí)別7個(gè)基因，AUC為0.968。結(jié)論 在模型選擇相合性和預(yù)測(cè)準(zhǔn)確性方面，piMOM方法與ISIS-SCAD和ISIS-MCP相比，具有優(yōu)勢(shì)，在一定意義上可有效控制假陽性率。

測(cè)序技術(shù)的發(fā)展催生大量組學(xué)數(shù)據(jù)，此類數(shù)據(jù)中，變量個(gè)數(shù)p遠(yuǎn)大于觀測(cè)個(gè)體數(shù)n(p=O(n))且變量間往往呈現(xiàn)高度相關(guān)性，稱為高維數(shù)據(jù)[1]。目前，解決高維數(shù)據(jù)變量篩選問題有三種基本思路：最優(yōu)子集選擇、懲罰以及貝葉斯方法。最優(yōu)子集選擇方法在低維數(shù)據(jù)中表現(xiàn)良好，但在高維數(shù)據(jù)中，估計(jì)結(jié)果不穩(wěn)定，加之其較大的計(jì)算量，不再適用。

懲罰類方法以懲罰項(xiàng)為依據(jù)，分為光滑平切絕對(duì)偏差[2](smoothly clipped absolute deviation,SCAD)、自適應(yīng)Lasso[3](adaptive lasso,aLasso)、最小最大凹懲罰[4](minimax concave penalty,MCP)等多種類型。在高維數(shù)據(jù)中，因計(jì)算量較大使懲罰類方法的運(yùn)用受到限制。基于該問題，F(xiàn)an和Lv[1]于2008年提出具有sure screening 性質(zhì)的(I)SIS(sure independence screening)方法，并通過模擬研究闡明了其應(yīng)用于高維數(shù)據(jù)分析時(shí)的良好效果。2016年，王彤團(tuán)隊(duì)[5]將ISIS+Lasso、ISIS+SCAD、ISIS+MCP三種方法分別應(yīng)用到模擬研究和實(shí)例分析中，探討了其應(yīng)用于高維數(shù)據(jù)二分類結(jié)局資料中實(shí)現(xiàn)變量選擇目的時(shí)的表現(xiàn)。

貝葉斯方法使用馬爾可夫鏈蒙特卡洛算法(Markov chain Monte Carlo algorithm,MCMC)，從后驗(yàn)分布中抽取待估計(jì)的參數(shù)值。2010年，基于貝葉斯假設(shè)檢驗(yàn)理論，Johnson和Rossell等人[6]將貝葉斯先驗(yàn)分為兩類：local 與non-local 先驗(yàn)。2010年，Johnson和Rossel等[6]引入乘積逆矩(product inverse moment,piMOM)和乘積矩(product moment,pMOM)。2013年Rossell等[7]提出乘積指數(shù)矩(product exponential moment,peMOM)，同時(shí)探究各貝葉斯方法用于高維數(shù)據(jù)變量選擇的性能。上述三種先驗(yàn)均為non-local先驗(yàn)，見圖1。該圖提示在高維組學(xué)數(shù)據(jù)變量選擇方面，基于non-local先驗(yàn)的貝葉斯方法優(yōu)于local先驗(yàn)。

本文首先闡釋基于non-local先驗(yàn)的貝葉斯變量篩選方法——piMOM，繼而通過模擬研究和實(shí)例分析探究其在高維數(shù)據(jù)中的應(yīng)用，從而比較piMOM、ISIS-SCAD、ISIS-MCP方法在高維數(shù)據(jù)中進(jìn)行變量篩選的性能優(yōu)劣。

圖1 備擇假設(shè)下模型參數(shù)先驗(yàn)分布的比較左：local先驗(yàn) 右：non-local先驗(yàn)

原理與方法

1.基于ISIS的懲罰類方法

該方法實(shí)現(xiàn)分為兩個(gè)階段：第一階段，從高維數(shù)據(jù)的多個(gè)自變量中識(shí)別出與結(jié)局變量相關(guān)性相對(duì)較強(qiáng)的變量。第二階段，基于懲罰類思想，從與結(jié)局變量強(qiáng)相關(guān)的變量中篩選出真正具有意義的自變量，本研究選用SCAD、MCP兩種懲罰類方法。其中SCAD[2]的懲罰項(xiàng)為：

(1)

MCP[4]的懲罰項(xiàng)為：

(2)

α，λ為調(diào)整參數(shù)。

2.基于non-local先驗(yàn)貝葉斯變量選擇方法：piMOM

應(yīng)用該方法進(jìn)行變量篩選時(shí)，需對(duì)模型參數(shù)先驗(yàn)分布、包括模型空間先驗(yàn)進(jìn)行合理設(shè)定，并使用恰當(dāng)方法以先驗(yàn)和似然獲得的模型為依據(jù)估計(jì)后驗(yàn)概率，這兩個(gè)過程是保證貝葉斯方法能夠合理進(jìn)行變量選擇的關(guān)鍵。

(1)模型參數(shù)先驗(yàn)分布及其超參數(shù)設(shè)置

本文采用的non-local先驗(yàn)為乘積逆矩先驗(yàn)，表達(dá)式為：

(3)

其中，τ,γ>0為piMOM的兩個(gè)超參數(shù)。τ為尺度參數(shù)，γ為形狀參數(shù)。上述兩個(gè)超參數(shù)分別決定先驗(yàn)函數(shù)0附近和兩端尾部的分布情況。某種意義上，所構(gòu)建模型中參數(shù)的最小值由尺度參數(shù)τ決定。針對(duì)“如何對(duì)τ值進(jìn)行合理選擇”這一問題，Nikooienejad[8]于2016年給出相關(guān)建議：數(shù)據(jù)經(jīng)標(biāo)準(zhǔn)化后，能使原假設(shè)下和備擇假設(shè)下概率密度函數(shù)交叉面積低于一定閾值(p-α)的最大τ值，即為合理τ值。合理選取該值能在有效控制模型的假陽性率(兩者密度函數(shù)交叉部分)的同時(shí)，保證模型具有較高的靈敏度，見圖2。

圖2 基于non-local先驗(yàn)貝葉斯變量選擇方法中超參數(shù)選擇的模擬圖

(2)模型空間先驗(yàn)

模型k的beta-binomial先驗(yàn)[9]表達(dá)式為：

(4)

(2)模型的后驗(yàn)概率

基于貝葉斯理論知識(shí)，后驗(yàn)概率公式可記為：

(5)

通過(5)式可獲得最高后驗(yàn)概率模型(highest posterior probability model,HPPM)。鑒于參數(shù)空間的高維屬性和后驗(yàn)概率密度函數(shù)的復(fù)雜屬性，本研究運(yùn)用MCMC算法實(shí)現(xiàn)HPPM模型的最大化。同時(shí)，對(duì)各邊際概率密度函數(shù)進(jìn)行普拉斯(Laplace)近似處理以達(dá)到降維目的。其中，對(duì)應(yīng)于模型k的邊際概率密度函數(shù)如下所示：

(6)

(7)

最大后驗(yàn)估計(jì)值可借助函數(shù)nlminb()完成計(jì)算。通過上述所示赫氏矩陣表達(dá)式可求得目標(biāo)函數(shù)f(yn,βk)的局部極小值，在此基礎(chǔ)上利用該局部極小值得到邊際概率密度函數(shù)最大值(mk(yn))，最終獲得后驗(yàn)概率。Johnson等人[10-11]于1998年改進(jìn)收斂性評(píng)價(jià)方法，提出校正成對(duì)診斷(modified coupling diagnostic)，本文即選用該方法對(duì)馬爾可夫鏈蒙特卡洛算法的收斂性進(jìn)行評(píng)價(jià)。

3.軟件實(shí)現(xiàn)

本文運(yùn)用R 3.3.2軟件進(jìn)行統(tǒng)計(jì)分析。采用Fan[12]2015年提供的SIS包實(shí)現(xiàn)兩類懲罰類方法ISIS-SCAD和ISIS-MCP，并采用十折交叉驗(yàn)證方法進(jìn)行調(diào)整參數(shù)選擇[13]，根據(jù)Fan等的建議[2]，對(duì)SCAD方法，令a=3.7，對(duì)MCP方法，令a=3。采用2012年Johnson等人[14]提供的mombf包以實(shí)現(xiàn)基于non-local先驗(yàn)(piMOM)的貝葉斯變量篩選方法。同時(shí)調(diào)用snowfall并行運(yùn)算包以達(dá)到提高貝葉斯方法運(yùn)算效率的目的。

模擬試驗(yàn)

1.模擬研究設(shè)計(jì)

本研究樣本量設(shè)置5個(gè)水平，分別為50，100，200，400和600。自變量個(gè)數(shù)設(shè)置兩個(gè)水平，分別為1000和3000，該設(shè)置滿足Li和Fan提出的高維數(shù)據(jù)條件[1]。以Minsuk Shin[15](2015)模擬情境設(shè)置作為參考，假定模型k包含5個(gè)有效(即真正具有意義)自變量，并設(shè)自變量x服從均值0，協(xié)方差為∑的多元正態(tài)分布，即xi～Np(0,∑)，其中∑可分為如下不同情況：

方案一：令各自變量間互不相關(guān)，即∑jj′=0，j≠j′且∑jj=1，1≤j，j′≤p。

方案二：令各自變量間呈復(fù)合對(duì)稱相關(guān)，即∑jj′=0.5，j≠j′且∑jj=1，1≤j，j′≤p。

方案三：令各自變量間呈自回歸相關(guān)，即∑jj′=0.5|j-j′|，j≠j′且∑jj=1，1≤j，j′≤p。

2.模型評(píng)價(jià)指標(biāo)

所構(gòu)建模型優(yōu)劣的評(píng)價(jià)從兩個(gè)方面進(jìn)行：模型選擇相合性和模型預(yù)測(cè)準(zhǔn)確性。

(1)模型選擇相合性

模型選擇相合性的定義為三種方法所構(gòu)建模型與真實(shí)模型相一致的程度。可使用假陽性數(shù)(false positive,FP)和真陽性數(shù)(true positive,TP)進(jìn)行評(píng)價(jià)。本研究中平均真陽性數(shù)越接近模擬情境中設(shè)置的真正有意義的變量數(shù)5，平均假陽性數(shù)越接近0，則模型選擇相合度越高。

(2)模型預(yù)測(cè)準(zhǔn)確性

模型預(yù)測(cè)準(zhǔn)確性定義為三種方法所構(gòu)建模型預(yù)測(cè)的結(jié)局與真實(shí)結(jié)局的一致程度?？墒褂没貧w系數(shù)均方誤差(mean squared error of regression coefficient,RMSE)、預(yù)測(cè)均方誤差(mean squared error of prediction,PMSE)、以及受試者工作特征曲線下面積(area under curve，AUC)進(jìn)行評(píng)價(jià)。PMSE和RMSE越小，則所構(gòu)建模型越接近真實(shí)模型，預(yù)測(cè)越準(zhǔn)確。AUC值越大，提示模型分辨能力越高，預(yù)測(cè)越準(zhǔn)確。

3.模擬結(jié)果

(1)超參數(shù)選擇的模擬結(jié)果

piMOM方法在方案一中的超參數(shù)選擇結(jié)果如表1，方案二和方案三超參數(shù)選擇結(jié)果表2，表3。

表1 方案一piMOM方法的超參數(shù)選擇

表2 方案二piMOM方法的超參數(shù)選擇

表3 方案三piMOM方法的超參數(shù)選擇

由以上超參數(shù)選擇結(jié)果可知：超參數(shù)τ和γ隨著數(shù)據(jù)維數(shù)的增高而增大；反之，即數(shù)據(jù)維數(shù)降低時(shí)，超參數(shù)τ和γ隨之逐漸減小。

(2)自變量選擇的模擬結(jié)果

方案一模擬結(jié)果見圖3(方案二和方案三見圖4，5)。

對(duì)于piMOM方法，協(xié)方差矩陣結(jié)構(gòu)和樣本量一定時(shí)，模型TP、FP、PMSE、RMSE、AUC以及所獲得模型總變量數(shù)隨自變量維數(shù)變化較小。協(xié)方差結(jié)構(gòu)和自變量維數(shù)一定時(shí)，樣本量增加，TP增加明顯；FP始終處于較低水平，相對(duì)穩(wěn)定；PMSE在樣本量達(dá)200之前下降較快，其后相對(duì)穩(wěn)定；RMSE相對(duì)穩(wěn)定，始終處于較低水平；AUC相對(duì)穩(wěn)定，均>0.8；所獲模型總變量數(shù)隨樣本量增加逐漸增加，直至獲得所有的真陽性變量。自變量維數(shù)和樣本量一定時(shí)，協(xié)方差結(jié)構(gòu)改變，預(yù)測(cè)模型性能指標(biāo)變化較小。

圖3 方案一的模擬結(jié)果

圖4 方案二的模擬結(jié)果

圖5 方案三的模擬結(jié)果

對(duì)懲罰類方法ISIS-SCAD和ISIS-MCP，協(xié)方差矩陣結(jié)構(gòu)和樣本量一定時(shí)，模型TP、PMSE、RMSE、AUC以及所獲得模型總變量數(shù)隨自變量維數(shù)變化較小。協(xié)方差結(jié)構(gòu)和自變量維數(shù)一定時(shí)，樣本量增加，TP明顯增加；FP呈現(xiàn)先上升后下降的趨勢(shì)，且均在篩得所有真陽性變量前取得最大值；PMSE和RMSE在樣本量達(dá)100之前下降較快，其后相對(duì)穩(wěn)定；AUC變化不大，始終>0.8；所獲得模型總變量數(shù)逐漸增加，直至獲得所有的真陽性變量。自變量維數(shù)和樣本量一定時(shí)，協(xié)方差結(jié)構(gòu)改變，預(yù)測(cè)模型效性能指標(biāo)變化較小。

給定某一特定高維情況，在模型選擇相合性方面，piMOM方法所獲得模型的TP與ISIS-SCAD、ISIS-MCP大致相等，F(xiàn)P小于兩類懲罰類方法。在模型選擇準(zhǔn)確性方面，piMOM方法所獲得模型的PMSE、RMSE以及AUC小于兩類懲罰類方法。在模型選擇穩(wěn)定性方面，piMOM方法各評(píng)價(jià)指標(biāo)相較于兩類懲罰類方法變化幅度較小。

實(shí)例分析

1.數(shù)據(jù)來源

本研究所使用數(shù)據(jù)集可在高通量基因表達(dá)數(shù)據(jù)庫(gene expression omnibus,GEO)中獲得，GEO數(shù)據(jù)庫鏈接為https：//www.ncbi.nlm.nih.gov/geo/，數(shù)據(jù)集相應(yīng)編號(hào)為GSE10846。該數(shù)據(jù)集包含性別、治療方案、國際預(yù)后指數(shù)(international prognostic index,IPI)、中位生存期以及基因表達(dá)信息，共420個(gè)樣本，54675個(gè)基因。因現(xiàn)有認(rèn)知中，DLBCL共分為3型：生發(fā)中心B細(xì)胞樣型(germinal center B cell-like,GCB)，預(yù)后較好；表達(dá)活化的外周血B細(xì)胞和漿細(xì)胞特征基因的活化外周血B細(xì)胞樣型(activated B cell-like,ABC)，預(yù)后較差；無明確特征的異源性類型，即第3型，預(yù)后較差[16-17]。因第3型不明確，為排除性定義，且該型樣本量?jī)H為總樣本量的1/6(70/420)，故本研究考慮將第3型排除，選擇DLBCL的兩個(gè)分子學(xué)亞型：ABC(167)和GCB(183)進(jìn)行探究。

2.研究方法

首先使用R軟件包Bioconductor中的limma包對(duì)原始數(shù)據(jù)進(jìn)行經(jīng)驗(yàn)貝葉斯差異表達(dá)分析[18-19]，繼而借助四分位數(shù)間距[20]對(duì)數(shù)據(jù)進(jìn)行降維處理，從而獲得包括350個(gè)樣本(ABC：167；GCB：183)，3237個(gè)基因進(jìn)行最終變量選擇的樣本集。本研究采用十折交叉驗(yàn)證方法識(shí)別有效基因，以保證結(jié)果穩(wěn)定性。經(jīng)10次篩選后，將頻數(shù)大于5的變量作為有效基因選入最終模型。

3.結(jié)果

piMOM方法篩選結(jié)果顯示：共有5個(gè)有意義的基因(TNFRSF13B，BATF，LMO2，SYTL4，PALD1)，AUC為0.996；ISIS-SCAD結(jié)果顯示：共有7個(gè)有意義的基因(MYBL1，CYB5R2，MAML3，S1PR2，ASB13，ZBTB46，BATF)，AUC為0.975；ISIS-MCP結(jié)果顯示：共有7個(gè)有意義的基因(CYB5R2，MAML3，S1PR2，ASB13，ACVR2A，TNFRSF13B，BATF)，AUC為0.968。詳見表4

表4 三種方法基因選擇結(jié)果比較

*：“_”為已有文獻(xiàn)報(bào)道與DLBCL分型有關(guān)的基因。

討 論

本文將piMOM方法應(yīng)用到高維數(shù)據(jù)中構(gòu)建二分類logistic回歸模型并實(shí)現(xiàn)變量篩選，繼而以模擬研究和實(shí)例分析兩種途徑探討其相對(duì)于兩種懲罰類方法ISIS-SCAD、ISIS-MCP的優(yōu)劣。

模擬研究表明：相比于ISIS-SCAD和ISIS-MCP，piMOM方法在模型選擇相合性方面能準(zhǔn)確篩選有效自變量，具有較小的假陽性率；模型預(yù)測(cè)準(zhǔn)確性方面其PMSE和RMSE均較小。這種結(jié)果是由貝葉斯方法和傳統(tǒng)懲罰類方法在較小非0模型參數(shù)附近不同的概率密度函數(shù)所導(dǎo)致：從貝葉斯角度看，兩類懲罰類方法SCAD、MCP方法傾向于使得較小非0模型參數(shù)附近先驗(yàn)概率更大，相當(dāng)于在模型參數(shù)上施加local先驗(yàn)；而piMOM方法則不同，其通過BIC模型賦予較小非0的模型參數(shù)以較小的先驗(yàn)概率，原理如下：BIC模型選擇的目標(biāo)函數(shù)為：

(8)

(9)

d>0，通過對(duì)(8)、(9)式的比較可知：以BIC為基礎(chǔ)，增加一個(gè)懲罰項(xiàng)后即可獲得piMOM目標(biāo)函數(shù)，這使得piMOM方法給予0附近(包含0)模型參數(shù)以無窮大懲罰，從而保證所建立終模型的假陽性率較小。

在模擬研究中，經(jīng)ISIS-SCAD和ISIS-MCP方法篩選所獲得模型，變量的假陽性數(shù)呈現(xiàn)先上升后下降的趨勢(shì)，推測(cè)原因如下：在自變量個(gè)數(shù)p不變的情況下，隨著樣本量的增加，數(shù)據(jù)提供的信息會(huì)越多，所能篩得的總的自變量個(gè)數(shù)相對(duì)增加，起初真陽性數(shù)和假陽性數(shù)都隨之增加，至獲得所有真陽性自變量前，模型中所獲得的真陽性變量漸能代表樣本信息，假陽性數(shù)繼而出現(xiàn)下降趨勢(shì)，直至獲得所有真陽性數(shù)時(shí)，假陽性數(shù)降至最低。而對(duì)于貝葉斯而言，則可以通過先驗(yàn)分布中超參數(shù)值的選擇控制假陽性數(shù)始終處于較低水平。

實(shí)例研究表明：基于non-local先驗(yàn)的貝葉斯變量篩選方法所構(gòu)建的模型能夠以較少的變量數(shù)，獲得與ISIS-SCAD、ISIS-MCP近乎相同的AUC。若將貝葉斯變量篩選的結(jié)果應(yīng)用于實(shí)踐，可縮小后期進(jìn)行實(shí)驗(yàn)驗(yàn)證的探索范圍，降低經(jīng)濟(jì)成本，為更好地進(jìn)行疾病機(jī)制研究提供便利。從這一角度上講，基于non-local先驗(yàn)的貝葉斯方法相較于懲罰類方法更可取。

本探究尚有不足之處：(1)貝葉斯變量篩選方法運(yùn)算量較大，將其應(yīng)用于高維組學(xué)數(shù)據(jù)時(shí)，該問題帶來的計(jì)算困難尤為突出，如何提高貝葉斯方法的計(jì)算效率，為該領(lǐng)域亟待解決的問題；(2)本文未將基于local先驗(yàn)的貝葉斯變量選擇方法與non-local先驗(yàn)同時(shí)應(yīng)用到二分類結(jié)局變量資料的高維數(shù)據(jù)中，以完成兩種先驗(yàn)之間的對(duì)比，關(guān)于貝葉斯方法變量選擇性能尚需作進(jìn)一步探究。