趙明濤，何曉群

（1.中國(guó)人民大學(xué)統(tǒng)計(jì)學(xué)院，北京 100872；2.哈爾濱理工大學(xué)應(yīng)用科學(xué)學(xué)院，哈爾濱 150001）

0 引言

縱向數(shù)據(jù)[1,2,3]涉及到對(duì)不同個(gè)體的重復(fù)觀測(cè)數(shù)據(jù)，其獨(dú)特的結(jié)構(gòu)使得個(gè)體內(nèi)的觀測(cè)趨于相關(guān)，如何處理這種個(gè)體內(nèi)的相關(guān)性是研究縱向數(shù)據(jù)不可回避的問(wèn)題，這不僅是縱向分析的難點(diǎn)，也正是諸多縱向研究的出發(fā)點(diǎn)。邊際建模[1,2,3]方法是一種典型的縱向分析方法，其研究的重點(diǎn)為響應(yīng)變量與協(xié)變量的總體效應(yīng)，而把個(gè)體內(nèi)的相關(guān)性視為討厭參數(shù)。邊際模型的優(yōu)點(diǎn)在于分別對(duì)均值和方差進(jìn)行建模，只要均值模型假設(shè)正確，無(wú)論方差模型假定是否正確，總能得到均值部分的相合估計(jì)[1,2,3]。作為擬似然[4,5]方法在縱向研究中的一種推廣，廣義估計(jì)方程（GEE）方法由其諸多優(yōu)點(diǎn)已廣為大家熟知。然而，GEE估計(jì)在其估計(jì)效率以及解釋性等其他方面也存在一些缺點(diǎn)[6]。

針對(duì)GEE存在的這些問(wèn)題，二次推斷函數(shù)[7]（QIF）則彌補(bǔ)了GEE的諸多不足。QIF通過(guò)對(duì)工作相關(guān)矩陣的逆矩陣進(jìn)行基矩陣展開(kāi)近似，則避免了對(duì)討厭參數(shù)的估計(jì)；然后構(gòu)造擴(kuò)展得分函數(shù)，進(jìn)而利用廣義矩方法（GMM）構(gòu)造估計(jì)的目標(biāo)函數(shù)，則不僅提高了估計(jì)的效率，得到了比GEE估計(jì)更加穩(wěn)健的結(jié)果，而且使得估計(jì)結(jié)果總會(huì)存在且具有良好的大樣本性質(zhì)（相合性、漸近正態(tài)性）。同時(shí)，由于QIF本身是一個(gè)形式明確的推斷函數(shù)，故而可以直接運(yùn)用于關(guān)于模型的假設(shè)檢驗(yàn)[7]。

在數(shù)值實(shí)現(xiàn)方面，修正二次推斷函數(shù)[8]（MQIF）則克服了在一些情況下QIF中最優(yōu)權(quán)矩陣矩陣奇異的情形。MQIF通過(guò)一個(gè)線性收縮估計(jì)量替代最優(yōu)權(quán)矩陣，使得在任何情況下最優(yōu)權(quán)矩陣都是可逆的，而且提高了估計(jì)的精度。常用的Newton-Raphson數(shù)值迭代算法，由于需要求解雅可比矩陣的逆矩陣，不僅計(jì)算量太大，更有可能會(huì)出現(xiàn)雅可比矩陣奇異的情況，使得迭代不收斂或者得到的估計(jì)偏差太大。割線法基于向量分割方法，不需要計(jì)算雅可比矩陣的逆矩陣，不僅大大減少了計(jì)算量，而且避免了雅可比矩陣奇異的情形。對(duì)于縱向數(shù)據(jù)的非參數(shù)模型，本文選用回歸樣條[9]的方法，通過(guò)一組基對(duì)其展開(kāi)近似，進(jìn)而構(gòu)造基函數(shù)系數(shù)的MQIF。在數(shù)值實(shí)現(xiàn)方面本文則選用割線法迭代求解。

1 估計(jì)過(guò)程

假設(shè)縱向數(shù)據(jù)集合

滿足模型yij=f(tij)+εij，其中tij與yij分別為第i個(gè)個(gè)體的第j次觀測(cè)時(shí)刻點(diǎn)以及相應(yīng)的響應(yīng)變量，εij為零均值隨機(jī)誤差。f(?)為未知平滑函數(shù)。下面給出模型的低階矩的假設(shè)條件[4,5,6]：

