引言

纵向数据经常出现在经济学和生物医学研究中，它通过对相同的个体进行多次重复测量得到，因而具有组内相关、组间独立的特点；同时由于重复测量，一个个体的突变往往会使数据产生一系列的异常值，因此如何处理组内相关性，提高估计的效率和稳健性，一直是纵向数据研究的重要内容。Koenker^[1-2]提出的分位数回归是研究纵向数据的一种重要方法，该方法不仅能实现稳健估计，还可描述响应变量的整个条件分布。纵向数据下的分位数回归最大的困难是估计相关矩阵，许多学者利用特定的相关结构^[3-5]，或结合经验似然(EL)^[4]、二次推断函数(QIF)^[5]、高斯copula函数以及改进Cholesky分解(MCD)法^[6]等方法将数据相关性纳入估计方程中，从而得到更有效的估计。

当数据相互独立时，Zou等^[7]提出了复合分位数回归(CQR)估计，该方法充分利用多个分位数信息，从而进一步提高了参数估计的效率，因此CQR方法被广泛应用于各种复杂数据中。在纵向数据下，Tang等^[8]利用EL和QIF方法构造权重，得到参数的加权CQR估计；Lv等^[6]利用MCD分解对相关数据的协方差矩阵进行估计，并与CQR方法相结合得到了参数的有效估计; Zhao等^[9]利用广义矩估计(GMM)方法^[10]将CQR和QIF方法相结合，得到了参数的复合分位数二次推断函数估计(CQRQIF)。

统计推断的另一个重要问题是对兴趣参数置信域的构造，通过EL得到的置信域具有一些很好的性质，如域保持性、变换不变性以及置信域的形状完全由数据决定等。Zhao等^[11-12]基于CQR和EL方法研究了线性模型和部分线性模型下参数的稳健估计，并利用经验对数似然比的渐近卡方分布得到参数的置信域，该方法无需估计渐近方差，且比一些已有方法更加有效和稳健，但是它们仅考虑独立数据的情形。

在纵向数据下，本文进一步研究基于CQR的线性模型的EL推断，利用QIF的思想将工作相关矩阵的逆表示为某些基矩阵的线性组合，从而将数据组内相关性纳入估计方程中；将CQR和EL方法相结合，得到了参数的复合分位数回归经验似然估计(CQREL)估计，证明了所得估计的渐近正态性；另外，无需预先估计渐近方差，利用检验统计量服从渐近卡方分布，得到兴趣参数的置信域；数值模拟结果进一步说明所得CQREL估计方法能更有效和更稳健地估计参数。

1 参数估计方法 1.1 CQREL方法

假设纵向数据有n个个体, y_i=(y_i1, …, y_imi)^T和X_i=(X_i1, …, X_imi)^T分别为个体i的响应变量和协变量，其中X_ij=(x_ij1, …, x_ijp)^T，i=1, …, n，j=1, …, m_i。假设m_i与i=1, …, n一致有界，则$M = \sum\limits_{i = 1}^n {{m_i}} $与n同阶。考虑如下线性模型

$ {y_{ij}} = \mathit{\boldsymbol{X}}_{ij}^{\rm{T}}\mathit{\boldsymbol{\beta }} + {\mathit{\boldsymbol{\varepsilon }}_{ij}} $

(1)

式中，β为未知的回归参数，ε_ij为随机误差，ε_i=(ε_i1, …, ε_imi)^T与X_i相互独立，且其均值为0，即E(ε_i)=0。假设纵向数据满足不同个体间相互独立，且相同个体内存在某种相关性；当忽略数据组内相关性时，令0 < τ₁ < τ₂ < … < τ_K, 可用复合分位数回归(CQR)^[13]方法估计参数

$ (\mathit{\boldsymbol{\tilde \beta }}, \mathit{\boldsymbol{\widetilde b}}) = \mathop {\arg \min }\limits_{\mathit{\boldsymbol{\beta , b}}} \sum\limits_{k = 1}^K {\left( {\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^{{m_i}} {{\rho _{{\tau _k}}}} } \left( {{y_{ij}} - {b_k} - } \right.} \right.} \left. {\mathit{\boldsymbol{X}}_{ij}^{\rm{T}}\mathit{\boldsymbol{\beta }}} \right) $

