田瑞琴，徐登可

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

0 引言

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

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

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

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

1 模型和缺失機制

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

(1)

這里考慮響應(yīng)變量Yij缺失，當(dāng)Rij=1，Yij可觀測；Rij=0，Yij缺失.進一步考慮Yij是單調(diào)隨機缺失，即當(dāng)Rij=0時，意味著Rik=0且k≥j.另外假定，對于所有個體初始時刻均可以觀測，即Ri1=1(i=1,…,n).對于缺失數(shù)據(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é)點個數(shù).B-樣條基具有有界支撐且計算穩(wěn)定[18].節(jié)點選擇是樣條光滑方法中很重要的一個部分，按照文[17]，本文節(jié)點的個數(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)廣義估計方程(IPW-GEE)方法.具體地，關(guān)于βn的IPW-GEE 函數(shù)為：

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

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

(4)

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

(5)

(6)

3 漸近性質(zhì)

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

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) 有界遠離0.

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

注2定理1和定理2的證明可參見文獻[12]和[23].

4 迭代算法

4.1 計算算法

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

(7)

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

進而，迭代計算具體步驟如下.

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

(8)

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

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

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

(9)

5 模擬研究

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

(10)

響應(yīng)變量Yij單調(diào)隨機缺失，并且假定每個個體在初始時刻都是可觀測的.缺失數(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 真實相關(guān)結(jié)構(gòu)為CS時高維部分線性模型 (10)的變量選擇結(jié)果Tab.1 Variable selections for high-dimensional partially linear model (10) when the true correlation structure is exchangeable

續(xù)表1

表 2 真實相關(guān)結(jié)構(gòu)為AR(1)時高維部分線性模型(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)即使在錯誤指定工作相關(guān)矩陣下，隨著ρ的增大，IPW-SGEE方法的變量選擇效果也越來越好.2)隨著缺失比例變小，IPW-SGEE、IPWSCAD和IPWLasso這3種變量選擇效果也變得越好.另外，與IPWSCAD和IPWLasso方法相比，IPW-SGEE變量選擇方法在各個指標(biāo)方面的表現(xiàn)效果都要較好一些.3)在正確指定相關(guān)結(jié)構(gòu)下，IPW-SGEE方法的變量選擇效果稍好于錯誤指定相關(guān)結(jié)構(gòu)下的效果.這也說明本文所提出的方法并不顯著地依賴工作相關(guān)結(jié)構(gòu).