其中ν(?)是已知的方差函數(shù)，進(jìn)一步假設(shè)不同個(gè)體的觀測(cè)之間相互獨(dú)立。

采用基函數(shù)展開(kāi)的方法近似未知平滑函數(shù)f(t)，此處不僅僅限于某一種基。不失一般性，選擇的一組基向量為：B(t)=[B1(t),…,Bk(t)]T，則可以得到f(t)=B(t)Tγ，其中γ為對(duì)應(yīng)的k×1維基函數(shù)系數(shù)向量。由此得到新的均值條件為

假設(shè)第i個(gè)個(gè)體內(nèi)觀測(cè)的協(xié)方差陣為，其中R(α)為工作相關(guān)矩陣，Ai為第i個(gè)個(gè)體觀測(cè)的邊際協(xié)方差陣，顯然Ai為一對(duì)角陣，進(jìn)而可以得到關(guān)于γ的GEE如下：

假設(shè)工作相關(guān)矩陣的逆矩陣R-1(α)有基矩陣展開(kāi)形式，將其代入（1），并構(gòu)造擴(kuò)展得分向量其中

為了避免出現(xiàn)C奇異的情況，根據(jù)GMM的思想，我們利用最優(yōu)權(quán)矩陣的一個(gè)可逆的線性收縮估計(jì)來(lái)替代C，進(jìn)而構(gòu)造關(guān)于基函數(shù)系數(shù)的γ的修正二次推斷函數(shù)（MQIF）。不妨假設(shè)

則得到：

且ei與ej相互獨(dú)立。不妨假設(shè)Ri=var(ei)，則對(duì)于任何固定設(shè)計(jì)，最優(yōu)權(quán)矩陣為，由此得到[14,15]E(C)=WN且

定義內(nèi)積＜H,K＞=其 中H∈?p×m，K∈?p×m，則矩陣范數(shù)為則得到WN的一個(gè)線性收縮估計(jì)為其中為單位陣。顯然，一定條件下，在二次期望損失意義下SN為WN的漸近最優(yōu)的相合估計(jì)。由于SN中的系數(shù)是未知的，故給出SN的估計(jì)量為其中

對(duì)于固定設(shè)計(jì)，在滿足一定條件下，有

由此構(gòu)造基函數(shù)系數(shù)γ的修正二次推斷函數(shù)為

得到γ的修正二次推斷函數(shù)估計(jì)為

顯然所構(gòu)造的修正二次推斷函數(shù)QN(γ)則避免了二次推斷函數(shù)中容易出現(xiàn)的權(quán)矩陣奇異的情形，使得權(quán)矩陣總是可逆的。

實(shí)際中通常利用Newton-Raphson迭代求解，然而在復(fù)雜數(shù)據(jù)情形下，由于需要計(jì)算雅可比矩陣的逆，使得迭代算法的計(jì)算量很大難于實(shí)現(xiàn)。不僅如此，在某些特定情況下也可能會(huì)出現(xiàn)雅可比矩陣奇異的情形，使得Newton-Raphson算法無(wú)法實(shí)現(xiàn)。為了避免這些問(wèn)題，我們利用割線法迭代求出γ的數(shù)值估計(jì)。割線法的迭代公式[16]為

第一步：給出γ的初始值以及收斂閥值ε0；

否則返回第二步迭代計(jì)算，直到收斂為止；

2 估計(jì)的漸近性質(zhì)

本節(jié)我們將研究上述估計(jì)過(guò)程的大樣本性質(zhì)，再給出大樣本性質(zhì)前，我們先給出一些必要的假設(shè)條件[13,14]。

條件1：最優(yōu)權(quán)矩陣WN→W0,a.s.，其中W0為一個(gè)可逆的常矩陣；

條件2：基函數(shù)系數(shù)γ是可識(shí)別的，也即存在唯一的γ0∈S，滿足均值模型假設(shè)E(gi(γ0))=0，其中S為基函數(shù)系數(shù)空間；

條件3：基函數(shù)系數(shù)空間S為一個(gè)緊空間，并且γ0為S的一個(gè)內(nèi)點(diǎn)；

條件4：E(gi(γ))關(guān)于γ連續(xù)；