(2)

式中，ρ_τk(u)=u(τ_k－1(u < 0))，1(·)为示性函数；b=(b₁, …, b_k, …, b_K)^T，b_k=F^－1(τ_k)，F(·)为ε_ij的分布函数。通常，取${\tau _k} = \frac{k}{{K + 1}}, k = 1, \cdots, K$。在R语言中可以利用cqrReg软件包得到估计$\left({\mathit{\boldsymbol{\widetilde \beta }}, \mathit{\boldsymbol{\widetilde b}}} \right)$。若记θ=(β^T, b^T)^T，则由式(2)可知，${\mathit{\boldsymbol{\tilde \theta }}}$满足估计方程

$ \sum\limits_{k = 1}^K {\sum\limits_{i = 1}^n {\mathit{\boldsymbol{Z}}_{ik}^{\rm{T}}} } {\mathit{\boldsymbol{S}}_{ik}}(\mathit{\boldsymbol{\theta }}) = 0 $

式中，Z_ik=(X_i, D_ik), D_ik为第k列全为1且其余列全为0的m_i×K矩阵；S_ik(θ)=ψ_τk(y_i－Z_ikθ), ψ_τk(u)=(ψ_τk(u₁), …, ψ_τk(u_s))^T，u=(u₁, …, u_s)^T为任意的s维向量，ψ_τk(u_i)=τ_k－1(u_i < 0)。

若记θ的真值为θ₀=(β₀^T, b₀^T)^T，则由文献[12]中定理2.1及其证明可知，${\mathit{\boldsymbol{\tilde \theta }}}$是θ₀的$\sqrt n $-相合估计，但因忽略了数据组内相关性，从而导致了估计效率的损失。为了提高估计的效率，需要在估计方程中纳入相关性，基于广义估计方程(GEE)方法^[13]和Jung^[14]的思路，可用方程(3)估计θ

$ \sum\limits_{k = 1}^K {\sum\limits_{i = 1}^n {\mathit{\boldsymbol{Z}}_{ik}^{\rm{T}}} } \mathit{\boldsymbol{V}}_{ik}^{ - 1}{\mathit{\boldsymbol{S}}_{ik}}(\mathit{\boldsymbol{\theta }}) = 0 $

(3)

式中，V_ik=Cov(S_ik(θ))=A_ik^1/2R_i(α)A_ik^1/2，R_i(α)为依赖于讨厌参数α的工作相关矩阵；A_ik=w_k^－1I_mi, w_k^－1=τ_k(1－τ_k), I_mi为m_i×m_i单位矩阵。由于参数α较难估计，受Qu等^[5]的启发，将R_i^－1(α)表示为某些基矩阵的线性组合，即$\mathit{\boldsymbol{R}}_i^{ - 1}(\mathit{\alpha }) = \sum\limits_{l = 1}^L {{a_l}} {\mathit{\boldsymbol{M}}_{il}}$，其中M_il为一些已知的对称基矩阵，a_l为未知常数，l=1, …, L。

若R_i(α)是可交换结构，则R_i^－1(α)可表示为a₁M_i1+a₂M_i2，其中，M_i1为单位矩阵，M_i2为主对角线均为0、其余均为1的矩阵。若R_i(α)是一阶自回归结构(AR1)，则R_i^－1(α)可表示为a₁M_i1+a₂M_i2+a₃M_i3，其中M_i1为单位矩阵，M_i2为主对角线两侧的对角线为1、其余均为0的矩阵，M_i3为(1, 1)和(m_i, m_i)元素为1、其余均为0的矩阵(更多细节可参见文献[5])。故此估计方程式(3)可以看成是如下U_i(θ)(i=1, …, n)的线性组合

