胡亞南，張?zhí)仗?，田茂?/p>

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

二值型響應與連續(xù)型響應聯(lián)合建模的變量選擇

胡亞南，張?zhí)仗?，田茂?/p>

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

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

EM算法；多元正態(tài)分布；Probit模型；聯(lián)合建模；LASSO；變量選擇

0　引言

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

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

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

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

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

1　模型

1.1潛變量和probit模型

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

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

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

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

1.2二重響應變量的聯(lián)合模型

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

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

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

1.3變量選擇

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

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

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

其中βjk是系數(shù)矩陣B中的元素，λ是調節(jié)參數(shù)。通過最優(yōu)化目標函數(shù)從而得到參數(shù)的估計

2　算法

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

2.1似然函數(shù)和條件分布

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

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

2.2計算步驟

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

（2）E-步：由于Y1i是潛變量，無法觀測到，但是其分布是知道的，因此利用它的條件分布，用期望值代替?；诋斍暗膮?shù)估計值根據(jù)條件分布用條件期望來代替潛變量Y1i；

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

2.3選擇調節(jié)參數(shù)

3　模特卡羅模擬

3.1模型設定

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

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

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

3.2估計量

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

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

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

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

3.3評價標準

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

定義參數(shù)的均方誤差為：

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

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

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

3.4結果分析

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

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

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

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

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

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

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

4　實證

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

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

表3　變量額匯總統(tǒng)計

表4　系數(shù)估計

5　結論

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

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

[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]肖堯,張達.匯率變動對經濟增長的影響——基于發(fā)達國家與發(fā)展中國家的對比分析[J].金融教學與研究,2009,(06).

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

（責任編輯/易永生）

021

A

1002-6487（2016）19-0004-05

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