基于自適應(yīng)Lasso的缺失數(shù)據(jù)圖模型選擇

鄭倩貞， 徐平峰

(長(zhǎng)春工業(yè)大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院， 吉林 長(zhǎng)春 130012)

0 引 言

在高斯圖模型中，協(xié)方差逆陣可解釋隨機(jī)變量間的條件獨(dú)立性。矩陣中的零元素表示對(duì)應(yīng)兩個(gè)隨機(jī)變量間條件獨(dú)立，因此,協(xié)方差逆陣中零元素的估計(jì)問(wèn)題是圖模型選擇中的關(guān)鍵問(wèn)題。懲罰似然方法可同時(shí)實(shí)現(xiàn)模型選擇和參數(shù)估計(jì)，其中Lasso懲罰是常用的求解稀疏估計(jì)的懲罰項(xiàng)，例如：Yuan M等[1]、Banerjee O等[2]、Friedman J等[3]、Stidler N等[4]、徐平峰等[5]都采用l1懲罰對(duì)數(shù)似然方法估計(jì)協(xié)方差逆陣,從而得到稀疏的圖結(jié)構(gòu)。文獻(xiàn)[5]基于gCoda方法將成分?jǐn)?shù)據(jù)的圖模型選擇問(wèn)題轉(zhuǎn)換為高斯圖模型的協(xié)方差逆陣估計(jì)問(wèn)題；文獻(xiàn)[4]研究的是數(shù)據(jù)存在隨機(jī)缺失時(shí)的高斯圖模型選擇問(wèn)題；而其余三者關(guān)注的均是數(shù)據(jù)完全觀測(cè)下的高斯圖模型選擇問(wèn)題，但三者采用的優(yōu)化算法不同。

研究表明，Lasso懲罰可得到稀疏估計(jì)，但其估計(jì)結(jié)果會(huì)產(chǎn)生偏差[6]。Fan J等[6]提出的SCAD懲罰和Zou H[7]提出的自適應(yīng)Lasso懲罰都能夠降低估計(jì)的偏差。Fan J等[8]將這兩種懲罰應(yīng)用于完全數(shù)據(jù)的圖模型選擇問(wèn)題以減弱協(xié)方差逆陣估計(jì)的偏差。但在實(shí)際高維復(fù)雜數(shù)據(jù)中，往往存在數(shù)據(jù)部分缺失的情況。針對(duì)數(shù)據(jù)完全隨機(jī)缺失情形，文中將采用自適應(yīng)Lasso懲罰對(duì)數(shù)似然方法對(duì)高斯圖模型選擇進(jìn)行研究，并與Stidler N等[1]的l1懲罰對(duì)數(shù)似然方法(MissGLasso)進(jìn)行比較。

1 自適應(yīng)MissGLasso

設(shè)x=(xobs,xmis)為來(lái)自多元高斯分布Np(μ,Ω-1)的n個(gè)獨(dú)立同分布樣本。令樣本中的觀測(cè)數(shù)據(jù)xobs=(xobs,1,…,xobs,n)，缺失數(shù)據(jù)xmis=(xmis,1,…,xmis,n)，其中xobs,i和xmis,i分別表示第i個(gè)樣本的觀測(cè)數(shù)據(jù)和缺失數(shù)據(jù)的集合，i=1,2,…,n。與文獻(xiàn)[4]相同，文中假設(shè)數(shù)據(jù)的缺失機(jī)制為可忽略的?；谟^測(cè)數(shù)據(jù)xobs，可給出觀測(cè)對(duì)數(shù)似然函數(shù)

l(μ,Ω;xobs)=

((Ω-1)obs,i)-1×(xobs,i-μobs,i))，

(1)

式中：μobs,i----第i個(gè)樣本中觀測(cè)變量的均值;

(Ω-1)obs,i----第i個(gè)樣本中觀測(cè)變量的協(xié)方差陣。

在懲罰似然的框架下，下述優(yōu)化問(wèn)題的解可作為均值μ和協(xié)方差逆陣Ω的估計(jì)，

(2)

式中：λ----正則化參數(shù);

ωjk----Ω中的(j,k)元素;

p(·)----定義在各元素ωjk上的懲罰函數(shù)，j,k=1,2,…,p。

若令懲罰函數(shù)p(|ωjk|)=|ωjk|，則可得到Lasso懲罰，那么由上式求得的解即為MissGLasso估計(jì)[4]。

