基于自適應Lasso的缺失數據圖模型選擇

鄭倩貞， 徐平峰

(長春工業大學 數學與統計學院， 吉林 長春 130012)

0 引 言

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

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

1 自適應MissGLasso

設x=(xobs,xmis)為來自多元高斯分布Np(μ,Ω-1)的n個獨立同分布樣本。令樣本中的觀測數據xobs=(xobs,1,…,xobs,n)，缺失數據xmis=(xmis,1,…,xmis,n)，其中xobs,i和xmis,i分別表示第i個樣本的觀測數據和缺失數據的集合，i=1,2,…,n。與文獻[4]相同，文中假設數據的缺失機制為可忽略的。基于觀測數據xobs，可給出觀測對數似然函數

l(μ,Ω;xobs)=

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

(1)

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

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

在懲罰似然的框架下，下述優化問題的解可作為均值μ和協方差逆陣Ω的估計，

(2)

式中：λ----正則化參數;

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

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

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

(3)

我們稱該估計為自適應MissGLasso(AdaMissGLasso)。針對上述優化問題，可結合EM算法[9]和GLasso算法[3]進行優化求解。

2 EM算法

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

1)E步。

給定觀測數據xobs和當前參數(μ(t),Ω(t))，計算完全數據懲罰對數似然函數的條件期望，記為Q函數。為簡單起見，記E(·|xobs,μ(t),Ω(t))為Et(·)。于是，Q函數可表示為

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

(4)

其中

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

(5)

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

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

給定觀測數據xobs,i和當前參數(μ(t),Ω(t))，xmis,i的條件分布為N(ci,σi)，其中

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

(6)

以及

式中：Rij=1----xij為觀測數據；

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

2)M步。

M步通過最小化Q(μ,Ω|μ(t),Ω(t))來更新參數估計(μ(t+1),Ω(t+1))，

(8)

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

m=1,2,…,k，

(9)

其中

3 模擬研究

在完全隨機缺失的機制下，文中對自適應MissGLasso方法和MissGLasso方法在高斯圖模型上的參數估計和模型選擇進行了模擬比較。現考慮以下兩個模型：

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

2)Ω=B+δI，其中B的對角線元素為0，非對角線元素之間相互獨立，取值為 0.5或0，且P(Bjk=0.5)=0.1。選擇恰當的δ，使得Ω的條件數為p。

對于以上兩個模型，設定隨機變量個數p=100,200,300，對兩個模型與隨機變量個數p的6種組合，均產生樣本量n=200的50個獨立的數據集，且在每個數據集上完全隨機刪除10%、20%、30%的數據。因此，每個完全觀測的數據集對應3個具有不同缺失比例的缺失數據集。為評價算法性能，文中比較以下指標。

Kullback-Leibler損失(kl)、真陽率(tpr)、陽性預測率(ppv)以及馬修斯相關系數(mcc)。

式中:Ω----真實的協方差逆陣;

tp----真陽類的個數;

tn----真陰類的個數;

fp----假陽類的個數;

fn----假陰類的個數。

基于模型1和模型2的不同參數設定各評價指標的均值和標準差分別見表1和表2。

表1 不同參數設定下模型1各評價指標的均值和(標準差)

表2 不同參數設定下模型2各評價指標的均值和(標準差)

由表中可以看出，對于模型1和模型2的所有設定，除了tpr以外，自適應MissGLasso方法均優于MissGLasso，并且模型1的tpr僅在缺失比例為30%時略低。對于模型2，自適應MissGLasso方法的mcc略低，但是兩者相差不大。kl為評價協方差逆陣估計效果的指標，不論模型1還是模型2，自適應MissGLasso的估計性能都更佳。綜上所述，自適應的方法更優。

4 實例分析

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

不同缺失程度下重疊邊數的箱線圖如圖1所示。

圖1 不同缺失比例下的重疊邊數

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

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

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

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

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

5 結 語

基于懲罰似然框架提出自適應MissGLasso方法以處理缺失數據情形下協方差逆陣的估計問題，采用EM算法和GLasso算法進行優化求解。該方法可用于數據完全隨機缺失的情況。通過不同模型下的模擬實驗結果顯示，自適應MissGLasso方法的圖模型選擇性能及協方差估計結果較MissGLasso方法更優。此外，對核黃素生產數據集的圖模型結構學習實驗同樣驗證了自適應MissGLasso方法具有良好的模型選擇性能。