$ {\mathit{\boldsymbol{U}}_i}\left( \mathit{\boldsymbol{\theta }} \right) = \sum\limits_{k = 1}^K {{w_k}} \left( {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{Z}}_{ik}^{\rm{T}}{\mathit{\boldsymbol{M}}_{i1}}}\\ \vdots \\ {\mathit{\boldsymbol{Z}}_{ik}^{\rm{T}}{\mathit{\boldsymbol{M}}_{iL}}} \end{array}} \right){\mathit{\boldsymbol{S}}_{ik}}(\mathit{\boldsymbol{\theta }}) = \sum\limits_{k = 1}^K {{w_k}} \mathit{\boldsymbol{G}}_{ik}^{\rm{T}}{\mathit{\boldsymbol{S}}_{ik}}(\mathit{\boldsymbol{\theta }}) $

式中G_ik=(M_i1^TZ_ik, …, M_iL^TZ_ik)。

给定R_i(α)，E(U_i(θ₀))=0，将CQR和EL方法相结合，重新组合U_i(θ₀)(i=1, …, n)的信息，得到θ的经验对数似然比函数为

$ {l_n}(\mathit{\boldsymbol{\theta }}) = - 2\max \left\{ {\sum\limits_{i = 1}^n {\lg } \left( {n{p_i}} \right)|{p_i} \ge 0, \sum\limits_{i = 1}^n {{p_i}} = } \right.{\rm{ 1, }}\left. {\sum\limits_{i = 1}^n {{p_i}} {\mathit{\boldsymbol{U}}_i}(\mathit{\boldsymbol{\theta }}) = 0} \right\} $

对给定的θ，当0位于{U_i(θ), i=1, …, n}的凸包之内时，l_n(θ)存在且唯一。利用拉格朗日乘子法可将l_n(θ)表示为

$ {l_n}(\mathit{\boldsymbol{\theta }}) = 2\sum\limits_{i = 1}^n {\lg } \left[ {1 + {\mathit{\boldsymbol{\lambda }}^{\rm{T}}}{\mathit{\boldsymbol{U}}_i}(\mathit{\boldsymbol{\theta }})} \right] $

式中，λ=λ(θ)为拉格朗日乘子，且满足

$ \sum\limits_{i = 1}^n {\left[ {\frac{{{\mathit{\boldsymbol{U}}_i}(\mathit{\boldsymbol{\theta }})}}{{1 + {\mathit{\boldsymbol{\lambda }}^{\rm{T}}}{\mathit{\boldsymbol{U}}_i}(\mathit{\boldsymbol{\theta }})}}} \right]} = 0 $

最小化l_n(θ)，得到θ的CQREL为

$ {\left( {{{\mathit{\boldsymbol{\widehat \beta }}}^{\rm{T}}}, {{\mathit{\boldsymbol{\widehat b}}}^{\rm{T}}}} \right)^{\rm{T}}} = \mathit{\boldsymbol{\widehat \theta }} = \mathop {\arg \min }\limits_\theta {l_n}(\mathit{\boldsymbol{\theta }}) $

(4)

值得注意的是，无需估计讨厌参数和未知常数a_l, l=1, …, L，可以在R语言中通过软件包emplik和optim来求解最小化问题式(4)。另外，与CQRQIF^[9]相比，本文的CQREL方法不需选择一个正定的权重矩阵，而且构造置信域的形状完全由数据本身决定，无需预先估计渐近方差，特别当渐近方差较难估计时，这一优势尤为重要。

1.2 估计的渐近性质

对任意的i=1, …, n，j=1, …, m_i，假设下列条件成立：

(1) Ω=E(U_i(θ₀)U_i(θ₀)^T)是正定矩阵；

(2) ε_ij的分布函数F(·)是二阶连续可微的，其密度函数f(·)有界，且f(·)>0；

(3) X_ij一致有界;

(4) β所在的参数空间为$\mathscr{B}$，且$\mathscr{B}$为R^p的紧子集，其真值β₀为$\mathscr{B}$的内点。

为了便于证明，引入以下记号：

① 令Δ_ik=f(b_0k)I_mi，Γ=(Γ₁, Γ₂)；

