文章快速检索     高级检索
  北京化工大学学报(自然科学版)  2018, Vol. 45 Issue (2): 100-104   DOI: 10.13543/j.bhxbzr.2018.02.017
0

引用本文  

关晓妮, 黄彬. 纵向数据下广义线性模型的稳健二次推断函数估计[J]. 北京化工大学学报(自然科学版), 2018, 45(2): 100-104. DOI: 10.13543/j.bhxbzr.2018.02.017.
GUAN XiaoNi, HUANG Bin. Robust quadratic inference functions for generalized linear models with longitudinal data[J]. Journal of Beijing University of Chemical Technology (Natural Science), 2018, 45(2): 100-104. DOI: 10.13543/j.bhxbzr.2018.02.017.

基金项目

国家自然科学基金(11471321)

第一作者

关晓妮, 女, 1990年生, 硕士生.

通信联系人

黄彬, E-mail:abinhuang@gmail.com

文章历史

收稿日期:2017-08-01
纵向数据下广义线性模型的稳健二次推断函数估计
关晓妮 , 黄彬     
北京化工大学 理学院, 北京 100029
摘要:针对纵向数据下的广义线性模型,为了有效控制离群点对估计的影响以及进一步提高估计的效率,利用二次推断函数(QIF)改进加权的指数得分函数,得到了模型参数有效且稳健的二次推断函数估计(ERQIF),并证明了在一定条件下所得估计的相合性和渐近正态性。数值计算结果进一步表明,当离群点存在或工作相关矩阵被错误指定时,所得估计有稳健的模拟结果。
关键词纵向数据    加权指数得分函数    二次推断函数法(QIF)    稳健估计    
Robust quadratic inference functions for generalized linear models with longitudinal data
GUAN XiaoNi , HUANG Bin     
Faculty of Science, Beijing University of Chemical Technology, Beijing 100029, China
Abstract: This paper presents an efficient and robust estimation method for generalized linear models with longitudinal data. By using a quadratic inferential function (QIF) to improve the weighted exponential score function, we can obtain an effective and robust quadratic inferential function(ERQIF). Under some regularity conditions, the resulting estimators are consistent and asymptotically normal distributed. Finally, simulation studies show that the proposed estimators have robust and efficient numerical results, even when many outliers are included and the working correlation matrix is misspecified.
Key words: longitudinal data    weighted exponential score function    quadratic inferential function(QIF)    robust estimation    
引言

纵向数据广泛应用在生物医药、流行病学和经济学等领域中,它通过对同一个体在不同的时间点重复测量得到, 具有组内相关、组间独立的特性。在实际的统计建模中,如果忽略这种特性,将会降低估计的效率,所以如何处理组内相关性是纵向数据研究不可避免的问题。Liang等[1]提出的广义估计方程(GEE)方法可以有效地处理组内相关问题,该方法利用含少量讨厌参数的工作相关矩阵建立估计方程,即使工作矩阵错误也可以得到回归系数的相合估计。但是GEE方法首先需要得到讨厌参数的相合估计,同时讨厌参数的相合估计可能不存在。为了提高估计的效率,Qu等[2]提出了二次推断函数法(QIF),利用某些基矩阵的线性组合来逼近工作相关矩阵的逆矩阵和避免估计讨厌参数,既使所得估计具有相合性,又提高了参数估计的效率。已有许多学者将GEE和QIF方法应用于不同的模型[3-5],但是这些应用方法大多类似于加权最小二乘估计,不具有稳健性。

在纵向数据中,由于重复测量,一个个体的突变往往会产生一系列的离群点,因此在纵向数据的研究中稳健估计是一个十分重要的问题[6-8]。为了获得更好的稳健性和提高估计的效率,Wang等[9]提出了基于指数得分函数的稳健回归估计方法, 通过引入一个调节参数来提高估计的稳健性和有效性。Lv等[10]将加权指数得分函数和GEE方法相结合,得到了纵向数据下广义线性模型参数的有效且稳健的估计(ERGEE)。但是估计讨厌参数必然会降低估计的效率。本文将指数得分函数和QIF方法相结合, 利用加权的指数得分函数有效控制协变量和因变量中离群点的影响,同时基于QIF方法建立稳健二次推断函数,得到了纵向数据下广义线性模型参数的ERQIF估计,并进一步提高了估计的效率。最后通过模拟计算进一步验证了所得估计的有限样本性质。

1 参数估计方法 1.1 ERQIF方法

