田瑞琴，徐登可

(1. 杭州師范大學(xué)理學(xué)院，浙江 杭州 311121； 2. 浙江農(nóng)林大學(xué)理學(xué)院，浙江 杭州 311300)

0 引言

縱向數(shù)據(jù)廣泛地產(chǎn)生于經(jīng)濟(jì)學(xué)、醫(yī)學(xué)、生物學(xué)和社會(huì)學(xué)等領(lǐng)域，它是指同一個(gè)體在不同時(shí)間點(diǎn)上觀測獲得的重復(fù)觀測數(shù)據(jù).由于數(shù)據(jù)結(jié)構(gòu)中的重復(fù)觀測性質(zhì)，容易導(dǎo)致其產(chǎn)生數(shù)據(jù)的缺失.例如，有些研究個(gè)體會(huì)在某次實(shí)驗(yàn)缺席或者在實(shí)驗(yàn)結(jié)束之前退出研究.最簡單的處理缺失數(shù)據(jù)的思想是直接對完全觀測數(shù)據(jù)進(jìn)行分析.但是，當(dāng)缺失數(shù)據(jù)的機(jī)制不是完全隨機(jī)缺失時(shí)，所得的統(tǒng)計(jì)推斷結(jié)果會(huì)出現(xiàn)較大的偏差.近年來對縱向缺失數(shù)據(jù)的研究發(fā)現(xiàn)，單調(diào)缺失和非單調(diào)缺失是兩種主要的缺失模式，其中單調(diào)缺失主要是指在實(shí)驗(yàn)中部分研究個(gè)體因?yàn)槟撤N原因在某一個(gè)時(shí)刻退出該實(shí)驗(yàn)后再也沒回來繼續(xù)實(shí)驗(yàn)，這樣獲得的缺失數(shù)據(jù)稱為單調(diào)缺失數(shù)據(jù)，否則稱之為非單調(diào)缺失數(shù)據(jù).例如，文[1]研究老年癡呆癥疾病，由于病人發(fā)展成癡呆患者而放棄了該項(xiàng)研究，所獲得的數(shù)據(jù)即為單調(diào)缺失數(shù)據(jù).目前針對縱向缺失數(shù)據(jù)模型的參數(shù)估計(jì)問題已有不少研究，但是大多集中在協(xié)變量維數(shù)是固定的情況下.Robins等[2]針對非單調(diào)缺失機(jī)制下縱向數(shù)據(jù)半?yún)?shù)模型，基于逆概率加權(quán)廣義估計(jì)方程(IPW-GEE)研究分析了模型的參數(shù)估計(jì)問題.Qu等[3]基于無偏估計(jì)函數(shù)方法對具有隨機(jī)缺失的相關(guān)數(shù)據(jù)模型提出了一種有效的參數(shù)估計(jì)方法.文[4-7]也對此類相關(guān)問題展開了研究.

另外，很多學(xué)者針對協(xié)變量是發(fā)散維的高維縱向數(shù)據(jù)回歸模型進(jìn)行了變量選擇研究.例如，Xu等[8]針對協(xié)變量發(fā)散維的相關(guān)二元數(shù)據(jù)模型，提出了一種基于懲罰加權(quán)最小二乘的變量選擇方法.Dziak[9]針對高維縱向數(shù)據(jù)模型提出了基于SCAD懲罰二次推斷函數(shù)的變量選擇方法.Yang等[10]針對高維部分線性模型提出了Dantzig變量選擇方法，最后能同時(shí)獲得參數(shù)部分和非參數(shù)部分的估計(jì).Wang[11]針對協(xié)變量是發(fā)散維的高維縱向數(shù)據(jù)模型，拓展了廣義估計(jì)方程分析方法的漸近理論.Tian等[12]針對高維數(shù)據(jù)下縱向部分線性回歸模型，基于光滑閾廣義估計(jì)方程方法研究分析了模型的變量選擇.綜上所述，目前針對帶有缺失數(shù)據(jù)或高維數(shù)據(jù)的縱向數(shù)據(jù)回歸模型的研究成果已經(jīng)比較豐富，但甚少涉及缺失數(shù)據(jù)下高維縱向回歸模型的變量選擇問題.現(xiàn)有的變量選擇方法大部分是基于懲罰函數(shù)(例如SCAD懲罰、Adaptive LASSO懲罰等)進(jìn)行壓縮估計(jì).由于懲罰函數(shù)在0點(diǎn)是奇異的，所以基于懲罰估計(jì)的變量選擇算法就要面臨去解決一個(gè)涉及非凸函數(shù)最優(yōu)化的問題.為了解決這個(gè)問題，Ueki[13]提出了一種基于光滑閾估計(jì)方程(SEE)的變量選擇方法.這種方法能自動(dòng)地剔除不重要的變量，同時(shí)得到重要變量的參數(shù)估計(jì)，并且此變量選擇算法中不涉及凸優(yōu)化的問題.Lai等[14]將這種光滑閾估計(jì)方程思想應(yīng)用到縱向數(shù)據(jù)單指標(biāo)模型上.Li等[15]和Tian等[16]基于光滑閾估計(jì)方程思想并結(jié)合廣義估計(jì)方程分別研究分析了縱向數(shù)據(jù)下廣義線性模型和縱向數(shù)據(jù)下變系數(shù)部分線性模型的變量選擇.