② ${\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}^{\rm{T}}} = \sum\limits_{k = 1}^K {{w_k}} E\left({\mathit{\boldsymbol{Z}}_{ik}^{\rm{T}}{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{ik}}{\mathit{\boldsymbol{M}}_{i1}}{\mathit{\boldsymbol{Z}}_{ik}}, \cdots, \mathit{\boldsymbol{Z}}_{ik}^{\rm{T}}{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{ik}}{\mathit{\boldsymbol{M}}_{iL}}\left. {{\mathit{\boldsymbol{Z}}_{ik}}} \right)} \right.$;

③ $\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_1^{\rm{T}} = \sum\limits_{k = 1}^K {{w_k}} E\left({\mathit{\boldsymbol{X}}_i^{\rm{T}}{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{ik}}{\mathit{\boldsymbol{M}}_{\mathit{i}{\rm{l}}}}{\mathit{\boldsymbol{Z}}_{ik}}, \cdots, \mathit{\boldsymbol{X}}_i^{\rm{T}}{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{ik}}{\mathit{\boldsymbol{M}}_{iL}}{\mathit{\boldsymbol{Z}}_{ik}}} \right)$;

④ $\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_2^{\rm{T}} = \sum\limits_{k = 1}^K {{w_k}} E\left({\mathit{\boldsymbol{D}}_{ik}^{\rm{T}}{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{ik}}{\mathit{\boldsymbol{M}}_{i1}}{\mathit{\boldsymbol{Z}}_{ik}}, \cdots, \mathit{\boldsymbol{D}}_{ik}^{\rm{T}}{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{ik}}{\mathit{\boldsymbol{M}}_{iL}}{\mathit{\boldsymbol{Z}}_{ik}}} \right)$。

定理1  假设条件式(1)~(4)成立，则当n→∞时，有$\sqrt{n}\left(\mathrm{\hat{ }\!\!\theta\!\!\text{ }}-{{\mathrm{ }\!\!\theta\!\!\text{ }}_{0}} \right)\xrightarrow{D}N(\bf{0}, \mathrm{ }\!\!\Sigma\!\!\text{ })$, 且${{l}_{n}}\left({{\mathrm{ }\!\!\theta\!\!\text{ }}_{0}} \right)\xrightarrow{D}{{\chi }^{2}}(p+K)$，其中$\mathrm{ }\!\!\Sigma\!\!\text{ }={{\left({{\mathrm{ }\!\!\Gamma\!\!\text{ }}^{\text{T}}}{{\mathrm{ }\!\!\Omega\!\!\text{ }}^{-1}}\mathrm{ }\!\!\Gamma\!\!\text{ } \right)}^{-1}}$。

证明：由文献[15]可知，${\mathit{\boldsymbol{\hat \theta }}}$是θ₀的相合估计，用‖·‖表示欧几里德范数。记${\mathit{\boldsymbol{D}}_n}(\mathit{\boldsymbol{\theta }}) = {n^{ - 1}}\sum\limits_{i = 1}^n {{\mathit{\boldsymbol{U}}_i}(\mathit{\boldsymbol{\theta }})} $, 由文献[16]得

$ \begin{array}{l} {\mathit{\boldsymbol{D}}_n}(\mathit{\boldsymbol{\theta }}) = {\mathit{\boldsymbol{D}}_n}\left( {{\mathit{\boldsymbol{\theta }}_0}} \right) - {n^{ - 1}}\mathop \sum \limits_{i = 1}^n \\ \left( {\begin{array}{*{20}{c}} {\sum\limits_{k = 1}^K {{w_k}} \mathit{\boldsymbol{Z}}_{ik}^{\rm{T}}{\mathit{\boldsymbol{M}}_{i{\rm{l}}}}{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{ik}}{\mathit{\boldsymbol{Z}}_{ik}}}\\ \cdots \\ {\sum\limits_{k = 1}^K {{w_k}} \mathit{\boldsymbol{Z}}_{ik}^{\rm{T}}{\mathit{\boldsymbol{M}}_{iL}}{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{ik}}{\mathit{\boldsymbol{Z}}_{ik}}} \end{array}} \right)\left( {\mathit{\boldsymbol{\theta }} - {\mathit{\boldsymbol{\theta }}_0}} \right) + {o_p}\left( {{n^{ - 1/2}}} \right) = {\mathit{\boldsymbol{D}}_n}\left( {{\mathit{\boldsymbol{\theta }}_0}} \right) - \\ \mathit{\boldsymbol{ \boldsymbol{\varGamma} }}\left( {\mathit{\boldsymbol{\theta }} - {\mathit{\boldsymbol{\theta }}_0}} \right) + {o_p}\left( {{n^{ - 1/2}}} \right) \end{array} $