條件5：(γ)關(guān)于γ的導(dǎo)數(shù)存在連續(xù)，且

下面給出基函數(shù)系數(shù)估計(jì)?的漸近性質(zhì)：

定理1：若條件（1～4）滿足，則通過(guò)（4）得到的γ?存在且

證明：顯然修正二次推斷函數(shù)有界，因此基函數(shù)系數(shù)γ的估計(jì)存在。不妨設(shè)U為γ0的任意一個(gè)鄰域。需證?不可能落在U的外邊。

由大數(shù)定律和條件（1）可知：

對(duì)任意一個(gè)緊集，有如下一致大數(shù)定律：

由此可以得到

定理2：若條件（1～5）均滿足，則通過(guò)（4）得到的基函數(shù)系數(shù)估計(jì)是漸近正態(tài)的

證明：由于?是由（4）得到的，所以可得到根據(jù)泰勒展開(kāi)得到：

其中γ?位于γ?和γ0之間。因此有

下面給出定理2的兩個(gè)推論：

推論1：基函數(shù)系數(shù)γ的修正二次推斷函數(shù)估計(jì)γ?的協(xié)方差矩陣的相合估計(jì)為

推論2：假設(shè)γ?的均值和協(xié)方差都存在，在定理2的條件都滿足的情況的估計(jì)為其協(xié)方差矩陣的估計(jì)為且f(t) 的100(1-α)%的置信區(qū)間為：

3 數(shù)值模擬

考慮模型[13]其中i=1,…,200。每個(gè)個(gè)體被觀測(cè)的時(shí)刻點(diǎn)集為{0,1,…,30}，每個(gè)真實(shí)的觀測(cè)時(shí)刻點(diǎn)與規(guī)定的時(shí)刻點(diǎn)之間的誤差服從均勻分布U(-0.5,0.5)，并且除了起始時(shí)刻點(diǎn)外，其他時(shí)刻點(diǎn)的觀測(cè)缺失的概率為0.6。假設(shè)個(gè)體內(nèi)的相關(guān)結(jié)構(gòu)為可交換結(jié)構(gòu)且ρ=0.8，εij～N(0,2)。對(duì)待上述兩個(gè)模型分別選用基B(t)=[1,t,t2,t3]T與B(t)=[1,t]T，利用（4）求解基函數(shù)系數(shù)γ的MQIF估計(jì)γ?，重復(fù)多遍，求均值，進(jìn)而得到擬合值。

圖1

圖2

圖1與圖2分別給出了真實(shí)的平滑函數(shù)y=15+20sin的曲線及其擬合曲線，其中黑色實(shí)線為真實(shí)曲線，紅色實(shí)線為擬合曲線。從圖上看出，在觀測(cè)時(shí)段內(nèi)，兩個(gè)真實(shí)的平滑函數(shù)與其MQIF估計(jì)的擬合曲線幾乎重合，擬合效果很好。由此說(shuō)明了文中所提方法的有效性以及實(shí)用性。

[1]王啟華，史寧中，耿直.現(xiàn)代統(tǒng)計(jì)研究基礎(chǔ)[M].北京：科學(xué)出版社,2010.

[2]Diggle P J,Liang,K Y ，Zeger.Analysis of Longitudinal Data[M].Lon?don：Oxford University Press，1994.

[3]Hedeker D,Gibbons D.Longitudinal Data Analysis[M].New York:John Wiley&Sons,2006.

[4]Wedderburn R W M.Quasi-likelihood Functions,Generalized Linear Models,and the Guass-Newton Method[J].Biometrika,1974,61(3).

[5]McCullagh P.Quasi-likelihood functions.The annals of statistics,1983,11(1).

[6]Crowder M.On the Use of a Working Correlation Matrix in Using Gen?eralized Linear Models for Repeated Measures[J].Biometrika,1995,(82).

[7]Qu A,Lindsay B G,and Li B.Improving Generalized Estimating Equations Using Quadratic Inference Functions[J].Biometrika,2000,(87).

[8]Han P S,Song Peter X.-K.A Note on Improving Quadratic Inference Functions Using a Linear Shrinkage Approach[J].Statistics and proba?bility letters,2011,(81).

[9]Wu H,Zhang J.Nonparametric Regression Methods for Longitudinal Data Analysis[M].New York：John Wiley&Sons,2006.