假定纵向数据中包含n个个体,对每个个体重复观测mi次, 总的观测次数$M = \sum\limits_{i = 1}^n {{m_i}} $。设Yi=(yi1, …, yimi)T为第i个个体的响应变量(i=1, …, n),Xi=(Xi1, Xi2, …, Ximi)T为其相应的协变量,其中Xij=(xij(1), xij(2), …, xij(p))T。本文重点讨论纵向数据下的广义线性模型(GLM)。假定yij的边际均值和边际方差函数分别为

$ {\mu _{ij}} = E\left( {{y_{ij}}\left| {{\mathit{\boldsymbol{x}}_{ij}}} \right.} \right) = g\left( {{\mathit{\boldsymbol{\beta }}^{\rm{T}}}{\mathit{\boldsymbol{x}}_{ij}}} \right) $
$ Var\left( {{y_{ij}}\left| {{\mathit{\boldsymbol{x}}_{ij}}} \right.} \right) = \phi v\left( {{\mu _{ij}}} \right)\;\;\;\;\;\left( {i = 1, \cdots ,n;j = 1, \cdots ,{m_i}} \right) $

其中β为未知的p×1回归系数,g(·)和v(·)分别是已知的逆连接函数和方差函数,ϕ为尺度参数。不失一般性,假设不同个体间的观测是相互独立的,且同一个个体内部的观测具有一定的相关关系。

Lv等[10]基于指数得分函数提出ERGEE方法,通过求解如式(1)所示的估计方程得到了参数β的基于GEE有效且稳健的估计${\mathit{\boldsymbol{\hat \beta }}_{{\rm{ERGEE}}}}$

$ \mathit{\boldsymbol{U}}_{\rm{n}}^{{\rm{ERGEE}}}\left( {\mathit{\boldsymbol{\beta }},\alpha } \right) = \sum\limits_{i = 1}^n {\mathit{\boldsymbol{D}}_i^{\rm{T}}\mathit{\boldsymbol{V}}_i^{ - 1}h_i^\gamma \left( {{\mu _i}\left( \mathit{\boldsymbol{\beta }} \right)} \right)} = 0 $ (1)

其中Di=∂μi/∂βni×p矩阵;Vi=Ri(α)Ai1/2Ri(α)是依赖于参数αmi×mi工作相关矩阵,且Ai=φdiag(v(μi1), …, v(μimi)), μi=(μi1, …, μimi)T; hiγ(μi(β))=Wi[ψγ(μi(β))-Ci(μi(β))], ψγ(μi)= ψγ(Ai-1/2(Yi-μi)), Ci(μi)=E[ψγ(μi)],ψγ(t)= $ - \frac{{2t}}{\gamma }{\rm{exp}}\left( { - {\mathit{t}^{\rm{2}}}\gamma } \right)$(γ>0)为调节参数。

权重矩阵Wi=diag(wi1, …, wimi)用于限制协变量中离群点的影响,可取Mallow s-type函数wij

$ {w_{ij}} = w\left( {{\mathit{\boldsymbol{x}}_{ij}}} \right) = \min \left\{ {1,{{\left( {\frac{{{b_0}}}{{{{\left( {{\mathit{\boldsymbol{x}}_{ij}} - {\mathit{\boldsymbol{m}}_x}} \right)}^{\rm{T}}}\mathit{\boldsymbol{S}}_x^{ - 1}\left( {{\mathit{\boldsymbol{x}}_{ij}} - {\mathit{\boldsymbol{m}}_x}} \right)}}} \right)}^{\rho /2}}} \right\} $

其中ρ≥1,b0是自由度为p的卡方分布的0.95分位点,mxSx分别表示关于Xij的位置参数和尺度参数的稳健估计,可以利用minimum covariance determinant(MCD)进行估计。

为了简单起见,假定重复观测次数是相同的,即mi=m<∞, 并令R(α)为公共的工作相关矩阵。当重复观测次数不同时,可采用Xue等[11]的方法加以处理。在实际计算中,为了得到${{\mathit{\boldsymbol{\hat \beta }}}_{{\rm{ERGEE}}}}$,ERGEE方法首先要对参数αϕ进行估计,同时Lv等[10]也指出若能得到参数αϕ$\sqrt {n - } $相合估计,则${{\mathit{\boldsymbol{\hat \beta }}}_{{\rm{ERGEE}}}}$β$\sqrt {n - } $相合估计。但即使在一些简单的情形下,讨厌参数α的相合估计也可能不存在;而且对参数α进行估计还会降低参数β估计的效率。为了弥补这一缺陷同时提高ERGEE估计的效率,根据Qu等[2]提出的QIF方法,本文利用一组基矩阵的线性组合来逼近R-1(α),即

$ {\mathit{\boldsymbol{R}}^{ - 1}} \approx {a_1}{\mathit{\boldsymbol{M}}_1} + {a_2}{\mathit{\boldsymbol{M}}_2} + \cdots + {a_k}{\mathit{\boldsymbol{M}}_k} $ (2)

