高仙立，姜玉英

（1.首都經(jīng)濟(jì)貿(mào)易大學(xué) 統(tǒng)計學(xué)院，北京 100070；2.北京印刷學(xué)院 基礎(chǔ)部，北京 102600）

縱向數(shù)據(jù)下半?yún)?shù)Logistic模型的變量選擇

高仙立1，姜玉英2

（1.首都經(jīng)濟(jì)貿(mào)易大學(xué) 統(tǒng)計學(xué)院，北京 100070；2.北京印刷學(xué)院 基礎(chǔ)部，北京 102600）

文章研究了縱向數(shù)據(jù)半?yún)?shù)Logistic回歸模型的估計問題，給出了模型中未知參數(shù)和未知函數(shù)的估計方法，探討了參數(shù)部分的變量選擇問題，并對不同的變量選擇方法進(jìn)行比較分析。從模擬結(jié)果可以看到，文中給出的方法具有很好的估計效果。

縱向數(shù)據(jù)；半?yún)?shù)Logistic模型；變量選擇

0 引言

縱向數(shù)據(jù)指同一研究對象在不同時點上的觀測數(shù)據(jù)，與計量經(jīng)濟(jì)學(xué)中廣泛應(yīng)用panel數(shù)據(jù)不同，縱向數(shù)據(jù)不要求所有研究對象的觀測時點相同，因而在實際應(yīng)用中更為普遍。對于縱向數(shù)據(jù)參數(shù)回歸模型的估計問題，目前已經(jīng)有了一些研究成果，如Diggle等（2002），F(xiàn)u和Wang（2016）等。當(dāng)前，隨著實驗技術(shù)、檢驗手段的日益提高，所獲數(shù)據(jù)的變量個數(shù)越來越多，變量選擇顯得尤為重要。關(guān)于變量選擇的研究一直是眾多統(tǒng)計學(xué)者的研究熱點，學(xué)者們提出了一系列標(biāo)志性的變量選擇方法，例如：嶺回歸、橋回歸、lasso、SCAD、Adaptive lasso、Adaptive Robust lasso等。但這些研究大都針對響應(yīng)變量為數(shù)量型變量的情形，針對分類數(shù)據(jù)的情形目前研究甚少。因此，本文將把變量選擇問題的研究領(lǐng)域拓展到分類數(shù)據(jù)中，研究縱向數(shù)據(jù)下Logistic回歸模型的變量選擇問題。

半?yún)?shù)Logistic回歸模型既保持了非參數(shù)模型的靈活性與參數(shù)模型的簡潔性，又可以對分類數(shù)據(jù)進(jìn)行建模。通過采用等間距樣條函數(shù)逼近非參數(shù)部分，將非參數(shù)問題轉(zhuǎn)化為參數(shù)問題，給出非參數(shù)部分的全局?jǐn)M合值。在此基礎(chǔ)上，采用極小化懲罰似然的方法，得到常系數(shù)部分和非參數(shù)部分的估計。本文的懲罰函數(shù)主要使用lasso與SCAD，從不同角度刻畫懲罰方法的優(yōu)劣。

1 模型與估計方法

考慮二分類變量Y和p維協(xié)變量X的縱向統(tǒng)計建模問題。假設(shè)獲得如下縱向數(shù)據(jù)：

這里Yij和Xij=(Xij1,…,Xijp)T分別為第i個個體在第j個觀測時間Tij時的觀測值。通常假定mi有界。

記Xi=(Xi1,…,Ximi)T，Ti=(Ti1,…,Timi)T，其中Xij為p維列向量，Ti為mi維列向量，建立縱向數(shù)據(jù)的半?yún)?shù)Logistic回歸模型如下：

其中β=(β1,…,βp)T為p維未知參數(shù)向量，g(t)為未知光滑函數(shù)。并記：

其中時間變量Tij的取值范圍在0到T之間，μ(x)為Logistic函數(shù)，即稱為聯(lián)接函數(shù)，且滿足嚴(yán)格單增。當(dāng)連接函數(shù)退化為恒等映射時，Logistic半?yún)?shù)回歸模型即退化為普通的半?yún)?shù)模型。

為了給出模型中未知參數(shù)和未知函數(shù)的估計，取B樣條基函數(shù)(B1(t),…,BL(t))T，其中L表示樣條基的維數(shù)，L=J+r+1，J為節(jié)點個數(shù)，r為B樣條基的次數(shù)，在本文中取立方樣條，即r=3。因此g(t)≈γTB(t)，其中γ為L維列向量，γ=(γ1,…,γL)T。從而模型可以近似表示為：

其懲罰似然函數(shù)為：

這里，λs為正則化參數(shù)，其決定對系數(shù)壓縮強(qiáng)度的大小。l(β,γ)為對數(shù)似然函數(shù)，且：

