Cox模型的SCAD變量選擇方法及其應(yīng)用

郎 英，苗新利，顧 凱

（楚雄師范學(xué)院數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，云南楚雄 675000）

英國(guó)著名統(tǒng)計(jì)學(xué)家Cox.aR[1]（1972）最早提出了Cox比例風(fēng)險(xiǎn)模型，它是生存分析中一種重要的半?yún)?shù)模型，可以處理分布無規(guī)律、存在刪失情況的數(shù)據(jù)。近年來Cox模型逐漸成為一種新的度量信貸違約風(fēng)險(xiǎn)的工具，它通過研究與貸款者違約相關(guān)的可能因素，建立預(yù)警模型。

在統(tǒng)計(jì)建模中，變量選擇是一個(gè)非常重要的環(huán)節(jié)。主要包含子集選擇和系數(shù)壓縮兩種方法。子集選擇法是在一個(gè)模型中，把所有變量表示成一個(gè)集合，然后根據(jù)相應(yīng)準(zhǔn)則從所有子集中選出一個(gè)最優(yōu)子集。然而，隨著變量個(gè)數(shù)的增加，子集選擇法不但計(jì)算量龐大且都缺乏一定的穩(wěn)定性。

系數(shù)壓縮方法的本質(zhì)是通過優(yōu)化一個(gè)“損失函數(shù)”+“懲罰函數(shù)”的形式來實(shí)現(xiàn)的。嶺回歸方法是一個(gè)連續(xù)的過程，可以收斂系數(shù)值，但卻不能將系數(shù)收斂至0，故不能實(shí)現(xiàn)變量選擇。Tibshiran（1996）在嶺回歸的基礎(chǔ)上進(jìn)行改進(jìn)，提出了Lasso（the Least Absolute Shrinkageand Selection Operator）方法[2]，能同時(shí)實(shí)現(xiàn)變量的選擇及其參數(shù)的估計(jì)，該懲罰函數(shù)有較好的穩(wěn)定性，但由于Lasso方法不具有Oracle性質(zhì)，后又有很多學(xué)者在Lasso基礎(chǔ)上進(jìn)行改進(jìn)。

把Cox模型運(yùn)用到信用風(fēng)險(xiǎn)度量研究中，以房地產(chǎn)上市公司為研究對(duì)象，結(jié)合上市公司的公開報(bào)表數(shù)據(jù)建立模型。用LLA算法和CD算法對(duì)帶有SCAD懲罰的對(duì)數(shù)似然函數(shù)進(jìn)行計(jì)算，再運(yùn)用BIC準(zhǔn)則選取合適的懲罰參數(shù)，旨在通過SCAD變量選擇方法得到影響上市公司生存期長(zhǎng)短的重要協(xié)變量，建立Cox模型并與基于傳統(tǒng)逐步回歸方法建立的Cox模型進(jìn)行比較。

1 模型

1.1 Cox模型

Cox模型主要是研究生存時(shí)間與協(xié)變量之間的關(guān)系。其一般形式可表示為：

其中h(t|z)表示在給定協(xié)變量z=(z1,z2,…,zk)T的條件下，生存時(shí)間t的風(fēng)險(xiǎn)函數(shù)，β=(β1,β2,…,βk)T是待估計(jì)的參數(shù)，且 βi＞0表示協(xié)變量zi為危險(xiǎn)因素，βi＜0表示zi為保護(hù)因素，h0(t)表示未知基準(zhǔn)風(fēng)險(xiǎn)函數(shù)。

假定有n個(gè)觀測(cè)樣本，觀測(cè)數(shù)據(jù)為(wi,δi,zi)(i=1,2,…,n)，其中ωi=min(Ti,Ci)為第i個(gè)個(gè)體的研究時(shí)間，Ti為研究對(duì)象的生存時(shí)間，Ci為隨機(jī)刪失時(shí)間，δi=I(Ti≤Ci)為示性函數(shù)。當(dāng)δi=1表示完整數(shù)據(jù)；當(dāng)δi=0表示刪失數(shù)據(jù)。zi=(zi1,zi2,…,zik)T為第i個(gè)個(gè)體的k維協(xié)變量。