其中M1是单位矩阵,M2,…, Mk是某些已知的对称的基矩阵[2]a1, a2,…, ak是未知常数。

将式(2)代入式(1)构造扩展得分向量为

$ {{\mathit{\boldsymbol{\bar U}}}_n}\left( \mathit{\boldsymbol{\beta }} \right) = \frac{1}{n}\sum\limits_{i = 1}^n {{\mathit{\boldsymbol{U}}_i}\left( \mathit{\boldsymbol{\beta }} \right)} = \frac{1}{n}\sum\limits_{i = 1}^n \\ {\left( {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{D}}_i^{\rm{T}}\mathit{\boldsymbol{A}}_i^{ - 1/2}{\mathit{\boldsymbol{M}}_1}h_i^\gamma \left( {{\mathit{\boldsymbol{\mu }}_i}\left( \mathit{\boldsymbol{\beta }} \right)} \right)}\\ \vdots \\ {\mathit{\boldsymbol{D}}_i^{\rm{T}}\mathit{\boldsymbol{A}}_i^{ - 1/2}{\mathit{\boldsymbol{M}}_k}h_i^\gamma \left( {{\mathit{\boldsymbol{\mu }}_i}\left( \mathit{\boldsymbol{\beta }} \right)} \right)} \end{array}} \right)} $

由于Un(β)的维数为kp,显然高于β的维数p,因此求解方程Un(β)=0会出现过度识别问题。基于GMM估计的思想,定义有效且稳健的ERQIF估计方程如式(3)

$ {Q_n}\left( \mathit{\boldsymbol{\beta }} \right) = n\mathit{\boldsymbol{\bar U}}_n^{\rm{T}}\left( \mathit{\boldsymbol{\beta }} \right)C_n^{ - 1}\left( \mathit{\boldsymbol{\beta }} \right){{\mathit{\boldsymbol{\bar U}}}_n}\left( \mathit{\boldsymbol{\beta }} \right) $ (3)

其中,${C_n}\left( \mathit{\boldsymbol{\beta }} \right) = \frac{1}{n}\sum {{\mathit{\boldsymbol{U}}_i}} \left( \mathit{\boldsymbol{\beta }} \right)\mathit{\boldsymbol{U}}_i^{\rm{T}}\left( \mathit{\boldsymbol{\beta }} \right)$

通过极小化目标函数Qn(·)得到β的估计

$ {{\mathit{\boldsymbol{\hat \beta }}}_{{\rm{ERQIF}}}} = \arg {\min _\mathit{\boldsymbol{\beta }}}{Q_n}\left( \mathit{\boldsymbol{\beta }} \right) $

利用如式(4)的Newto n-Raphson迭代算法得到估计${{\mathit{\boldsymbol{\hat \beta }}}_{{\rm{ERQIF}}}}$[2]

$ {{\mathit{\boldsymbol{\hat \beta }}}^{\left( {k + 1} \right)}} = {{\mathit{\boldsymbol{\hat \beta }}}^{\left( k \right)}} - {\left[ {{\nabla ^2}{Q_n}\left( {{{\mathit{\boldsymbol{\hat \beta }}}^{\left( k \right)}}} \right)} \right]^{ - 1}}\nabla {Q_n}\left( {{{\mathit{\boldsymbol{\hat \beta }}}^{\left( k \right)}}} \right) $ (4)

其中,∇Qn(β)和∇2Qn(β)分别为Qn(β)对β的一阶和二阶导数,且满足

$ \left\{ \begin{array}{l} \nabla {Q_n}\left( \mathit{\boldsymbol{\beta }} \right) = 2n\nabla \mathit{\boldsymbol{\bar U}}_n^{\rm{T}}\left( \mathit{\boldsymbol{\beta }} \right)C_n^{ - 1}\left( \mathit{\boldsymbol{\beta }} \right){{\mathit{\boldsymbol{\bar U}}}_n}\left( \mathit{\boldsymbol{\beta }} \right) - \\ \;\;\;\;\;\;n\mathit{\boldsymbol{\bar U}}_n^{\rm{T}}\left( \mathit{\boldsymbol{\beta }} \right)\nabla C_n^{ - 1}\left( \mathit{\boldsymbol{\beta }} \right){{\mathit{\boldsymbol{\bar U}}}_n}\left( \mathit{\boldsymbol{\beta }} \right)\\ {\nabla ^2}{Q_n}\left( \mathit{\boldsymbol{\beta }} \right) = 2n\nabla \mathit{\boldsymbol{\bar U}}_n^{\rm{T}}\left( \mathit{\boldsymbol{\beta }} \right)C_n^{ - 1}\left( \mathit{\boldsymbol{\beta }} \right)\nabla {{\mathit{\boldsymbol{\bar U}}}_n}\left( \mathit{\boldsymbol{\beta }} \right) + {O_p}\left( 1 \right) \end{array} \right. $
1.2 估计的渐近性质

β0为参数β的真值,β所处的参数空间记为Θ。首先给出一些假设条件[4](本文中Eβ0(·)均表示随机变量关于真值β0的数学期望):

1) 参数空间Θ是可识别的,若ββ0, 则Eβ0[U1(β)]≠0;

