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

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

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

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

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

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

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

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

0　引言

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

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

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

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

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

1　模型

1.1潛變量和probit模型

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

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

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

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

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

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

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

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

1.3變量選擇

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

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

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

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

2　算法

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

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

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

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

2.2計算步驟

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

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

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

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

3　模特卡羅模擬

3.1模型設(shè)定

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

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

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

3.2估計量

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

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

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

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

3.3評價標(biāo)準(zhǔn)

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

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

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

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

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

3.4結(jié)果分析

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

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

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

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

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

表1　當(dāng)樣本量為50時的模擬結(jié)果

表2　當(dāng)樣本量為100時的模擬結(jié)果

4　實證

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

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

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

表4　系數(shù)估計

5　結(jié)論

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

在很多實際問題中，會觀測到連續(xù)變量、二值變量等不同類型的數(shù)據(jù)，本文構(gòu)造了多元線性回歸來刻畫連續(xù)型響應(yīng)變量和二值響應(yīng)變量的聯(lián)合建模，并重點研究了變量選擇問題。對二值響應(yīng)引入了服從正態(tài)分布的潛變量，從而把問題轉(zhuǎn)化為多重響應(yīng)的多元線性回歸的變量選擇。本文在構(gòu)造目標(biāo)函數(shù)時，考慮了不同的響應(yīng)變量之間的協(xié)方差，選取了l1懲罰，通過交叉驗證的方法來選擇調(diào)節(jié)參數(shù)。在模擬研究中，考慮不同的樣本量以及不同結(jié)構(gòu)的設(shè)計陣和協(xié)方差矩陣，并與其他方法比較，本文提出的方法利用不同的響應(yīng)變量之間的信息，提高預(yù)測的精確性。表現(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]肖堯,張達(dá).匯率變動對經(jīng)濟(jì)增長的影響——基于發(fā)達(dá)國家與發(fā)展中國家的對比分析[J].金融教學(xué)與研究,2009,(06).

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

（責(zé)任編輯/易永生）

021

A

1002-6487（2016）19-0004-05

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