則Cox模型的偏似然函數(shù)可以表示為：

故偏對(duì)數(shù)似然函數(shù)為：

1.2 基于SCAD方法的Cox模型

Fan and Li（2001）在Lasso方法的基礎(chǔ)上進(jìn)行了相應(yīng)的改進(jìn)，提出了一種非凸懲罰函數(shù)絕對(duì)值偏差的光滑剪切懲罰(The Smooth integration of Counting and Absolute deviation，簡(jiǎn)稱SCAD)[3]，與Lasso方法對(duì)所有系數(shù)都進(jìn)行壓縮不同，SCAD方法只對(duì)模型中不顯著的變量系數(shù)進(jìn)行壓縮，同時(shí)保留顯著變量。將(3)式的相反數(shù)極小化，再加上SCAD罰函數(shù)，并表示成無約束條件形式：

由于SCAD罰函數(shù)是一個(gè)分段函數(shù)，且在(-∞,0)?(0,+∞)上是連續(xù)可微的，故 ρλ定義為：

這里a=3.7[4],λ＞0為調(diào)整參數(shù)。則Cox模型的SCAD估計(jì)量可表示為如下形式：

用BIC準(zhǔn)則求解λ。其定義如下：

其中L為Cox模型的部分偏對(duì)數(shù)似然函數(shù)，b為模型中不為零的變量的個(gè)數(shù)。BIC準(zhǔn)則的基本思想是在所有可能的自變量集合中，找到BIC值最小的一組自變量組合。

1.3 SCAD估計(jì)量的迭代算法

SCAD懲罰函數(shù)在原點(diǎn)沒有連續(xù)的二階導(dǎo)函數(shù)。局部線性逼近算法LLA是針對(duì)非凸懲罰函數(shù)在原點(diǎn)不可導(dǎo)而提出的，針對(duì)非凸懲罰函數(shù)提出了坐標(biāo)下降算法，其計(jì)算速度比LLA算法更快。首先采用LLA算法對(duì)SCAD罰函數(shù)進(jìn)行局部線性逼近[5]。假設(shè)初始值為 β0,若 βj0近似等于0,則令βj0=0;若 βj0≠0,則SCAD罰函數(shù) pλ(| |βj)可用如下線性函數(shù)來近似，當(dāng) βj≈βj0時(shí)有

參考Zhang and Liu[6]對(duì)Cox模型的極大似然函數(shù)的近似處理，通過牛頓迭代更新公式，把Cox模型的對(duì)數(shù)似然函數(shù)轉(zhuǎn)化為線性回歸問題。分別定義梯度向量 ?l(β)和海瑟矩陣 ?2l(β)為：

然后根據(jù)cholesky分解可得

令

則根據(jù)泰勒展開對(duì)l(β)進(jìn)行局部二次逼近有：

結(jié)合對(duì)數(shù)似然函數(shù)的二階近似和SCAD懲罰函數(shù)的LLA近似，原函數(shù)可轉(zhuǎn)化為如下形式

將上式化簡(jiǎn)，刪除一些常數(shù)項(xiàng)后可轉(zhuǎn)化為：

其中,

此時(shí)原問題就轉(zhuǎn)化成了一個(gè)基于最小二乘的自適應(yīng)Lasso估計(jì)問題[7]。這時(shí)可采用坐標(biāo)下降（CD）算法來計(jì)算SCAD估計(jì)量的值[8]，具體算法步驟如下：

Step1，算出負(fù)對(duì)數(shù)似然函數(shù)的最小二乘估計(jì)值，令這個(gè)值為初始估計(jì)值0；

2 實(shí)證分析

2.1 樣本選取及數(shù)據(jù)的預(yù)處理

信用風(fēng)險(xiǎn)（credit risk）又稱違約風(fēng)險(xiǎn)，由于銀行主要是通過向企業(yè)發(fā)放貸款來獲得本金及利息，故信用風(fēng)險(xiǎn)將是商業(yè)銀行所面臨的最大風(fēng)險(xiǎn)。以滬、深A(yù)股的所有房地產(chǎn)上市公司為研究對(duì)象，鑒于違約公司的數(shù)據(jù)不易獲得，故把因財(cái)務(wù)狀況異常而被貼別處理的（ST和*ST）房地產(chǎn)上市公司視為違約樣本，沒有被特別處理的視為非違約樣本[9]。

