赵明涛，何晓群

（1.中国人民大学统计学院，北京 100872；2.哈尔滨理工大学应用科学学院，哈尔滨 150001）

0 引言

纵向数据[1,2,3]涉及到对不同个体的重复观测数据，其独特的结构使得个体内的观测趋于相关，如何处理这种个体内的相关性是研究纵向数据不可回避的问题，这不仅是纵向分析的难点，也正是诸多纵向研究的出发点。边际建模[1,2,3]方法是一种典型的纵向分析方法，其研究的重点为响应变量与协变量的总体效应，而把个体内的相关性视为讨厌参数。边际模型的优点在于分别对均值和方差进行建模，只要均值模型假设正确，无论方差模型假定是否正确，总能得到均值部分的相合估计[1,2,3]。作为拟似然[4,5]方法在纵向研究中的一种推广，广义估计方程（GEE）方法由其诸多优点已广为大家熟知。然而，GEE估计在其估计效率以及解释性等其他方面也存在一些缺点[6]。

针对GEE存在的这些问题，二次推断函数[7]（QIF）则弥补了GEE的诸多不足。QIF通过对工作相关矩阵的逆矩阵进行基矩阵展开近似，则避免了对讨厌参数的估计；然后构造扩展得分函数，进而利用广义矩方法（GMM）构造估计的目标函数，则不仅提高了估计的效率，得到了比GEE估计更加稳健的结果，而且使得估计结果总会存在且具有良好的大样本性质（相合性、渐近正态性）。同时，由于QIF本身是一个形式明确的推断函数，故而可以直接运用于关于模型的假设检验[7]。

在数值实现方面，修正二次推断函数[8]（MQIF）则克服了在一些情况下QIF中最优权矩阵矩阵奇异的情形。MQIF通过一个线性收缩估计量替代最优权矩阵，使得在任何情况下最优权矩阵都是可逆的，而且提高了估计的精度。常用的Newton-Raphson数值迭代算法，由于需要求解雅可比矩阵的逆矩阵，不仅计算量太大，更有可能会出现雅可比矩阵奇异的情况，使得迭代不收敛或者得到的估计偏差太大。割线法基于向量分割方法，不需要计算雅可比矩阵的逆矩阵，不仅大大减少了计算量，而且避免了雅可比矩阵奇异的情形。对于纵向数据的非参数模型，本文选用回归样条[9]的方法，通过一组基对其展开近似，进而构造基函数系数的MQIF。在数值实现方面本文则选用割线法迭代求解。

1 估计过程

假设纵向数据集合

满足模型yij=f(tij)+εij，其中tij与yij分别为第i个个体的第j次观测时刻点以及相应的响应变量，εij为零均值随机误差。f(⋅)为未知平滑函数。下面给出模型的低阶矩的假设条件[4,5,6]：

其中ν(⋅)是已知的方差函数，进一步假设不同个体的观测之间相互独立。

采用基函数展开的方法近似未知平滑函数f(t)，此处不仅仅限于某一种基。不失一般性，选择的一组基向量为：B(t)=[B1(t),…,Bk(t)]T，则可以得到f(t)=B(t)Tγ，其中γ为对应的k×1维基函数系数向量。由此得到新的均值条件为

假设第i个个体内观测的协方差阵为，其中R(α)为工作相关矩阵，Ai为第i个个体观测的边际协方差阵，显然Ai为一对角阵，进而可以得到关于γ的GEE如下：

假设工作相关矩阵的逆矩阵R-1(α)有基矩阵展开形式，将其代入（1），并构造扩展得分向量其中

为了避免出现C奇异的情况，根据GMM的思想，我们利用最优权矩阵的一个可逆的线性收缩估计来替代C，进而构造关于基函数系数的γ的修正二次推断函数（MQIF）。不妨假设

则得到：

且ei与ej相互独立。不妨假设Ri=var(ei)，则对于任何固定设计，最优权矩阵为，由此得到[14,15]E(C)=WN且

