二值型響應與連續型響應聯合建模的變量選擇

胡亞南，張陶陶，田茂再

（1.中國人民大學統計學院；2.中國人民大學統計研究中心，北京100875）

二值型響應與連續型響應聯合建模的變量選擇

胡亞南，張陶陶，田茂再

（1.中國人民大學統計學院；2.中國人民大學統計研究中心，北京100875）

由于多重響應變量之間可能存在相關性，文章考慮對二值型響應變量和連續型響應變量進行聯合建模。利用probit模型，對二值響應引入了具有正態分布的潛變量，從而對多重響應建立線性回歸模型，能得到二值變量和連續變量的聯合分布。然后考慮回歸系數會存在稀疏性，通過對似然函數加懲罰，從而對二重響應的回歸系數和協方差矩陣的逆矩陣進行估計，達到參數估計和變量選擇的目標。文中目標函數基于l1懲罰。數值模擬和實證分析展示了所提出方法的良好性質。

EM算法；多元正態分布；Probit模型；聯合建模；LASSO；變量選擇

0　引言

在一個具體問題中，人們所關心的響應變量可能是在不同尺度下測量，既有二值響應變量，也有連續型響應變量。要看協變量對響應變量的影響，此時如果對二值變量和連續變量分別進行建模，會忽略不同的響應變量之間的相關性，從而丟失有用的信息。基于這種考慮，對二值響應變量和連續型響應變量進行聯合建模，從而有效利用樣本觀測信息，并更好的刻畫響應變量之間的相關性。

有很多關于二值變量和連續變量的聯合建模的方法被提出來。這些方法往往引入潛變量，對二值響應變量使用了probit模型，它引入服從正態分布的潛變量進行建模。Catalano&Ryan(1992)利用潛變量概念推導出了連續變量和離散變量的聯合分布，并把模型應用到聚類數據。他們把二重響應變量的聯合分布寫成連續變量的隨機效應模型與離散變量的probit模型乘積的形式。采用廣義估計方程(GEE)的方法來估計參數。Albert&Chib(1993)根據數據增廣的思想，利用精確貝葉斯的方法對類別響應變量進行建模分析。把二值響應的probit回歸模型看作連續型潛變量的正態回歸。由于潛變量能從合適的截斷正態分布中產生，一旦潛變量的實現值已知，那么參數的后驗分布可以從標準的線性模型結果中得出。數據增廣的方法為分析二值回歸模型提供了一般的結構。作者把probit模型應用到無序的多項響應變量和有序的多項響應變量，用貝葉斯的方法進行估計推斷。Dunson(2000)提出了一個靈活的方法來對混合變量進行貝葉斯分析。通過利用廣義線性模型來描述潛變量的聯合分布，模型能適用于更廣泛的數據結構。在文章所提出的結構下，新的模型可以推廣到聯合的二值變量、分類變量和連續變量。連續型響應變量和類別響應變量聯合建模的一個難點在于缺乏自然的多元分布。Gueorguieva&Agresti(2001)提出了相關的probit模型來對聚類的二值型響應和連續型響應進行聯合建模，他們對二值響應引入了服從正態分布的潛變量，并對這樣一個相關的probit模型進行研究。作者對引入的潛變量和連續變量同時建立線性混合效應模型，采用MCEM算法估計參數。Liu etal.(2009)對縱向的二值和連續過程進行聯合建模，并應用到戒煙試驗中，這兩個過程的相依性由無限制的回歸系數所刻畫；作者采用貝葉斯變量選擇來估計參數，尋找稀疏模型。Holstetal.(2015)處理這類聯合建模的問題，引入了潛變量和線性潛變量模型，提出了極大似然的估計方法，并且能分析含左刪失、右刪失的觀測數據。

對二值型響應和連續型響應的聯合建模，同時作回歸分析，然而回歸中的變量選擇是統計研究的熱點問題。選擇稀疏模型，不但能提高預測的精確性，而且更好解釋。隨著大數據時代來臨，高維數據越來越普遍，諸如最優子集等傳統的變量選擇方法在面臨這些數據時，由于計算量太大，往往無法滿足需求。基于懲罰函數的變量選擇方法越來越受到統計學者的關注。這類方法是在最小二乘或極大似然目標函數上加上或者減去懲罰函數而得到新的目標函數，然后最優化目標函數，進而得到參數的估計。這種方法的優點在于參數估計和變量選擇同時進行，大大提高了計算速度。