所用數(shù)據(jù)均來源于wind數(shù)據(jù)庫(kù)。根據(jù)證監(jiān)會(huì)行業(yè)分類，2015年一共有126家房地產(chǎn)上市公司。樣本觀測(cè)時(shí)間為：2005年1月至2015年12月，以年為時(shí)間尺度。選擇樣本公司T-3年的數(shù)據(jù)建立模型，最后選擇了121個(gè)樣本公司，在建模時(shí)按相應(yīng)比例對(duì)樣本進(jìn)行了分組，其中訓(xùn)練樣本83個(gè)（完整數(shù)據(jù)22個(gè)，刪失數(shù)據(jù)61個(gè)），檢驗(yàn)樣本38個(gè)（完整數(shù)據(jù)11個(gè)，刪失數(shù)據(jù)27個(gè)）。

考慮到數(shù)據(jù)的可獲得性，共選取了28個(gè)初始變量，主要有反應(yīng)公司每股情、盈利能力、償債能力、發(fā)展能力、運(yùn)營(yíng)能力、及資本結(jié)構(gòu)等幾個(gè)方面的指標(biāo)，具體指標(biāo)見表1。

采用了均值插補(bǔ)的方法對(duì)缺失數(shù)據(jù)進(jìn)行處理，為了消除不同量綱對(duì)模型結(jié)果解釋性的干擾，在SPSS中采用Z-Score方法對(duì)數(shù)據(jù)進(jìn)行標(biāo)了準(zhǔn)化處理。

表1 指標(biāo)符號(hào)定義

2.2 模型的建立

通過計(jì)算協(xié)變量的相關(guān)系數(shù)矩陣，發(fā)現(xiàn)z13與z10的相關(guān)系數(shù)大于0.7,且z13與z23的相關(guān)系數(shù)大于了0.9,說明變量之間存在嚴(yán)重的共線性。本文采用了基于SCAD懲罰的方法，同時(shí)進(jìn)行變量選擇和參數(shù)估計(jì)。

在用SCAD罰函數(shù)進(jìn)行變量選擇時(shí)，可調(diào)用R軟件中的scadcox(.)函數(shù)，再通過BIC準(zhǔn)則求得懲罰參數(shù)λ，結(jié)果如圖1，圖中橫坐標(biāo)為λ的對(duì)數(shù)值，從圖中可以看到篩選出的變量個(gè)數(shù)隨λ取值的不同而不同，得到最小的BIC值下lnλ=-3.2，對(duì)應(yīng)λ值為0.04,這時(shí)選擇了10個(gè)變量，按照其絕對(duì)值大小分別為：z1,z5,z20,z28,z3,z18,z21,z27,z15,z26。SCAD算法的系數(shù)解路徑如圖2。

圖1 λ的對(duì)數(shù)值

圖2 SCAD系數(shù)解路徑

通過調(diào)用R語(yǔ)言中的coxph(.)函數(shù)，得到進(jìn)入模型的10個(gè)變量及其對(duì)應(yīng)的回歸參數(shù)如表2所示：

表2 SCAD變量選擇結(jié)果

由表2最后一列可得z26的p值大于0.1，變量z26沒有通過顯著性檢驗(yàn)，即對(duì)模型沒有顯著影響，故建模時(shí)刪除變量z26。對(duì)所選的9個(gè)變量做共線性檢驗(yàn)，得到所有變量的VIF值均比較小，也即是說變量之間不存在多重共線性影響，故用這9個(gè)變量建立Cox模型，得

這時(shí)通過代入相應(yīng)的樣本數(shù)據(jù)可計(jì)算出每個(gè)房地產(chǎn)上市公司的生存概率。

2.3 模型檢驗(yàn)