定义内积＜H,K＞=其 中H∈ℝp×m，K∈ℝp×m，则矩阵范数为则得到WN的一个线性收缩估计为其中为单位阵。显然，一定条件下，在二次期望损失意义下SN为WN的渐近最优的相合估计。由于SN中的系数是未知的，故给出SN的估计量为其中

对于固定设计，在满足一定条件下，有

由此构造基函数系数γ的修正二次推断函数为

得到γ的修正二次推断函数估计为

显然所构造的修正二次推断函数QN(γ)则避免了二次推断函数中容易出现的权矩阵奇异的情形，使得权矩阵总是可逆的。

实际中通常利用Newton-Raphson迭代求解，然而在复杂数据情形下，由于需要计算雅可比矩阵的逆，使得迭代算法的计算量很大难于实现。不仅如此，在某些特定情况下也可能会出现雅可比矩阵奇异的情形，使得Newton-Raphson算法无法实现。为了避免这些问题，我们利用割线法迭代求出γ的数值估计。割线法的迭代公式[16]为

第一步：给出γ的初始值以及收敛阀值ε0；

否则返回第二步迭代计算，直到收敛为止；

2 估计的渐近性质

本节我们将研究上述估计过程的大样本性质，再给出大样本性质前，我们先给出一些必要的假设条件[13,14]。

条件1：最优权矩阵WN→W0,a.s.，其中W0为一个可逆的常矩阵；

条件2：基函数系数γ是可识别的，也即存在唯一的γ0∈S，满足均值模型假设E(gi(γ0))=0，其中S为基函数系数空间；

条件3：基函数系数空间S为一个紧空间，并且γ0为S的一个内点；

条件4：E(gi(γ))关于γ连续；

条件5：(γ)关于γ的导数存在连续，且

下面给出基函数系数估计̂的渐近性质：

定理1：若条件（1～4）满足，则通过（4）得到的γ̂存在且

证明：显然修正二次推断函数有界，因此基函数系数γ的估计存在。不妨设U为γ0的任意一个邻域。需证̂不可能落在U的外边。

由大数定律和条件（1）可知：

对任意一个紧集，有如下一致大数定律：

由此可以得到

定理2：若条件（1～5）均满足，则通过（4）得到的基函数系数估计是渐近正态的

证明：由于̂是由（4）得到的，所以可得到根据泰勒展开得到：

其中γ͂位于γ̂和γ0之间。因此有

下面给出定理2的两个推论：

推论1：基函数系数γ的修正二次推断函数估计γ̂的协方差矩阵的相合估计为

推论2：假设γ̂的均值和协方差都存在，在定理2的条件都满足的情况的估计为其协方差矩阵的估计为且f(t) 的100(1-α)%的置信区间为：

3 数值模拟

考虑模型[13]其中i=1,…,200。每个个体被观测的时刻点集为{0,1,…,30}，每个真实的观测时刻点与规定的时刻点之间的误差服从均匀分布U(-0.5,0.5)，并且除了起始时刻点外，其他时刻点的观测缺失的概率为0.6。假设个体内的相关结构为可交换结构且ρ=0.8，εij～N(0,2)。对待上述两个模型分别选用基B(t)=[1,t,t2,t3]T与B(t)=[1,t]T，利用（4）求解基函数系数γ的MQIF估计γ̂，重复多遍，求均值，进而得到拟合值。

图1

图2

图1与图2分别给出了真实的平滑函数y=15+20sin的曲线及其拟合曲线，其中黑色实线为真实曲线，红色实线为拟合曲线。从图上看出，在观测时段内，两个真实的平滑函数与其MQIF估计的拟合曲线几乎重合，拟合效果很好。由此说明了文中所提方法的有效性以及实用性。

[1]王启华，史宁中，耿直.现代统计研究基础[M].北京：科学出版社,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.