式（6）中pλ(|·|)為懲罰函數(shù)，本文中應(yīng)用的懲罰函數(shù)有：①lasso懲罰pλ(|θ|)=λ|θ|，②SCAD懲罰：

通常a=3.7，對于正則化參數(shù)以及節(jié)點個數(shù)的選取，本文中采用K折交叉驗證的方法．將研究對象分成大致相等的K組，令k[i]為含有第i個個體的那一組。記與為去掉k[i]組后對β與γ的估計。則K折交叉驗證得分為：

令λs=λ,s=1,…,p，極小化上述CV得分選擇正則化參數(shù)λ與節(jié)點個數(shù)J。式（8）中，ωi為權(quán)重，結(jié)合縱向數(shù)據(jù)的特點，取ωi=1mi。式（8）給出的剔除個體的K折交叉驗證法與傳統(tǒng)的交叉驗證法的不同之處在于，同一個體的所有不同時間點的測量數(shù)據(jù)被同時刪除，這樣處理的好處在于不會破壞縱向數(shù)據(jù)的組內(nèi)相關(guān)性結(jié)構(gòu)。本文選擇固定節(jié)點位置的方式，當(dāng)然也可用CV方法通過數(shù)據(jù)驅(qū)動選擇節(jié)點位置，只是這樣處理將導(dǎo)致計算的復(fù)雜度極度上升，故建議在估計過程中采用Z.Huang等（2004）中所采用的等間距樣條方法，只選擇節(jié)點個數(shù)，從而大大降低計算量。

2 數(shù)據(jù)模擬

對于非參數(shù)的部分，給出兩種估計方法，即采用固定樣條節(jié)點個數(shù)，即節(jié)點個數(shù)約為n10和用CV方法選擇節(jié)點個數(shù)?？紤]縱向數(shù)據(jù)半?yún)?shù)Logistic模型：

其中β=(5,3,8,2,0,0,0,0,0,0)T，對于協(xié)變量的選取，分兩種情況進(jìn)行討論，以說明所提方法在不同情形下的數(shù)據(jù)表現(xiàn)。

情形1：協(xié)變量X1,…,X10相互獨立，來自于均勻分布U(-1,1)，此情況是為了檢驗均勻設(shè)計且獨立設(shè)計的情況下本文所提估計方法的數(shù)據(jù)表現(xiàn)。

情形2：協(xié)變量X1,…,X10相互獨立，來自于混合正態(tài)分布，其中混合有三個正態(tài)分布N(0,0.1)，N(-0.9,0.1)與N(0.9,0.1)。樣本來自這三個總體的概率分別為2 5,3 10,3 10，此情況是為了檢驗聚集設(shè)計的情況下估計量的數(shù)據(jù)表現(xiàn)。

針對不同的樣本量n=100，n=200，n=300均生成200個數(shù)據(jù)集，用以研究收斂速度。在數(shù)據(jù)模擬中，觀測時間是生成的，但存在隨機(jī)缺失。觀測時間生成方式為：每一個研究對象均有一個時間表{0,2,4,6,8,10,12}，并且每一個觀測時刻除時刻0之外均有20%的概率缺失。再對沒有缺失的時刻加上一個均勻分布U(-1,1)的隨機(jī)擾動，構(gòu)成最終的觀測時刻表。這樣生成的觀測時間，每個個體的觀測時刻均不相同，符合縱向數(shù)據(jù)的特點。

模型1：g(t)=sin(t)，在這種情況下，非參數(shù)的部分不能用線性模型近似，其隨時間呈非線性變化。

模型2：g(t)=tsin(t)12，相比模型1，此種模型非參數(shù)部分隨時間變化的特點不同于上一個模型，此模型隨時間的推移其變化愈發(fā)劇烈。

模型3：g(t)≡0，此模型為參數(shù)模型，本文中將應(yīng)用g(t)≡0的先驗信息進(jìn)行估計，并在不使用先驗信息的情況下采用本文的估計方法對參數(shù)部分進(jìn)行估計，用以檢驗當(dāng)沒有采用先驗信息（即模型設(shè)定并不精確）的條件下的估計方法的穩(wěn)健性。

模擬結(jié)果見表1至表6：

表1 模型1的變量選擇結(jié)果

表2 模型2的變量選擇結(jié)果

表3 模型3的變量選擇結(jié)果

表4 模型3有先驗信息(即已知g(t)≡0)時的變量選擇結(jié)果

表5 模型1采用CV方法選擇節(jié)點個數(shù)時的變量選擇結(jié)果(樣本量為100)

表6 模型1在樣本量為300時的表現(xiàn)