2) 对任意βΘEβ0[U1(β)]存在且有限,同时关于β连续;

3) Θ是紧空间,且β0Θ的一个内点;

4) 存在β0的某个邻域O0,当βO0时,有

$ {C_n}\left( \mathit{\boldsymbol{\beta }} \right)\xrightarrow{{a.s.}}{E_{{\mathit{\boldsymbol{\beta }}_0}}}\left[ {{\mathit{\boldsymbol{U}}_1}\left( \mathit{\boldsymbol{\beta }} \right)\mathit{\boldsymbol{U}}_1^{\rm{T}}\left( \mathit{\boldsymbol{\beta }} \right)} \right] \triangleq {\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{{\mathit{\boldsymbol{\beta }}_0}}}\left( \mathit{\boldsymbol{\beta }} \right) $

其中,Σβ0(β)为正定矩阵,且关于βO0连续;

5) Un(β)关于β的一阶导数存在且连续,且当${\mathit{\boldsymbol{\beta }}^{ * }}\xrightarrow{P}{\mathit{\boldsymbol{\beta }}_0}$时,有

$ \frac{{\partial {{\mathit{\boldsymbol{\bar U}}}_n}\left( {{\mathit{\boldsymbol{\beta }}^ * }} \right)}}{{\partial \mathit{\boldsymbol{\beta }}}}\xrightarrow{P}\mathit{\boldsymbol{B}}\left( {{\mathit{\boldsymbol{\beta }}_0}} \right) = {E_{{\mathit{\boldsymbol{\beta }}_0}}}\left[ {\frac{{\partial {{\mathit{\boldsymbol{\bar U}}}_n}\left( {{\mathit{\boldsymbol{\beta }}_0}} \right)}}{{\partial \mathit{\boldsymbol{\beta }}}}} \right] $

基于条件1)~4),根据文献[12]和文献[2]中相关定理的证明,可得定理1。

定理1  在条件1)~4)下,通过式(3)得到的${\mathit{\boldsymbol{\hat \beta }}}_\rm{ERQIF}$是存在的,且当n→∞时,有${\mathit{\boldsymbol{\hat \beta }}_{\rm{ERQIF}}\xrightarrow{a.s.}\mathit{\boldsymbol{\beta }}_{0}}$,即得到了β0的相合估计${{\mathit{\boldsymbol{\hat \beta }}}_{{\rm{ERQIF}}}}$

定理2  在条件1)~5)下,若Un(β)和Cn(β)的二阶导数有有限的均值和方差,则当n→∞时,有

$ \sqrt n \left( {{{\mathit{\boldsymbol{\hat \beta }}}_{{\rm{ERQIF}}}} - {\mathit{\boldsymbol{\beta }}_0}} \right)\xrightarrow{D}{N_p}\left( {0,{\mathit{\boldsymbol{J}}^{ - 1}}\left( {{\mathit{\boldsymbol{\beta }}_0}} \right)} \right) $

即证明了所得估计的渐进正态性。其中,J(β0)=BT(β0)Σβ0-1(β0)B(β0)。

证明  因${{\mathit{\boldsymbol{\hat \beta }}}_{{\rm{ERQIF}}}} = \arg \;\;\mathop {\min }\limits_\beta \;{Q_n}\left( \mathit{\boldsymbol{\beta }} \right)$,则∇Qn(${{\mathit{\boldsymbol{\hat \beta }}}_{{\rm{ERQIF}}}}$)=0。由Taylor公式得到

$ {\bf{0}} = \nabla {Q_n}\left( {{{\mathit{\boldsymbol{\hat \beta }}}_{{\rm{ERQIF}}}}} \right) = \nabla {Q_n}\left( {{\mathit{\boldsymbol{\beta }}_{\rm{0}}}} \right) + {\nabla ^2}{Q_n}\left( {\mathit{\boldsymbol{\tilde \beta }}} \right)\left( {{{\mathit{\boldsymbol{\hat \beta }}}_{{\rm{ERQIF}}}} - {\mathit{\boldsymbol{\beta }}_0}} \right) $