(5)

且式(5)在Θ_n={θ:‖θ－θ₀‖=O(n^－1/2)}上一致成立。由文献[17]有

$ \mathit{\boldsymbol{\lambda }}(\mathit{\boldsymbol{\theta }}) = {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}^{ - 1}}{\mathit{\boldsymbol{D}}_n}(\mathit{\boldsymbol{\theta }}) + {o_p}\left( {{n^{ - 1/2}}} \right) $

(6)

且式(6)在Θ_n上一致成立。由条件式(1)~(4)和式(5)、(6)，将l_n(θ)泰勒展开得到

$ \begin{array}{l} {l_n}(\mathit{\boldsymbol{\theta }}) = 2\left\{ {\sum\limits_{i = 1}^n \mathit{\boldsymbol{\lambda }} {{(\mathit{\boldsymbol{\theta }})}^{\rm{T}}}{\mathit{\boldsymbol{U}}_i}(\mathit{\boldsymbol{\theta }}) - \frac{1}{2}\sum\limits_{i = 1}^n ( \mathit{\boldsymbol{\lambda }}} \right.\\ \left. {{{\left. {{{(\mathit{\boldsymbol{\theta }})}^{\rm{T}}}{\mathit{\boldsymbol{U}}_i}(\mathit{\boldsymbol{\theta }})} \right)}^2} + {o_p}(1)} \right\} = n{\mathit{\boldsymbol{D}}_n}{(\mathit{\boldsymbol{\theta }})^{\rm{T}}}{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}^{ - 1}}{\mathit{\boldsymbol{D}}_n}(\mathit{\boldsymbol{\theta }}) + \\ {o_p}(1) = n{\left( {{\mathit{\boldsymbol{D}}_n}\left( {{\mathit{\boldsymbol{\theta }}_0}} \right) - \mathit{\boldsymbol{ \boldsymbol{\varGamma} }}\left( {\mathit{\boldsymbol{\theta }} - {\mathit{\boldsymbol{\theta }}_0}} \right)} \right)^{\rm{T}}}{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}^{ - 1}}\left( {{\mathit{\boldsymbol{D}}_n}\left( {{\mathit{\boldsymbol{\theta }}_0}} \right) - } \right.\\ \left. {\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}\left( {\mathit{\boldsymbol{\theta }} - {\mathit{\boldsymbol{\theta }}_0}} \right)} \right) + {o_p}(1) = n{\left( {\mathit{\boldsymbol{\theta }} - {\mathit{\boldsymbol{\theta }}_0}} \right)^{\rm{T}}}{\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}^{\rm{T}}}{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}^{ - 1}}\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}\left( {\mathit{\boldsymbol{\theta }} - {\mathit{\boldsymbol{\theta }}_0}} \right) - \\ 2n{\left( {\mathit{\boldsymbol{\theta }} - {\mathit{\boldsymbol{\theta }}_0}} \right)^{\rm{T}}}{\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}^{\rm{T}}}{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}^{ - 1}}{\mathit{\boldsymbol{D}}_n}\left( {{\mathit{\boldsymbol{\theta }}_0}} \right) + n{\mathit{\boldsymbol{D}}_n}{\left( {{\mathit{\boldsymbol{\theta }}_0}} \right)^{\rm{T}}}{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}^{ - 1}}{\mathit{\boldsymbol{D}}_n}\left( {{\mathit{\boldsymbol{\theta }}_0}} \right) + \\ {o_p}(1) \end{array} $

(7)

由条件式(1)~(4)可知，式(5)、(6)与式(7)在Θ_n上一致成立，且有