Tibshirani(1996)提出了lasso的方法，通過對回歸系數作l1范數的懲罰，壓縮系數，把一些絕對值較小的系數壓縮為0，從而達到估計參數和變量系數的目的。lasso方法克服了傳統變量選擇方法的不足，在統計領域受到了極大關注。繼而，lasso開始應用到其他模型中，并且也有很多文獻對lasso進行改進。Tibshirani(1997)把lasso方法應用到生存分析領域，對Cox比例風險模型做變量選擇。Zou(2006)對lasso方法做了改進，提出了自適應lasso，即對不同的回歸系數施加不同的權重的懲罰，所得到的估計量具有良好的性質，并且這種方法具有0 racle性質。這些研究都是把lasso應用到單一響應變量的情形。Turlach (2005)把lasso擴展到多重響應變量的情形，通過對回歸系數加懲罰，選擇共同的解釋變量。對于多重響應的變量選擇問題，Simon etal.(2013)提出了區塊降速算法來求解加group懲罰的目標函數，得到了回歸系數的系數估計，但是沒有考慮響應變量之間的相依性。Friedman etal.(2008)以多元正態分布為研究對象，用圖lasso的方法得到協方差矩陣的逆矩陣的估計，簡化了概率圖模型的結構。Rothman etal.(2010)研究了多重響應的回歸分析，既構建回歸系數矩陣的稀疏估計量，同時又考慮了響應變量之間的相關性，通過最優化加懲罰的似然函數，得到回歸系數和協方差結構的估計。

在實證分析部分，本文主要研究了國內生產總值（連續型響應變量）和是否為發達國家（二值響應變量）的聯合建模。國內生產總值和是否為發達國家，作為衡量國家經濟發展和評價綜合國力的重要指標，都是經濟研究中的重要課題，但在以往的研究中，多以其一為響應變量，考慮其影響因素進行建模。冶濤（2012）以固定資產投資總額、財政收入等六個解釋變量建立GDP的多元回歸模型，肖堯等（2009）研究匯率變動對經濟增長的影響在發達國家與發展中國家的對比分析。兩者分別的研究已經較為成熟，但本文考慮到國內生產總值和國家發達水平之間的相關性，對兩者進行聯合建模；然后,對二值響應引入服從正態分布的潛變量,然后對連續性變量和潛變量的聯合分布進行建模。為了得到回歸系數的稀疏估計,同時利用響應變量之間的相依性信息,我們對目標函數加自適應lasso的懲罰。由于二值型響應變量和連續型響應變量之間會存在響應性，我們考慮聯合建模；然后，對二值響應引入服從正態分布的潛變量，然后對連續性變量和潛變量的聯合分布進行建模。為了得到回歸系數的稀疏估計，同時利用響應變量之間的相依性信息，我們對目標函數加自適應lasso的懲罰。

1　模型

1.1潛變量和probit模型

在二重響應變量的聯合建模中，二值型響應變量的存在很大程度上增加了建模的難度，此時潛變量提供了一個實用且直觀的方法來對離散型響應變量進行建模。本文二值型響應變量的probit模型引入了潛變量，即模型事先假定一個不可觀測的連續型隨機變量存在，并且潛變量超過一定的門限值時，假設一個二值事件發生。在對多重響應建模之前，我們首先回顧一下單變量的情形。考慮到線性模型

其中Y1i表示響應變量，是協變量，其中β0表示截距項。?i是誤差項，且?i～N(0,σ2)。

當觀測數據Zi是二值響應變量，與潛變量Y1i滿足關系：由線性模型(1)誤差項的假設，可以得出Zi滿足probit模型

Φ(?)是正態分布的累積分布函數。一則，潛變量在具體應用中的體現；二則從統計角度來看，潛變量具有非常吸引人的地方，因為這樣的假設下，二值響應的正態模型，有一個非常方便的形式。

1.2二重響應變量的聯合模型

假設有n個觀測，第i個觀測的響應變量為(Zi,Y2i)，其中Zi是二值型變量，我們引入潛變量Y1i，則Zi與Y1i滿足方程(2)；Y1i是連續變量。對二值型變量和連續型變量聯合建模，那么相關的probit模型如下：

如果ρ=0,則Σ退化為對角矩陣。對二值響應Zi，引入了潛變量了Y1i來建模，即使用了probit模型.記Yi=(Y1i,Y2i),對方程(4),可改寫為

其中B=(β1,β2)，是系數矩陣。

1.3變量選擇

統計學習中，有兩個基本的目標，一則預測的精確性；二則找到相關的協變量，從而方便解釋.當真實的模型有稀疏表示時，變量選擇尤為重要。從式(4)中可以看出，當引入潛在變量之后，連續型響應和二值型響應的聯合建模問題，轉變為多重響應的線性回歸模型，令X是n×p的設計矩陣，Y的n×2的響應變量矩陣，由于誤差?1,…,?n是獨立同分布于N2(0,Σ)，那么在給定X的情況下，多重響應的協方差陣為Σ。