其中${\mathit{\boldsymbol{\tilde \beta }}}$位于${{\mathit{\boldsymbol{\hat \beta }}}_{{\rm{ERQIF}}}}$β0之间。从而有

$ \sqrt n \left( {{{\mathit{\boldsymbol{\hat \beta }}}_{{\rm{ERQIF}}}} - {\mathit{\boldsymbol{\beta }}_0}} \right) = - {\left( {\frac{1}{{2n}}{\nabla ^2}{Q_n}\left( {\mathit{\boldsymbol{\tilde \beta }}} \right)} \right)^{ - 1}}\frac{1}{{2\sqrt n }}\nabla {Q_n}\left( {{\mathit{\boldsymbol{\beta }}_{\rm{0}}}} \right) $

其中

$ \nabla {Q_n}\left( {{\mathit{\boldsymbol{\beta }}_{\rm{0}}}} \right) = 2n\nabla \mathit{\boldsymbol{\bar U}}_n^{\rm{T}}\left( {{\mathit{\boldsymbol{\beta }}_0}} \right)C_n^{ - 1}\left( {{\mathit{\boldsymbol{\beta }}_0}} \right)\nabla {{\mathit{\boldsymbol{\bar U}}}_n}\left( {{\mathit{\boldsymbol{\beta }}_0}} \right) + {O_p}\left( 1 \right) $
$ {\nabla ^2}{Q_n}\left( {\mathit{\boldsymbol{\tilde \beta }}} \right) = 2n\nabla \mathit{\boldsymbol{\bar U}}_n^{\rm{T}}\left( {\mathit{\boldsymbol{\tilde \beta }}} \right)C_n^{ - 1}\left( {\mathit{\boldsymbol{\tilde \beta }}} \right)\nabla {{\mathit{\boldsymbol{\bar U}}}_n}\left( {\mathit{\boldsymbol{\tilde \beta }}} \right) + {O_p}\left( 1 \right) $

由条件4)和5)以及$\sqrt n {{\mathit{\boldsymbol{\bar U}}}_n}\left( {{\mathit{\boldsymbol{\beta }}_0}} \right)\xrightarrow{D}{N_{kp}}\left( {{\bf{0}}, {\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{{\mathit{\boldsymbol{\beta }}_0}}}\left( \mathit{\boldsymbol{\beta }}_0 \right)} \right)$,可得

$ \begin{gathered} \frac{1}{{2\sqrt n }}{Q_n}\left( {{\mathit{\boldsymbol{\beta }}_{\rm{0}}}} \right) = \nabla \mathit{\boldsymbol{\bar U}}_n^{\rm{T}}\left( {{\mathit{\boldsymbol{\beta }}_0}} \right)C_n^{ - 1}\left( {{\mathit{\boldsymbol{\beta }}_0}} \right)\sqrt n {{\mathit{\boldsymbol{\bar U}}}_n}\left( {{\mathit{\boldsymbol{\beta }}_0}} \right) + \hfill \\ {Q_p}\left( {{n^{ - 1/2}}} \right) = {\mathit{\boldsymbol{B}}^{\rm{T}}}\left( {{\mathit{\boldsymbol{\beta }}_0}} \right)\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{{\mathit{\boldsymbol{\beta }}_0}}^{ - 1}\left( {{\mathit{\boldsymbol{\beta }}_0}} \right)\sqrt n {{\mathit{\boldsymbol{\bar U}}}_n}\left( {{\mathit{\boldsymbol{\beta }}_0}} \right) + {Q_p}\left( {{n^{ - 1/2}}} \right)\xrightarrow{D} \hfill \\ {N_p}\left( {{\bf{0}},\mathit{\boldsymbol{J}}\left( {{\mathit{\boldsymbol{\beta }}_0}} \right)} \right)\frac{1}{{2n}}{\nabla ^2}{Q_n}\left( {\mathit{\boldsymbol{\tilde \beta }}} \right) = \nabla \mathit{\boldsymbol{\bar U}}_n^{\rm{T}}\left( {\mathit{\boldsymbol{\tilde \beta }}} \right)C_n^{ - 1}\left( {\mathit{\boldsymbol{\tilde \beta }}} \right)\nabla {{\mathit{\boldsymbol{\bar U}}}_n} \hfill \\ \left( {\mathit{\boldsymbol{\tilde \beta }}} \right) + {O_p}\left( 1 \right)\xrightarrow{P}\mathit{\boldsymbol{J}}\left( {{\mathit{\boldsymbol{\beta }}_0}} \right) \hfill \\ \end{gathered} $