在表1至表6中，每個格內(nèi)的數(shù)據(jù)是200次模擬數(shù)據(jù)的估計均值，格中數(shù)據(jù)的下標(biāo)表示估計結(jié)果的標(biāo)準(zhǔn)差。從6個表中可以看到，對于SCAD估計方法，估計結(jié)果與真實結(jié)果相差很小，SCAD估計方法的偏差小于估計量方差的12，因此可以認(rèn)為SCAD方法是無偏的。SCAD方法的標(biāo)準(zhǔn)差約為非零參數(shù)的110，說明估計非常準(zhǔn)確，特別是表6當(dāng)中，對于非零參數(shù)部分，SCAD估計值與真實值之間的差值小于1%，與真實值已相當(dāng)接近。對于lasso方法，其偏度大于一個標(biāo)準(zhǔn)差，說明lasso方法有偏，這與普通均值模型中l(wèi)asso有偏而SCAD無偏的情況相同。觀察數(shù)據(jù)可以發(fā)現(xiàn)：對于有正則信息的最簡單情況（表4），不論是均勻設(shè)計還是聚集設(shè)計，在樣本量相同的情況下對于非零系數(shù)而言，SCAD懲罰的偏度要小于lasso，方差大于lasso。但SCAD較lasso的方差增大量小于其偏度的減小量，因此從均方誤差的角度講，SCAD懲罰的均方誤差小于lasso。對于零系數(shù)而言，兩種懲罰方法在偏度和方差上的表現(xiàn)相當(dāng)，SCAD方法的方差略大。當(dāng)樣本量相同時，不論lasso還是SCAD，均勻設(shè)計與聚集設(shè)計在偏度上的表現(xiàn)相當(dāng)，但聚集設(shè)計的方差較大。隨著樣本量的增大，兩種懲罰方法的偏差和方差均減小。表3中的數(shù)據(jù)表現(xiàn)與表4接近，表4中偏差和方差較表3來說，雖然可以看出一定程度的減小，但在采用SCAD時，在沒有正則信息的情況下，估計量的偏差非常小，小于真實參數(shù)的110，因此可以認(rèn)為當(dāng)g(t)≡0時，采用本文所述估計方法得到的估計效果與正確設(shè)定非參數(shù)部分形式時的估計效果相差并不太大，說明文中所提方法具有很好的穩(wěn)健性。

從表1與表6中還可以發(fā)現(xiàn)，不論是lasso方法還是SCAD方法，對于非零參數(shù)，聚集設(shè)計的方差與偏差均較大，因而聚集設(shè)計的均方誤差較大，但隨著樣本量的增大，兩種設(shè)計下的方差與偏差大小分別趨向于相同。隨著樣本量的增大，lasso方法與SCAD方法的標(biāo)準(zhǔn)差趨于相同，但SCAD方法的方差始終大于lasso。對于零參數(shù)，兩種方法的偏度相同，均小于標(biāo)準(zhǔn)差的110，因此可以認(rèn)為對于零參數(shù)來說，兩種方法均是無偏的。與非零參數(shù)相同的是SCAD的方差始終大于lasso，且隨著樣本量的增加，兩種懲罰方法的標(biāo)準(zhǔn)差趨于相同，均勻設(shè)計與聚集設(shè)計下的標(biāo)準(zhǔn)差趨于相同，但與非零參數(shù)的表現(xiàn)不同的是，零參數(shù)下均勻設(shè)計的標(biāo)準(zhǔn)差始終大于混合設(shè)計。故對非零參數(shù)而言，SCAD方法的均方誤差對不同的樣本量一致較小，對零參數(shù)一致較大。非零參數(shù)下的標(biāo)準(zhǔn)差是零參數(shù)下標(biāo)準(zhǔn)差的兩倍，說明本文所提兩種正則化方法確實起到了壓縮零參數(shù)的作用。比較表2自身的數(shù)據(jù)，表中估計量的表現(xiàn)與表1相同，說明不論非參數(shù)部分是怎樣的形式，文中所提估計量的數(shù)據(jù)表現(xiàn)均有相似的特點。

縱向比較表1、表2與表3、表4，發(fā)現(xiàn)對于lasso懲罰，不論理論模型是否為非參數(shù)模型，懲罰所導(dǎo)致的估計結(jié)果的偏差和方差大致不變，再次說明本文所提出的估計方法十分穩(wěn)健。

比較表1與表5發(fā)現(xiàn)，不論是否采用CV方法懲罰節(jié)點個數(shù)，估計效果大致相同，CV方法的偏差較大，但方差較小，CV方法的均方誤差小，但采用CV方法懲罰節(jié)點個數(shù)將急劇增加計算上的復(fù)雜度，為了節(jié)省估計時間，不建議采用交叉驗證法選擇節(jié)點個數(shù)，而是采用通常的方法，即節(jié)點個數(shù)約為樣本數(shù)量的十分之一。