(3)

我們稱(chēng)該估計(jì)為自適應(yīng)MissGLasso(AdaMissGLasso)。針對(duì)上述優(yōu)化問(wèn)題，可結(jié)合EM算法[9]和GLasso算法[3]進(jìn)行優(yōu)化求解。

2 EM算法

EM算法常用于處理缺失數(shù)據(jù)下的參數(shù)估計(jì)問(wèn)題，該算法通過(guò)迭代交替進(jìn)行E(期望)步和M(最大化)步，直至算法收斂。E步是在給定觀測(cè)數(shù)據(jù)和當(dāng)前參數(shù)的情況下，計(jì)算完全數(shù)據(jù)的對(duì)數(shù)似然函數(shù)的條件期望；M步是關(guān)于待估參數(shù)極大化E步中的完全似然條件期望。令第t次迭代的估計(jì)值(μ(t),Ω(t))作為當(dāng)前參數(shù)，其中(μ(0),Ω(0))代表EM算法的初始值，設(shè)為均值插補(bǔ)后的GLasso估計(jì)。第(t+1)次迭代如下:

1)E步。

給定觀測(cè)數(shù)據(jù)xobs和當(dāng)前參數(shù)(μ(t),Ω(t))，計(jì)算完全數(shù)據(jù)懲罰對(duì)數(shù)似然函數(shù)的條件期望，記為Q函數(shù)。為簡(jiǎn)單起見(jiàn)，記E(·|xobs,μ(t),Ω(t))為Et(·)。于是，Q函數(shù)可表示為

Q(μ,Ω|μ(t),Ω(t))=

(4)

其中

與Fan J等[8]相同，文中取γ=0.5。

(5)

因此，E步只需計(jì)算Et(xij)和Et(xijxik)。

對(duì)第i個(gè)樣本xi=(xobs,i,xmis,i)，μ和Ω可分塊表示為:

給定觀測(cè)數(shù)據(jù)xobs,i和當(dāng)前參數(shù)(μ(t),Ω(t))，xmis,i的條件分布為N(ci,σi)，其中

于是，可推得[4,10-11]

(6)

以及

式中：Rij=1----xij為觀測(cè)數(shù)據(jù)；

Rij=0----xij缺失，i=1,2,…,n，j,k=1,2,…,p。

2)M步。

M步通過(guò)最小化Q(μ,Ω|μ(t),Ω(t))來(lái)更新參數(shù)估計(jì)(μ(t+1),Ω(t+1))，

(8)

式(8)中Ω(t+1)的優(yōu)化問(wèn)題可直接由GLasso算法求解。將(μ(t+1),Ω(t+1))作為下一次迭代的當(dāng)前參數(shù)，重復(fù)E步和M步直至算法收斂。

m=1,2,…,k，

(9)

其中

3 模擬研究

在完全隨機(jī)缺失的機(jī)制下，文中對(duì)自適應(yīng)MissGLasso方法和MissGLasso方法在高斯圖模型上的參數(shù)估計(jì)和模型選擇進(jìn)行了模擬比較?，F(xiàn)考慮以下兩個(gè)模型：

1)AR(1)，(Ω-1)jk=0.7|j-k|。

2)Ω=B+δI，其中B的對(duì)角線(xiàn)元素為0，非對(duì)角線(xiàn)元素之間相互獨(dú)立，取值為 0.5或0，且P(Bjk=0.5)=0.1。選擇恰當(dāng)?shù)摩?，使得Ω的條件數(shù)為p。

對(duì)于以上兩個(gè)模型，設(shè)定隨機(jī)變量個(gè)數(shù)p=100,200,300，對(duì)兩個(gè)模型與隨機(jī)變量個(gè)數(shù)p的6種組合，均產(chǎn)生樣本量n=200的50個(gè)獨(dú)立的數(shù)據(jù)集，且在每個(gè)數(shù)據(jù)集上完全隨機(jī)刪除10%、20%、30%的數(shù)據(jù)。因此，每個(gè)完全觀測(cè)的數(shù)據(jù)集對(duì)應(yīng)3個(gè)具有不同缺失比例的缺失數(shù)據(jù)集。為評(píng)價(jià)算法性能，文中比較以下指標(biāo)。

Kullback-Leibler損失(kl)、真陽(yáng)率(tpr)、陽(yáng)性預(yù)測(cè)率(ppv)以及馬修斯相關(guān)系數(shù)(mcc)。