因此可得当n→∞时,$\sqrt n \left( {{{\mathit{\boldsymbol{\hat \beta }}}_{{\rm{ERQIF}}}} - {\mathit{\boldsymbol{\beta }}_0}} \right)\xrightarrow{D}{N_p}$(0, J-1(β0))。

推论1  在定理2的条件下,可得${{\mathit{\boldsymbol{\hat \beta }}}_{{\rm{ERQIF}}}}$渐近协方差的相合估计${\rm{\hat Cov}}\left( {{{\mathit{\boldsymbol{\hat \beta }}}_{{\rm{ERQIF}}}}} \right) = {\left( {\mathit{\boldsymbol{\hat \beta }}_n^{\rm{T}}\mathit{\boldsymbol{ \boldsymbol{\hat \varSigma} }}_n^{ - 1}{{\mathit{\boldsymbol{\hat \beta }}}_n}} \right)^{ - 1}}$,其中

$ {{\mathit{\boldsymbol{\hat B}}}_n} = \sum\limits_{i = 1}^n {{{\left( {{{\mathit{\boldsymbol{\hat B}}}_{i1}}, \cdots ,{{\mathit{\boldsymbol{\hat B}}}_{ik}}} \right)}^{\rm{T}}}} $
$ {{\mathit{\boldsymbol{ \boldsymbol{\hat \varSigma} }}}_n} = \sum\limits_{i = 1}^n {{{\left( {{{\mathit{\boldsymbol{ \boldsymbol{\hat \varSigma} }}}_{i1}}, \cdots ,{{\mathit{\boldsymbol{ \boldsymbol{\hat \varSigma} }}}_{ik}}} \right)}^{\rm{T}}}} $
$ \mathit{\boldsymbol{ \boldsymbol{\hat \varSigma} }}_{is}^{\rm{T}} = \mathit{\boldsymbol{D}}_i^{\rm{T}}\mathit{\boldsymbol{A}}_i^{ - 1/2}{\mathit{\boldsymbol{M}}_s}h_i^\gamma \left( {{\mu _i}\left( {\mathit{\boldsymbol{\hat \beta }}} \right)} \right)h_i^\gamma {\left( {{\mu _i}\left( {\mathit{\boldsymbol{\hat \beta }}} \right)} \right)^{\rm{T}}}{\mathit{\boldsymbol{M}}_s}\mathit{\boldsymbol{A}}_i^{ - 1/2}{\mathit{\boldsymbol{D}}_i} $
$ \mathit{\boldsymbol{\hat B}}_{is}^{\rm{T}} = \mathit{\boldsymbol{D}}_i^{\rm{T}}\mathit{\boldsymbol{A}}_i^{ - 1/2}{\mathit{\boldsymbol{M}}_s}\dot h_i^\gamma \left( {{\mu _i}\left( {\mathit{\boldsymbol{\hat \beta }}} \right)} \right){\mathit{\boldsymbol{D}}_i}\left( {s = 1, \cdots ,k} \right) $
$ \dot h_i^\gamma \left( {{\mathit{\boldsymbol{\mu }}_i}\left( {\mathit{\boldsymbol{\hat \beta }}} \right)} \right) = \partial h_i^\gamma \left( {{\mathit{\boldsymbol{\mu }}_i}\left( {\mathit{\boldsymbol{\hat \beta }}} \right)} \right)/\partial {\mathit{\boldsymbol{\mu }}_i}\left| {_{{\mathit{\boldsymbol{\mu }}_i} = {\mathit{\boldsymbol{\mu }}_i}\left( {\mathit{\boldsymbol{\hat \beta }}} \right)}} \right. $
$ \mathit{\boldsymbol{\hat \beta }} = {{\mathit{\boldsymbol{\hat \beta }}}_{{\rm{ERQIF}}}} $
2 ERQIF算法 2.1 尺度参数和调节参数

要通过式(3)得到ERQIF估计${{\mathit{\boldsymbol{\hat \beta }}}_{{\rm{ERQIF}}}}$,首先需要估计尺度参数ϕ,并且对调节参数γ进行选择,γ的选取对估计的有效性和稳健性有至关重要的影响。假设${\mathit{\boldsymbol{\tilde \beta }}}$β当前的估计,通过绝对偏差中位获取ϕ的稳健估计

$ \hat \phi = {\left\{ {1.483{\rm{medi}}\left\{ {\left| {{{\hat \xi }_{ij}} - {\rm{medi}}\left( {{{\hat \xi }_{ij}}} \right)} \right|} \right\}} \right\}^2} $

其中${{\hat \xi }_{ij}}$=Aij-1/2(Yij-μij(${\mathit{\boldsymbol{\tilde \beta }}}$)),Aij是矩阵Ai对角线上第j个元素。同时通过最小化${\rm{\hat Cov}}\left( {\mathit{\boldsymbol{\hat \beta }}} \right)$的行列式[13]得到当前最优的γ