由于不同的響應之間存在相關性，所以E所對應的協方差陣不是對角矩陣。為了方便起見，記Ω=Σ-1，那么模型(5)的對數似然表示：

由于要對多重回歸模型作變量選擇，則考慮如下對系數懲罰的似然函數

其中βjk是系數矩陣B中的元素，λ是調節參數。通過最優化目標函數從而得到參數的估計

2　算法

由于對二值響應引入了潛變量，而潛變量是不可觀測的。Dempster etal.(1977)提出了EM(expectation-maximization)算法。對于解決含缺失數據、潛變量等不完整數據，EM算法是一種行之有效的方法.Gueorguieva&Agresti (2001)在處理聚類的二值響應和連續型響應聯合建模問題時，利用改進的EM算法得到相關probit模型的極大似然估計。本文也采用EM算法，用潛變量的條件期望代替潛變量。

2.1似然函數和條件分布

利用數據增廣技術，引入潛變量Y1i，根據方程(2)和(4),得到完全數據(Y1i,Y2i,Zi)的概率密度函數

在已知聯合分布(10)情形下，公式(11)可由簡單的數值積分求解。

2.2計算步驟

要得到參數的估計，需要對目標函數(8)最大化。Rothman etal.(2010)給出了優化過程中的計算細節，這里只給出大致的計算步驟：

（2）E-步：由于Y1i是潛變量，無法觀測到，但是其分布是知道的，因此利用它的條件分布，用期望值代替。基于當前的參數估計值根據條件分布用條件期望來代替潛變量Y1i；

考慮到EM算法是尋找的局部最優解，迭代過程需要設定模型參數的初始值其初始值的選取，利用分別建模的方法。

2.3選擇調節參數

3　模特卡羅模擬

3.1模型設定

本文并未考慮對系數矩陣B的選擇，只是簡單設置滿足稀疏性。令

生成n×p的協變量矩X，每個行的觀測Xi獨立同分布于Np(0,ΣX)，其中所有協變量的邊際方差為1。誤差矩陣的行向量?i來自于正態分布其中

根據方程(4)，可以得到響應(Y1i,Y2i)的取值，由于Y1i是潛變量，根據方程(2)，當Y1i≥0時，Zi=1;當Y1i＜0時，Zi=0。則(Zi,Y2i)是要進行分析的觀測。

3.2估計量

首先，采用數據增廣技術，對二值響應Zi引入潛變量Y1i，然后用用條件期望E(Y1i|Y2i,Zi)來代替潛變量Y1i。在此基礎上，為了比較所提出方法的表現，我們設置了對照模型。模擬中所展示的模型有：

模型1：對二值響應和連續型響應分別建模，不考慮變量選擇；

模型2：對二值響應和連續型響應同時建模，不考慮協方差結構；

模型3：對二值響應和連續型變量同時建模，考慮協方差結構。

3.3評價標準

我們從參數的均方誤差和稀疏指標選擇兩個角度來評估模型的好壞。

定義參數的均方誤差為：

βj表示連續型響應變量或者二值響應變量所對應的回歸系數向量，表示第i次重復所得到的參數估計向量，d1+1表示參數的維數。MSE就小，說明模型的估計效果越好。

度量稀疏性的指標包括敏感性(Sensitivity)和特異性(Specificity)：

其中#表示計數。敏感性和特異性在0到1之間，越接近1，說明變量選擇的效果越好。

3.4結果分析

（1）模型1利用一般線性模型擬合，并沒有作變量選擇，其結果展示敏感性全為0，特異性全為1，沒有把有效的變量篩選出來，這在預期之中。

（2）在其他設置不變的情況下，隨著樣本量的增加，模型1-模型3的MSE變小；整體看來，模型2和模型3的特異性和敏感性變大，即變量選擇的效果越好。

（3）在設置相同的情況下，連續型響應部分和二值響應部分，模型3的MSE較小，敏感性和特異性較大，這也反映了本文所提方法的優越性及合理性。一種我們所提出的方法在所有準則下表現最好。這表明，對所有組聯合建模有助于提高預測精確性和估計。

（4）整體上看，連續型響應部分和二值響應部分比較，二值響應部分的敏感性和特異性較小，即變量選擇效果較差。

（5）在其他設置不變的情況下，隨著ρe的增加，模型3的敏感性和特異性增加。

表1　當樣本量為50時的模擬結果

表2　當樣本量為100時的模擬結果

4　實證

本文實證分析部分主要研究了國內生產總值和國家發達程度的聯合建模。本文的數據來源為《國際統計年鑒》，除響應變量為國內生產總值和發達程度以外，還考慮的協變量包括：對外直接投資、外商直接投資、貨物進口總額、貨物出口總額、資本形成率、居民消費率、發電量七個經濟類指標，以及森林資源、淡水資源、國土面積、二氧化碳排放量四個環境資源類指標。對所有變量信息匯總如表3。