本文針對單調(diào)缺失數(shù)據(jù)下縱向高維部分線性回歸模型，結(jié)合逆概率加權(quán)廣義估計(jì)方程和Ueki[13]的變量選擇思想，提出了一種基于逆概率加權(quán)光滑閾廣義估計(jì)方程(IPW-SGEE)的變量選擇方法，其中非參數(shù)部分使用樣條方法逼近.該方法具有下列3個(gè)特點(diǎn)：第一，因?yàn)椴簧婕巴购瘮?shù)最優(yōu)化問題，所以算法實(shí)現(xiàn)相對簡便，且獲得估計(jì)具有Oracle性質(zhì).第二，將Ueki提出的光滑閾估計(jì)方程方法[13]進(jìn)行了推廣，使其應(yīng)用到協(xié)變量是發(fā)散維且響應(yīng)變量是單調(diào)缺失的縱向數(shù)據(jù)半?yún)?shù)回歸模型中.第三，與Wang[11]方法比較，本文方法分析了半?yún)?shù)回歸模型.半?yún)?shù)回歸模型既具有參數(shù)模型的優(yōu)點(diǎn)，也具有非參數(shù)回歸模型的優(yōu)點(diǎn)，因此有更廣泛的實(shí)際應(yīng)用價(jià)值.

本文首先介紹模型和缺失機(jī)制，以及基于逆概率加權(quán)光滑閾估計(jì)方程的變量選擇方法.然后在一些正則條件下研究給出該IPW-SGEE變量選擇方法的漸近理論性質(zhì)，并詳細(xì)列出獲得IPW-SGEE估計(jì)的迭代計(jì)算算法.最后通過隨機(jī)模擬研究驗(yàn)證了該方法的有限樣本性質(zhì).

1 模型和缺失機(jī)制

對于高維縱向數(shù)據(jù)，考慮如下部分線性模型：

(1)

這里考慮響應(yīng)變量Yij缺失，當(dāng)Rij=1，Yij可觀測；Rij=0，Yij缺失.進(jìn)一步考慮Yij是單調(diào)隨機(jī)缺失，即當(dāng)Rij=0時(shí)，意味著Rik=0且k≥j.另外假定，對于所有個(gè)體初始時(shí)刻均可以觀測，即Ri1=1(i=1,…,n).對于缺失數(shù)據(jù)問題一般考慮隨機(jī)缺失機(jī)制，即假定

(2)

2 IPW-SGEE方法

類似于He等[17]，利用樣條回歸逼近非參數(shù)分量f(·).具體地，設(shè)B(u)=(B1(u),…,BL(u))T是階數(shù)為M的B-樣條基函數(shù)，其中L=K+M，K為內(nèi)節(jié)點(diǎn)個(gè)數(shù).B-樣條基具有有界支撐且計(jì)算穩(wěn)定[18].節(jié)點(diǎn)選擇是樣條光滑方法中很重要的一個(gè)部分，按照文[17]，本文節(jié)點(diǎn)的個(gè)數(shù)取為N1/5的整數(shù)部分.f(Tij)可以由B(Tij)Tθ逼近，其中θ∈RL是樣條回歸系數(shù)向量.由此，回歸模型 (1)可表示為：

(3)

(WTA-1W)-1WTA-1(Y-Xβn),

在單調(diào)缺失框架下，按照Robins等[2]的思想，采用逆概率加權(quán)廣義估計(jì)方程(IPW-GEE)方法.具體地，關(guān)于βn的IPW-GEE 函數(shù)為：

其中Πi=diag{Rij/πij(α),j=1,…,mi}，用于校正由缺失數(shù)據(jù)造成的偏差.然而，本文的主要目標(biāo)是同時(shí)估計(jì)和選擇重要的協(xié)變量，并且為了避免懲罰估計(jì)所帶來的非凸函數(shù)最優(yōu)化的問題，基于SEE思想，關(guān)于參數(shù)βn構(gòu)造如下的逆概率加權(quán)光滑閾廣義估計(jì)方程(IPW-SGEE)：