$ {\gamma _{{\rm{opt}}}} = \arg \mathop {\min }\limits_\gamma \left( {\det \left( {{\rm{\hat Cov}}\left( {\mathit{\boldsymbol{\tilde \beta }}} \right)} \right)} \right) $

其中(${\rm{\hat Cov}}\left( {\mathit{\boldsymbol{\tilde \beta }}} \right)$)由推论1给出。

2.2 估计算法

通过一个迭代算法得到ERQIF估计,具体步骤如下:

1) 给定β的初始估计${{\mathit{\boldsymbol{\hat \beta }}}^{\left( 0 \right)}}$以及迭代收敛阈值ε,并令k=0;

2) 利用当前估计${{\mathit{\boldsymbol{\hat \beta }}}^{\left( k \right)}}$来估计尺度参数${{\hat \phi }^{\left( k \right)}}$和调节参数γopt(k);

3) 利用公式(4)迭代求解${{\mathit{\boldsymbol{\hat \beta }}}^{\left( {k + 1} \right)}}$;

4) 若$\left\| {{{\mathit{\boldsymbol{\hat \beta }}}^{\left( {k + 1} \right)}} - {{\mathit{\boldsymbol{\hat \beta }}}^{\left( k \right)}}} \right\| < \varepsilon $, 则终止计算,并令${{\mathit{\boldsymbol{\hat \beta }}}_{{\rm{ERQIF}}}}$= ${{\mathit{\boldsymbol{\hat \beta }}}^{\left( {k + 1} \right)}}$;否则令kk+1, 返回步骤2)。

在数值计算中,利用独立工作相关矩阵的GEE方法[1]得到初始估计${{\mathit{\boldsymbol{\hat \beta }}}^{\left( 0 \right)}}$,并取ε=10-3,计算结果表明一般在迭代50步内即可收敛。

3 数值模拟

为了得到估计的有限样本性质,本文通过随机模拟考察数据出现离群点时估计的稳健性和工作相关矩阵被正确或错误指定对估计效率的影响。用本文提出的ERQIF方法与GEE、ERGEE方法分别计算3种参数估计的样本偏差(Bias)和标准差(SD)并进行比较。设样本容量n=100,重复观测次数m=5, 模拟次数为500。在模拟计算中,分别指定组内工作相关矩阵Ri(α)为Exch和AR(1)两种结构。

3.1 线性模型

当响应变量是连续型变量时,假设模型满足:

$ {y_{ij}} = \sum\limits_{k = 1}^3 {x_{ij}^{\left( k \right)}{\mathit{\boldsymbol{\beta }}_{0k}}} + {\varepsilon _{ij}}\left( {i = 1, \cdots ,n;j = 1, \cdots ,m} \right) $ (5)

其中,β01=0.7,β02=0.7,β03=-0.4;xij(k)~N(0, 1);Cov(xij(k), xij(l))=0.5|k-l|k, l=1, 2, 3。随机误差εi=(εi1, …, εi5)T的真实工作相关矩阵是Exch结构,相关系数α=0.7。

分别考虑两种情况:

εij~N(0, 1);

εij~0.9N(0, 1)+0.1N(0, 100)。

用3种方法计算得到两种情况下两种结构的估计值,结果分别如表 1表 2所示。

下载CSV 表 1 情况①下参数的估计结果 Table 1 Parameter estimation results for ①
下载CSV 表 2 情况②下参数的估计结果 Table 2 Parameter estimation results for ②

表 1可以看出,不管相关结构是否正确指定,3种方法的Bias差异不大,说明3种方法均能得到参数的相合估计。由GEE和ERGEE方法相比较可知,ERQIF方法的SD均较小,特别是当工作相关结构被错误指定时SD有较大降低,说明ERQIF方法可避免对讨厌参数的估计,同时提高估计的效率。

表 2可以看出,GEE方法的SD值受离群点的影响很大,说明该方法不能得到参数的稳健估计。通过引入加权指数得分函数得到的ERGEE和ERQIF方法可以有效控制离群点的影响,从而保证了估计的稳健性。同时从SD的结果可以看出,不管相关结构是否正确指定,ERQIF方法均优于ERGEE方法,其估计效率有了显著提高。

3.2 逻辑回归模型

当响应变量是0-1型变量时,假设模型满足:

$ \ln \left( {{\mu _{ij}}/\left( {1 - {\mu _{ij}}} \right)} \right) = {x_{ij}}{\beta _0}\left( {i = 1, \cdots ,n;j = 1, \cdots ,m} \right) $ (6)