在用Cox方法建模時(shí)需進(jìn)行PH假定檢驗(yàn)，選擇了schoenfeld提出的殘差圖檢驗(yàn)方法[10]。其定義為

若模型滿足PH假定，則E(rik)≈0,這時(shí)對(duì)ti作圖時(shí)應(yīng)以0為中心隨機(jī)分布；若所繪圖形傾斜趨勢(shì)較大，說明不滿足PH假定。9個(gè)變量的非零斜度統(tǒng)計(jì)量值可通過調(diào)用R中的cox.phz(.)函數(shù)得到，結(jié)果如表3所示：

表3 PH假定檢驗(yàn)

表2中的rho表示生存時(shí)間與schoenfeld殘差的相關(guān)系數(shù)，chiq表示卡方檢驗(yàn)的值，從表中可看出每個(gè)變量的p值均大于0.1,且全局檢驗(yàn)的p值也大于0.1,說明該模型滿足PH假定。

對(duì)參數(shù)β采用似然比檢驗(yàn)。令模型參數(shù)β=(β1,β2,…βp)T,則根據(jù) LR 原理在原假設(shè)H0:β=(β1,β2,…,βk)T=0,(k＜p)成立的條件下，似然比 LR可表示為：

其中L()為似然函數(shù)，k為加入的解釋變量的個(gè)數(shù)。通過調(diào)用R軟件中的summary()函數(shù)，得到模型參數(shù)的似然比值如表4。

2.4 模型預(yù)測(cè)

通過以上建立好的模型對(duì)所選的11家ST公司及27家非ST公司進(jìn)行預(yù)測(cè)，并與傳統(tǒng)逐步回歸方法建立的Cox模型進(jìn)行比較，可得判別結(jié)果如表5。

表4 基于SCAD的模型顯著性檢驗(yàn)

從表5可看出基于SCAD方法的Cox模型的擬合效果是較好的，在38個(gè)檢驗(yàn)樣本中正確識(shí)別了32個(gè)，總體正確率為0.842，而逐步回歸方法總體正確率只有0.736。

表5 模型預(yù)測(cè)效果

綜上可知，Cox模型在信用風(fēng)險(xiǎn)度量中有較高的應(yīng)用價(jià)值，不僅可以動(dòng)態(tài)預(yù)測(cè)企業(yè)的財(cái)務(wù)狀況，還可以預(yù)測(cè)企業(yè)的生存時(shí)間。基于SCAD變量選擇的Cox模型對(duì)房地產(chǎn)上市公司的信用風(fēng)險(xiǎn)有較高的識(shí)別準(zhǔn)確率，可以為商業(yè)銀行的信用風(fēng)險(xiǎn)度量提供良好的依據(jù)，有一定的參考價(jià)值。

3 結(jié)論

把Cox模型運(yùn)用到信用風(fēng)險(xiǎn)度量中，根據(jù)我國(guó)房地產(chǎn)上市公司的特點(diǎn)，選取影響公司經(jīng)營(yíng)狀況的財(cái)務(wù)指標(biāo)建立模型。由于傳統(tǒng)Cox比例風(fēng)險(xiǎn)模型主要應(yīng)用于維數(shù)較低的數(shù)據(jù)環(huán)境中，當(dāng)變量維數(shù)較高時(shí)所得結(jié)果可能有偏差。針對(duì)這一問題，把SCAD懲罰變量選擇方法運(yùn)用到Cox模型中，介紹了Cox模型下SCAD方法的定義，構(gòu)造了β的懲罰估計(jì)量，給出了一種有效的迭代算法。

根據(jù)所選樣本，建立了基于SCAD懲罰的Cox模型，通過PH假定檢驗(yàn)及參數(shù)顯著性檢驗(yàn)，得到所建模型滿足要求。最后利用判斷準(zhǔn)確率對(duì)比分析了基于SCAD方法的Cox模型及傳統(tǒng)逐步回歸方法的Cox模型，結(jié)果表明基于SCAD方法建立的Cox模型的預(yù)測(cè)準(zhǔn)確率高于逐步回歸方法建立的模型，更適合我國(guó)房地產(chǎn)上市公司的信用風(fēng)險(xiǎn)研究。