(Ipn-Δ)Un(βn,α,ρ)+Δβn=0,

(4)

其中Δ是對角線元素為δ=(δj)j=1,…,pn的對角矩陣，Ipn是pn維的單位陣.注意到，式(4)中第j(j=1,…,pn)個(gè)方程，如果δj=1，那么有βnj=0.這就說明式(4)可以得到稀疏解.

(5)

(6)

3 漸近性質(zhì)

1)參數(shù)向量α是緊集空間Γ上的內(nèi)點(diǎn).

2)非參數(shù)函數(shù)f(·)是r階可導(dǎo)，其中r≥2.

4)存在正常數(shù)c1和c2使得

其中λmin(或者λmax)表示矩陣的最小(或者最大)特征值.

5)對所有的pn，存在正常數(shù)c3和c4使得

其中λmin(或者λmax)表示矩陣的最小(或者最大)特征值.

6)P(Rij=1|Yi,Xi,Ti) 有界遠(yuǎn)離0.

ii)漸近正態(tài)性，即

注2定理1和定理2的證明可參見文獻(xiàn)[12]和[23].

4 迭代算法

4.1 計(jì)算算法

為得到式(6)中參數(shù)βn的IPW-SGEE估計(jì)，需要得到相關(guān)參數(shù)α和ρ的相合估計(jì).首先，通過極大似然估計(jì)方法得到缺失模型參數(shù)α的估計(jì).然后，使用迭代算法估計(jì)βn和ρ.按照文[19]，利用校正矩方法估計(jì)ρ.定義如下皮爾遜殘差

(7)

當(dāng)相關(guān)結(jié)構(gòu)是AR(1)時(shí)，工作矩陣R(ρ)的第(t,t+1)個(gè)元素的估計(jì)為

進(jìn)而，迭代計(jì)算具體步驟如下.

步驟4通過式(8)獲得βn的更新估計(jì)：

(8)

步驟5重復(fù)步驟3～4直到滿足收斂條件，并且記最終的估計(jì)為IPW-SGEE估計(jì).

4.2 調(diào)整參數(shù)的選擇

使用BIC準(zhǔn)則[20]選擇調(diào)整參數(shù)(λ,γ)，其中具體的BIC準(zhǔn)則定義如下：

(9)

5 模擬研究

以下主要通過隨機(jī)模擬實(shí)驗(yàn)來驗(yàn)證上文所述逆概率加權(quán)光滑閾廣義估計(jì)方程方法的有限樣本性質(zhì).產(chǎn)生數(shù)據(jù)的模型如下：

(10)

響應(yīng)變量Yij單調(diào)隨機(jī)缺失，并且假定每個(gè)個(gè)體在初始時(shí)刻都是可觀測的.缺失數(shù)據(jù)模型為logit(λij)=α0+α1Yi,j-1,其中λij=P(Rij=1|Ri,j-1=1,Yi,j-1)，參數(shù)真值α1=-3，α0的真值分別取5和1，相應(yīng)的缺失比例分別約為26%和40%.

表1 真實(shí)相關(guān)結(jié)構(gòu)為CS時(shí)高維部分線性模型 (10)的變量選擇結(jié)果Tab.1 Variable selections for high-dimensional partially linear model (10) when the true correlation structure is exchangeable

續(xù)表1

表 2 真實(shí)相關(guān)結(jié)構(gòu)為AR(1)時(shí)高維部分線性模型(10)的變量選擇結(jié)果Tab.2 Variable selections for high-dimensional partially linear model (10) when the true correlation structure is AR(1)

續(xù)表 2

由表1和表2可得：1)即使在錯(cuò)誤指定工作相關(guān)矩陣下，隨著ρ的增大，IPW-SGEE方法的變量選擇效果也越來越好.2)隨著缺失比例變小，IPW-SGEE、IPWSCAD和IPWLasso這3種變量選擇效果也變得越好.另外，與IPWSCAD和IPWLasso方法相比，IPW-SGEE變量選擇方法在各個(gè)指標(biāo)方面的表現(xiàn)效果都要較好一些.3)在正確指定相關(guān)結(jié)構(gòu)下，IPW-SGEE方法的變量選擇效果稍好于錯(cuò)誤指定相關(guān)結(jié)構(gòu)下的效果.這也說明本文所提出的方法并不顯著地依賴工作相關(guān)結(jié)構(gòu).