其中xij独立地服从[0, 1]上的均匀分布,β0=0.4,Yi的真实相关矩阵Ri(α)为可交换结构,且α=0.2。在模拟计算中,分别指定Ri(α)为Exch和AR(1)两种结构;参考文献[14]生成不同结构下的yij;为了研究稳健性,对部分数据进行扰动,即将Xij(i=1, …, 10)的每个值都减去5,计算得到3种方法下参数的偏差和标准差如表 3所示。

下载CSV 表 3 参数β0的估计结果 Table 3 Parameter estimation results for β0

表 3可以看出,3.2节和3.1节的模拟结果相似。在不同情况下,ERQIF方法在Bias与SD上均比其他方法占优势,特别当数据有扰动时,ERGEE与ERQIF方法比非稳健GEE方法的Bias有显著降低;且ERQIF较ERGEE的SD更小,且在工作相关结构被错误指定的情况下,ERQIF也能得到更稳健的估计,估计效率更高。

总之,从3.1节和3.2节的模拟结果可以得出,当数据出现离群点或工作相关矩阵被错误指定时,ERQIF方法均有更好的稳健性和有效性。

4 结束语

本文将加权指数得分函数和QIF方法相结合,得到了纵向数据下广义线性模型的ERQIF估计,证明了估计的渐近性质,并通过数值计算得到了稳健的模拟结果。由于本文仅研究了模型参数估计问题,未来可望进一步研究各种半参数模型的稳健估计和变量选择问题,并对估计的稳健性进行理论分析。

参考文献
[1]
Liang K Y, Zeger S L. Longitudinal data analysis using generalized linear models[J]. Biometrika, 1986, 73(1): 13-22. DOI:10.1093/biomet/73.1.13
[2]
Qu A, Lindsay B G, Li B. Improving generalised estimating equations using quadratic inference functions[J]. Biometrika, 2000, 87(4): 823-836. DOI:10.1093/biomet/87.4.823
[3]
Lian H, Liang H, Wang L. Generalized additive partial linear models for clustered data with diverging number of covariates using GEE[J]. Statistica Sinica, 2014, 24: 173-196.
[4]
Qu A, Li R Z. Quadratic inference functions for varying-coefficient models with longitudinal data[J]. Biometrics, 2006, 62(2): 379-391. DOI:10.1111/j.1541-0420.2005.00490.x
[5]
Tian R Q, Xue L G, Liu C L. Penalized quadratic inference functions for semiparametric varying coefficient partially linear models with longitudinal data[J]. Journal of Multivariate Analysis, 2014, 132: 94-110. DOI:10.1016/j.jmva.2014.07.015
[6]
Fan Y L, Qin G Y, Zhu Z Y. Variable selection in robust regression models for longitudinal data[J]. Journal of Multivariate Analysis, 2012, 109: 156-167. DOI:10.1016/j.jmva.2012.03.007
[7]
Zheng X Y, Fung W K, Zhu Z Y. Robust estimation in joint mean-covariance regression model for longitudinal data[J]. Annals of the Institute of Statistical Mathematics, 2013, 65(4): 617-638. DOI:10.1007/s10463-012-0383-8
[8]
Qin G Y, Zhu Z Y, Fung W K. Robust estimation of generalized partially linear model for longitudinal data with dropouts[J]. Annals of the Institute of Statistical Mathematics, 2016, 68(5): 977-1000. DOI:10.1007/s10463-015-0519-8
[9]
Wang X Q, Jiang Y L, Huang M, et al. Robust variable selection with exponential squared loss[J]. Journal of the American Statistical Association, 2013, 108(502): 632-643. DOI:10.1080/01621459.2013.766613
[10]
Lv J, Yang H, Guo C H. An efficient and robust variable selection method for longitudinal generalized linear models[J]. Computational Statistics & Data Analysis, 2015, 82: 74-88.
[11]
Xue L, Qu A N, Zhou J H. Consistent model selection for marginal generalized additive model for correlated data[J]. Journal of the American Statistical Association, 2010, 105(492): 1518-1530. DOI:10.1198/jasa.2010.tm10128
[12]
Hansen L P. Large sample properties of generalized method of moments estimators[J]. Econometrica, 1982, 50(4): 1029-1054. DOI:10.2307/1912775
[13]
Yao W X, Lindsay B G, Li R Z. Local modal regression[J]. Journal of Nonparametric Statistics, 2012, 24(3): 647-663. DOI:10.1080/10485252.2012.678848
[14]
Oman S D. Easily simulated multivariate binary distributions with given positive and negative correlations[J]. Computational Statistics & Data Analysis, 2009, 53(4): 999-1005.