式中:Ω----真實(shí)的協(xié)方差逆陣;

tp----真陽(yáng)類(lèi)的個(gè)數(shù);

tn----真陰類(lèi)的個(gè)數(shù);

fp----假陽(yáng)類(lèi)的個(gè)數(shù);

fn----假陰類(lèi)的個(gè)數(shù)。

基于模型1和模型2的不同參數(shù)設(shè)定各評(píng)價(jià)指標(biāo)的均值和標(biāo)準(zhǔn)差分別見(jiàn)表1和表2。

表1 不同參數(shù)設(shè)定下模型1各評(píng)價(jià)指標(biāo)的均值和(標(biāo)準(zhǔn)差)

表2 不同參數(shù)設(shè)定下模型2各評(píng)價(jià)指標(biāo)的均值和(標(biāo)準(zhǔn)差)

由表中可以看出，對(duì)于模型1和模型2的所有設(shè)定，除了tpr以外，自適應(yīng)MissGLasso方法均優(yōu)于MissGLasso，并且模型1的tpr僅在缺失比例為30%時(shí)略低。對(duì)于模型2，自適應(yīng)MissGLasso方法的mcc略低，但是兩者相差不大。kl為評(píng)價(jià)協(xié)方差逆陣估計(jì)效果的指標(biāo)，不論模型1還是模型2，自適應(yīng)MissGLasso的估計(jì)性能都更佳。綜上所述，自適應(yīng)的方法更優(yōu)。

4 實(shí)例分析

基于自適應(yīng)MissGLasso方法，文中對(duì)枯草芽孢桿菌生產(chǎn)核黃素(維生素B2)的數(shù)據(jù)集進(jìn)行了圖模型分析。該數(shù)據(jù)集由DSM(瑞士)提供，可在文獻(xiàn)[12]的補(bǔ)充材料中下載。數(shù)據(jù)為完全觀測(cè)，包含71個(gè)樣本，4 088個(gè)協(xié)變量以及1個(gè)響應(yīng)變量。協(xié)變量為4 088個(gè)基因表達(dá)水平的對(duì)數(shù)；響應(yīng)變量為核黃素生產(chǎn)率的對(duì)數(shù)。為簡(jiǎn)單展現(xiàn)自適應(yīng)MissGLasso方法的圖模型選擇效果，文中僅對(duì)經(jīng)驗(yàn)方差最大的100個(gè)基因，以及測(cè)量核黃素生產(chǎn)率的響應(yīng)變量作圖模型推斷，即數(shù)據(jù)集riboflavinv100.csv[12]。該數(shù)據(jù)集樣本量n=71，隨機(jī)變量個(gè)數(shù)p=101。

不同缺失程度下重疊邊數(shù)的箱線(xiàn)圖如圖1所示。

圖1 不同缺失比例下的重疊邊數(shù)

從圖1中可以看出，數(shù)據(jù)缺失的比例越高，重疊的邊數(shù)越少。即使數(shù)據(jù)集中缺失30%的數(shù)據(jù)，自適應(yīng)MissGLasso方法仍能識(shí)別出約30條重疊的邊。

不同缺失比例下的平均稀疏矩陣如圖2所示。

由圖2可以看出，該矩陣的(j,k)元素為50次模擬中協(xié)方差逆陣估計(jì)值(j,k)元素的非零頻率。

圖2說(shuō)明缺失比例越高,選出的邊越少，圖模型越稀疏。

圖2 不同缺失比例下的平均稀疏矩陣

5 結(jié) 語(yǔ)

基于懲罰似然框架提出自適應(yīng)MissGLasso方法以處理缺失數(shù)據(jù)情形下協(xié)方差逆陣的估計(jì)問(wèn)題，采用EM算法和GLasso算法進(jìn)行優(yōu)化求解。該方法可用于數(shù)據(jù)完全隨機(jī)缺失的情況。通過(guò)不同模型下的模擬實(shí)驗(yàn)結(jié)果顯示，自適應(yīng)MissGLasso方法的圖模型選擇性能及協(xié)方差估計(jì)結(jié)果較MissGLasso方法更優(yōu)。此外，對(duì)核黃素生產(chǎn)數(shù)據(jù)集的圖模型結(jié)構(gòu)學(xué)習(xí)實(shí)驗(yàn)同樣驗(yàn)證了自適應(yīng)MissGLasso方法具有良好的模型選擇性能。