運用本文方法對于二重響應變量進行聯合建模，得到參數估計結果如表4。可以看出，6個協變量對發達程度有顯著的影響，包括：對外直接投資、森林資源、國土面積、二氧化碳排放量、貨物進口總額、貨物出口總額；3個協變量對國內生產總值有顯著的影響，包括：對外直接投資、發電量和貨物進口總額。

表3　變量額匯總統計

表4　系數估計

5　結論

多元回歸是解決實際問題的一個常用工具。許多多元回歸技術是為單個響應的情況設計的。對于多重響應變量的情況，一個通常的方法是應用單個響應變量的回歸技術，分別對每個響應變量作回歸分析。盡管這樣很簡單、也很流行，但是這樣處理，會忽略不同的響應變量之間的聯合信息。

在很多實際問題中，會觀測到連續變量、二值變量等不同類型的數據，本文構造了多元線性回歸來刻畫連續型響應變量和二值響應變量的聯合建模，并重點研究了變量選擇問題。對二值響應引入了服從正態分布的潛變量，從而把問題轉化為多重響應的多元線性回歸的變量選擇。本文在構造目標函數時，考慮了不同的響應變量之間的協方差，選取了l1懲罰，通過交叉驗證的方法來選擇調節參數。在模擬研究中，考慮不同的樣本量以及不同結構的設計陣和協方差矩陣，并與其他方法比較，本文提出的方法利用不同的響應變量之間的信息，提高預測的精確性。表現出了一定的優勢。

[1]Albert JH,Chib S.Bayesian Analysis of Binary and Polychotomous Response Data[J].Journal of the American Statistical Association, 1993,88(422).

[2]Catalano P J,Ryan L M.Bivariate Latent Variable Models for Clus?tered Discrete and Continuous Outcomes[J].Journal of the American Statistical Association,1992,87(419).

[3]Dunson D B.Bayesian Latent Variable Models for Clustered Mixed Outcomes[J].Journalof the Royal Statistical Society.Series B,Statis?ticalMethodology,2000.

[4]Friedman J,Hastie T,Tibshirani R.Sparse Inverse Covariance Esti?mationWith the Graphical Lasso[J].Biostatistics,2008,9(3).

[5]Gueorguieva R V,AgrestiA.A Correlated ProbitModel for JointMod?eling of Clustered Binary and Continuous Responses[J].Journal of the American StatisticalAssociation,2001,96(455).

[6]Holst K K,Budtz-Jorgensen E,Knudsen GM.A IatentVariableMod?elWith Mixed Binary and Continuous Response Variables[J].Staist?ics,2015.

[7]Liu X,DanielsM J,Marcus B.JointModels for the Association of Lon?gitudinal Binary and Continuous Processes With Application to a Smoking Cessation Trial[J].Journal of the American Statistical Asso?ciation,2012.

[8]Rothman A J,Levina E,Zhu J.Sparse Multivariate Regression With Covariance Estimation[J].Journal of Computational and Graphical Statistics,2010,19(4).

[9]Simon N,Friedman J,Hastie T.A Blockwise Descent Algorithm for Group-penalized Multiresponse and Multinomial Regression[J].Sta?tistics,2013.

[10]Tibshirani R.Regression Shrinkage and Selection via the Lasso[J]. Journalof the Royal Statistical Society.Series B(Methodological).

[11]Tibshirani R.The Lasso Method for Variable Selection in the Cox Model[J].Statistics in Medicine,1997,16(4).

[12]Turlach B A,VenablesW N,Wright S J.Simultaneous Variable Se?lection[J].Technometrics,2005,47(3).

[13]Zou H.The Adaptive Lasso and Its Oracle Properties[J].Journal of the American StatisticalAssociation,2006,101(476).

[14]肖堯,張達.匯率變動對經濟增長的影響——基于發達國家與發展中國家的對比分析[J].金融教學與研究,2009,(06).

[15]冶濤.國內生產總值影響因素實證分析——以新疆GDP增長因素為例[J].金融經濟,2012,(16).

（責任編輯/易永生）

021

A

1002-6487（2016）19-0004-05

國家自然科學基金資助項目(11271368)；國家社會科學基金重點項目(13AZD064)；教育部哲學社會科學研究重大課題攻關項目(15JZD015)；北京市社會科學基金重大項目(15ZDA17)；教育部高等學校博士學科點專項科研基金(20130004110007)；教育部人文社會科學重點研究基地重大項目(15JJD910001)；中國人民大學科學研究基金資助項目(15XNL008)