$ \sqrt n \left( {\mathit{\boldsymbol{\widehat \theta }} - {\mathit{\boldsymbol{\theta }}_0}} \right) = {\left( {{\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}^{\rm{T}}}{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}^{ - 1}}\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}} \right)^{ - 1}}{\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}^{\rm{T}}}{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}^{ - 1}}\sqrt n {\mathit{\boldsymbol{D}}_n}\left( {{\mathit{\boldsymbol{\theta }}_0}} \right) + {o_p}(1) $

(8)

另外有

$ \sqrt{n}{{\boldsymbol{D}}_{n}}\left( {{\mathrm{ }\!\!{\boldsymbol{\theta }}\!\!\text{ }}_{0}} \right)={{n}^{-\frac{1}{2}}}\sum\limits_{i=1}^{n}{{{\boldsymbol{U}}_{i}}}\left( {{\mathrm{ }\!\!{\boldsymbol{\theta }}\!\!\text{ }}_{0}} \right)\xrightarrow{D}(\bf{0}, \mathrm{ }\!\!{\boldsymbol{\Omega }}\!\!\text{ }) $

(9)

则当n→∞时, 有

$ \sqrt n \left( {\mathit{\boldsymbol{\widehat \theta }} - {\mathit{\boldsymbol{\theta }}_0}} \right)\xrightarrow{D}N\left( {{\bf{0}},{{\left( {{\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}^{\rm{T}}}{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}^{ - 1}}\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}} \right)}^{ - 1}}} \right) $

由式(7)、(9)可得

$ {l_n}\left( {{\mathit{\boldsymbol{\theta }}_0}} \right) = \sqrt n {\mathit{\boldsymbol{D}}_n}{\left( {{\mathit{\boldsymbol{\theta }}_0}} \right)^{\rm{T}}}{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}^{ - 1}}\sqrt n {\mathit{\boldsymbol{D}}_n}\left( {{\mathit{\boldsymbol{\theta }}_0}} \right) + {o_p}(1)\xrightarrow{D}{\chi ^2}(p + K) $

值得注意的是，若文献[9]中式(4)取A_ik=w_k^－1I_mi，则由定理1可知，本文的CQREL与CQRQIF^[9]所得估计有相同的渐近协方差阵。

本文感兴趣的参数是β，所以可以利用profile经验对数似然比检验统计量对β进行检验或构造β的置信域。欲检验H₀:β=β₀, 需定义检验统计量

$ {W_n}\left( {{\mathit{\boldsymbol{\beta }}_0}} \right) = {l_n}\left( {{\mathit{\boldsymbol{\beta }}_0},\mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over b} }}} \right) - {l_n}\left( {{\mathit{\boldsymbol{\beta }}_0},\mathit{\boldsymbol{\widehat b}}} \right) $

式中，$\mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over b} }} = \mathop {\arg \min {l_n}}\limits_b \left({{\mathit{\boldsymbol{\beta }}_0}, \mathit{\boldsymbol{b}}} \right)$。

定理2  假设条件式(1)~(4)成立，H₀:β=β₀，当H₀为真时，有${W_n}\left({{\mathit{\boldsymbol{\beta }}_0}} \right)\xrightarrow{D}{\chi ^2}(p)$。

证明：记A=Ω^－1/2(I_(p+K)L－ΓΣ^－1Γ^TΩ^－1)Ω^1/2，A₂=Ω^－1/2(I_(p+K)L－Γ₂Σ₂^－1Γ₂^TΩ^－1)Ω^1/2，Σ₂=Γ₂^TΩ^－1Γ₂；由式(7)、(8)及文献[17]得

$ {l_n}(\mathit{\boldsymbol{\widehat \beta }}, \mathit{\boldsymbol{\widehat b}}) = {\left[ {\sqrt n {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}^{ - 1/2}}{\mathit{\boldsymbol{D}}_n}\left( {{\mathit{\boldsymbol{\theta }}_0}} \right)} \right]^{\rm{T}}}\mathit{\boldsymbol{A}}\left[ {\sqrt n {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}^{ - 1/2}}{\mathit{\boldsymbol{D}}_n}} \right.\left. {\left( {{\mathit{\boldsymbol{\theta }}_0}} \right)} \right] + {o_p}(1) $