此外，對于協(xié)變量存在多重共線性的情況，估計量的方差將變大。

3 結(jié)論

本文討論了縱向數(shù)據(jù)半?yún)?shù)Logistic模型估計，通過非參數(shù)函數(shù)的估計和參數(shù)部分的變量選擇方法，與數(shù)值模擬發(fā)現(xiàn)，所提方法在有限樣本下具有優(yōu)良的性質(zhì)。需要說明的是，在進(jìn)行變量選擇時，發(fā)現(xiàn)對于系數(shù)不為零的變量，不論協(xié)變量是聚集設(shè)計還是均勻設(shè)計，采用SCAD方法時，估計量的標(biāo)準(zhǔn)差以約為的速度縮小，而lasso方法的標(biāo)準(zhǔn)差縮小速度較慢。聚集設(shè)計的估計量方差普遍大于均勻設(shè)計。SCAD方法的偏差與均方誤小于lasso，方差普遍大于lasso，SCAD方法的偏差隨著樣本量的增加而減小。本文所提方法具有很好的穩(wěn)健性，當(dāng)理論模型為參數(shù)模型，而估計時設(shè)定為本文所提模型時，估計量的表現(xiàn)與正確設(shè)定模型形式時相當(dāng)。等間距樣條節(jié)點個數(shù)選為n10是合理的，與采用CV方法選擇節(jié)點個數(shù)時的表現(xiàn)相似，但CV方法造成計算上的負(fù)擔(dān)，故不建議用CV方法。

[1]P.J.Diggle,P.J.Heagerty,K.L.Liang,S.L.Zeger,Analysis of Longitudinal Data,Oxford University Press,Oxford,2002.

[2]Liya Fu,You-Gan Wang.Efficient Parameter Estimation via Gaussian Copulas for Quantile Regression with Longitudinal Data[J].Journal of Multivariate Analysis1(2016)43.

[3]Jianqing Fan and Runze Li.New Estimation and Model Selection Procedures for Semiparametric Modeling in Longitudinal Data Analysis[J].Journal of American Statistical Association,Vol.99,No.467(2004).

[4]Jianhua Z.Huang,Colin O.Wu and Lan Zhou.Polynomial Spline Estimation and Inference for Varying Coefficient Models with Longitudinal Data[J].Statistica Sinica(2004),14.

[5]Jianqing Fan,Tao Huang and Runze Li.Analysis of Longitudinal Data with Semiparametric Estimation of Covariance Function[J].Journal of American Statistical Association(2007)Vol.102.,No.478

[6]Hong Z,Hu Y and Lian H.Variable Selection for High-Dimensional Varying Coefficient Partially Linear Models via Nonconcave Penalty[J].Metrika,Volume 76,2013,7.

[7]Jianqing Fan and Runze Li.Variable Selection via Nonconcave Penalized Likelihood and its Oracle Properties[J].Journal of American Statistical Association,96,(2001)456.

[8]Hui Zou.The Adaptive Lasso and Its Oracle Properties[J].Journal of American Statistical Association,Vol.101,No.476,(2006).

[9]Verhasselt A.Generalized Varying Coefficient Models:A Smooth Variable Selection Technique[J].Statistica Sinica(2014),24.

[10]Jun Dong,Jason P.Estes,Gang Li and Damla Senturk.A Two-step Estimation Approach for Logistic Varying Coefficient Modeling of Longitudinal Data[J].Journal of Statistical Planning and Inference(2016),174.

Variable Selection for Semi-parametric Logistic Regression Model with Longitudinal Data

Gao Xianli1,Jiang Yuying2

(1.School of Statistics,Capital University of Economics and Business,Beijing 100070,China;2.Department of Basic Courses,Beijing Institute of Graphic Communication,Beijing 102600,China)

This paper studies the estimation problems of semi-parametric logistic regression model with longitudinal data,and presents a method of estimating the unknown parameter and the unknown function.The paper also investigates the variable selection of the parametric components,and makes a comparison among different methods of variable selection.Simulation results show that the method proposed in this paper has very good estimation effects.

longitudinal data;semi-parametric logistic regression model;variable selection.

O212.7

A

1002-6487（2017）20-0026-04

國家社會科學(xué)基金資助項目（10CTJ001）；首都經(jīng)濟(jì)貿(mào)易大學(xué)研究生科技創(chuàng)新重點項目；北京印刷學(xué)院校級資助項目（Eb201606）

高仙立（1991—），男，天津人，博士研究生，研究方向：應(yīng)用統(tǒng)計。

（通訊作者）姜玉英（1976—），女，山東濟(jì)寧人，副教授，研究方向：經(jīng)濟(jì)統(tǒng)計分析。

(責(zé)任編輯/亦 民)