$ {l_n}\left( {{\mathit{\boldsymbol{\beta }}_0}, \mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over b} }}} \right) = {\left[ {\sqrt n {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}^{ - 1/2}}{\mathit{\boldsymbol{D}}_n}\left( {{\mathit{\boldsymbol{\theta }}_0}} \right)} \right]^{\rm{T}}}{\mathit{\boldsymbol{A}}_2}\left[ {\sqrt n {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}^{ - 1/2}}{\mathit{\boldsymbol{D}}_n}} \right.\left. {\left( {{\mathit{\boldsymbol{\theta }}_0}} \right)} \right] + {o_p}(1) $

从而有

$ \begin{array}{l} {W_n}\left( {{\mathit{\boldsymbol{\beta }}_0}} \right) = {\left[ {\sqrt n {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}^{ - 1/2}}{\mathit{\boldsymbol{D}}_n}\left( {{\mathit{\boldsymbol{\theta }}_0}} \right)} \right]^{\rm{T}}}\left( {{\mathit{\boldsymbol{A}}_2} - \mathit{\boldsymbol{A}}} \right)\\ [\sqrt n \left. {{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}^{ - 1/2}}{\mathit{\boldsymbol{D}}_n}\left( {{\mathit{\boldsymbol{\theta }}_0}} \right)} \right] + {o_p}(1) = {\left[ {\sqrt n {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}^{ - 1/2}}{\mathit{\boldsymbol{D}}_n}\left( {{\mathit{\boldsymbol{\theta }}_0}} \right)} \right]^{\rm{T}}}\\ {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}^{ - 1/2}}\left( {\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}^{ - 1}}{\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}^{\rm{T}}} - {\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_2}\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_2^{ - 1}\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_2^{\rm{T}}} \right){\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}^{ - 1/2}}\left[ {\sqrt n {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}^{ - 1/2}}{\mathit{\boldsymbol{D}}_n}\left( {{\mathit{\boldsymbol{\theta }}_0}} \right)} \right] + {o_p}(1) \end{array} $

(10)

由$\sqrt{n} \boldsymbol{\Omega}^{-1 / 2} \boldsymbol{D}_{n}\left(\boldsymbol{\theta}_{0}\right) \stackrel{D}{\longrightarrow}\left(\bf{0}, \boldsymbol{I}_{(p+K) L}\right)$得到

$ \mathit{\boldsymbol{I}}{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}^{ - 1}}{\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}^{\rm{T}}} \ge \left( {{\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_1},{\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_2}} \right)\left( {\begin{array}{*{20}{c}} {\bf{0}}&{\bf{0}}\\ {\bf{0}}&{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_2^{ - 1}} \end{array}} \right)\left( {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_1^{\rm{T}}}\\ {\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_2^{\rm{T}}} \end{array}} \right) = {\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_2}\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_2^{ - 1}\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_2^{\rm{T}} $

式(10)中，Ω^-1/2ΓΣ^－1Γ^TΩ^-1/2与Ω^-1/2Γ₂Σ₂^－1Γ₂^TΩ^-1/2为对称的幂等矩阵，且迹分别为p+K和K，根据文献[18]中的结论，有${W_n}\left({{\mathit{\boldsymbol{\beta }}_0}} \right)\xrightarrow{D}{\chi ^2}(p)$。

由定理2可知，β的置信度为1－α的置信域可近似表示为{β:β∈$\mathscr{B}$, W_n(β)≤χ_1-α²(p)}, 其中χ_1-α²(p)为分布χ²(p)的1－α分位数。

2 数值模拟

考虑以下4种不同的误差分布情况进行模拟。

① 多元正态分布。假设ε_i满足多元正态分布，即ε_i~N(0, Σ_ε(ρ))(Σ_ε(ρ)为一阶自回归结构，相关系数ρ=0.7)。

② 多元t分布。假设ε_i满足自由度为2的t分布，即ε_i~t₂(0, Σ_ε(ρ)), Σ_ε(ρ)同情形①。

③ 有异常值的多元正态分布。从情形①产生的数据中分别随机选取2%的y_ij各加减20。

④ 有异常值的多元t分布。从情形②产生的数据中分别随机选取2%的y_ij各加减20。

模拟中样本容量n分别取100、200，为简单起见，取m_i=m=5，模拟次数为300。数据由如下模型生成

y_ij=X_ij^Tβ+ε_ij

式中，i=1, …, n，j=1, …, m，β=(β₁, β₂)^T，β₁=－1，β₂=2，X_ij=(x_ij1, x_ij2)^T，x_ij1~N(0, 1)，x_ij2~N(0, 1)。

将本文提出的CQREL方法分别与最小二乘经验似然估计(LSEL)和CQRQIF方法进行比较，采用均方误差根(RMSE)评估不同估计${\mathrm{\hat{ }\!\!\beta\!\!\text{ }}}$的表现

$ {x_{{\rm{RusE }}}} = \sqrt {{{(\mathit{\boldsymbol{\widehat \beta }} - \mathit{\boldsymbol{\beta }})}^{\rm{T}}}(\mathit{\boldsymbol{\widehat \beta }} - \mathit{\boldsymbol{\beta }})} $

模拟中，将工作相关矩阵R_i分别指定为独立(WI)、可交换(EX)和一阶自回归(AR1)结构。对基于CQR的估计方法，均选取9个不同的分位数τ₁=0.1, …, τ₉=0.9。

表 1~4给出了4种情形下β估计的RMSE均值(Mean)和标准差(SD)。限于篇幅，表中只列出了n=100的情形，当n=200时也得到了类似的结果, 且RMSE有更小的Mean和SD，同样验证了3种估计的相合性。

由表 1~4可以看出，对于3种方法，忽略了数据相关结构(WI)的估计结果均较差，因此研究纵向数据时考虑组内的相关结构是很有必要的。当误差服从正态分布(表 1)时，LSEL的结果优于CQREL和CQRQIF；当误差服从非正态分布(表 2)或数据存在异常值时(表 3、4)，CQREL和CQRQIF因其RMSE有更小的Mean和SD值，明显优于LSEL。说明CQREL和CQRQIF方法所得估计具有更好的有效性和稳健性。CQREL和CQREL两种方法具有相近的结果，原因主要是两种方法的估计有相同的渐近协方差阵。

接下来考虑参数β的区间估计。基于300次模拟结果，表 5给出了4种情形下(Case1、Case2、Case3、Case4)95%置信域的平均覆盖率(CP)。从表中可以看出，当考虑相关性时，无论哪种情形，基于CQR的两种方法的CP值均接近95%；而LSEL只在情形①下表现较好，当误差服从非正态分布或数据存在异常值时，其CP值均偏离95%。

为了比较CQREL和CQRQIF的置信域性质，图 1给出了情形①下两种方法的置信域。从图中可以看出，基于AR1相关结构的估计比基于EX相关结构有更小的置信区域，并且在CP值都接近95%的情况下，CQREL方法所得的置信域在精度上明显高于CQRQIF。

综上所述，当误差服从非正态分布或数据存在异常值时，本文的CQREL方法比LSEL方法有更好的估计效率和稳健性；另外，尽管CQREL与CQRQIF方法所得估计RMSE的Mean、SD和置信域的CP均有相近的结果，但是CQREL能得到精度更高的置信域，从而比CQRQIF有更好的有限样本性质。

3 结束语

在纵向数据下，本文将CQR和EL方法相结合，得到了模型参数更加有效且稳健的估计。同时在理论上证明了所得估计和经验对数似然比统计量的渐近性质。借助QIF的思想考虑数据相关性，在估计过程中无需估计工作相关矩阵中的讨厌参数，并且所构造的置信域的形状完全由数据本身决定，无需预先估计渐近方差。通过数值模拟得到了更加有效和稳健的参数估计结果，并且在平均覆盖概率接近95%的情况下，构造出精度更高的置信域。

本文仅究了纵向数据下线性模型的统计推断问题，未来可将该方法推广到部分线性变系数模型、单指标模型等半参数模型中，同时考虑响应变量缺失的情况。