SCAD-L2正则化下广义估计方程的性质分析

Analysis of the Properties of Generalized Estimating Equations under SCAD-L2 Regularization

摘要: 广义估计方程(Generalized estimating equations, GEE)因能够有效处理个体内相关性而在纵向数据分析中得到广泛应用。然而，当纵向数据的协变量高度相关时，传统变量选择方法往往面临变量选择不稳定的问题。本文将SCAD-L2正则化项融入GEE框架中，以实现变量选择与参数估计的双重优化。随后，本文提出一种适用于多重共线性纵向数据的初始值选择策略，即使用L₂惩罚下的GEE估计值作为计算初始值。最后，本文研究了SCAD-L2惩罚下GEE估计的大样本渐近性质。模拟实验表明，该方法在纵向数据的参数估计与变量选择中显著优于现有方法，为复杂相关结构下的纵向数据建模提供了有效方法。

Abstract: Generalized estimating equations (GEE) have been widely used in longitudinal data analysis due to their ability to effectively account for within-subject correlation. However, when covariates in longitudinal data are highly correlated, traditional variable selection methods often suffer from instability. In this paper, we incorporate the SCAD-L2 regularization into the GEE framework to simultaneously optimize variable selection and parameter estimation. We then propose an initial value selection strategy for longitudinal data with multicollinearity, which uses the GEE estimator under an L2 penalty as the starting value for computation. Finally, we investigate the large-sample asymptotic properties of the SCAD-L2 penalized GEE estimator. Simulation studies show that the proposed method substantially outperforms existing approaches in both parameter estimation and variable selection, providing an effective tool for modeling longitudinal data with complex correlation structures.

1. 引言

随着数据来源的不断丰富，数据维度呈现爆炸式增长，尤其是在生物医学研究、基因组学和公共卫生等领域，研究数据越来越多地表现出高维纵向的特征。为了处理对同一受试者进行重复测量所产生的复杂相关结构，Liang等[1]提出了广义估计方程(generalized estimating equations, GEE)，该方法已成为纵向数据分析的基础工具，具有较强灵活性，被广泛应用于各领域的纵向数据研究中[2]。Wang [3]针对具有大规模协变量的聚类二元数据，系统建立了GEE估计量的存在性、一致性与渐近正态性。Xie等[4]则进一步将GEE的渐近理论推广到每个受试者观测次数趋于无穷的情形，拓展了其理论的适用范围。

在高维纵向数据中，有效的特征选择能够降低模型复杂度与计算负担，有助于更深入地理解数据趋势，从而提高预测精度，并揭示时间动态性与个体异质性[5] [6]。常见的变量选择方法包括过滤法、包裹法和嵌入法[5] [7]，其中嵌入法能够将变量选择直接融入模型估计过程，在高维建模中得到广泛应用。纵向数据的特征选择通常基于GEE框架引入正则化项(如LASSO [8]、SCAD [9]或自适应LASSO [10])，以同步实现特征选择与参数估计[11]。在此基础上，Xu等[12]针对参数维数发散条件下的相关二元数据，引入收缩惩罚并系统分析了正则化GEE估计量的一致性与变量选择性质，为高维纵向二元数据的特征选择提供了重要的理论依据。

近年来，“大n，发散p”(large-n, diverging-p)框架下的高维理论研究受到广泛关注[13]，其主要考察协变量维度随样本量趋于无穷时的变量选择问题。Wang等[14]针对参数维数随样本量增长的广义线性模型，研究了Bridge惩罚下估计量一致性与变量选择一致性，为高维正则化方法在发散维度情形下的理论分析提供了重要参考。Wang等[6]研究了“大n，发散p”框架下的惩罚化GEE (penalized GEE, PGEE)，证明其具有良好的大样本渐近性质与Oracle性质[9] [13]，但在多重共线性情形下难以平衡稀疏性与稳定性。

为了解决多重共线性问题，Zou等[15]提出的Elastic net通过结合L₁与L₂惩罚促进相关变量的联合选择，随后SCAD-L2 [16]与自适应Elastic net [17]等扩展方法相继被提出，在保持估计稳定性的同时提升了变量选择一致性，使其更适用于强相关协变量下的高维稀疏建模。在此基础上，Blommaert等[18]在正态情形且协变量维数固定的设定下，通过引入非凸惩罚以应对纵向数据的共线性问题。此外，Lin等[19]提出的GPGEE将组惩罚机制引入GEE，可以实现已知分组信息时纵向数据的变量选择，但该方法无法主动识别纵向数据中隐含的特征关联，在实际应用中存在局限。

为此，本文首先将SCAD-L2正则化引入GEE框架，SCAD-L2正则化中SCAD惩罚能够在实现变量选择的同时减少对大系数的惩罚偏差，从而提高变量选择的准确性；L₂惩罚则增强了模型在多重共线性下的数值稳定性，有助于获得更为稳健的参数估计。因此，SCAD-L2惩罚能够在稀疏性与稳定性之间取得更好的平衡，特别适用于协变量高度相关的高维纵向数据。随后，本文在“大n，发散p”框架下，研究了该方法的参数估计一致性与Oracle性质。最后通过模拟研究验证了该方法在估计精度、变量选择准确性及稳定性等方面的综合优势。

2. 模型与方法

2.1. 基于SCAD-L2的高维GEE方法

2.1.1. 方法框架与理论基础

GEE是一种用于纵向数据的估计方法，其核心思想是通过刻画响应变量的边际均值函数，并结合工作相关矩阵近似描述组内相关结构，从而获得稳健且有效的参数估计。对于第$i$ 个个体的第$j$ 次观测，记响应变量为$Y_{ij}$ ，协变量为$X_{ij}$ ，其中$i=1,\cdots ,n$ ，$j=1,\cdots ,m_i$ 。这里，$m_i$ 表示个体$i$ 的重复测量次数，为简化理论推导，通常假设各个体的测量次数相同，即$m_i=m$ ；协变量的维数记为$p$ ；因此可定义个体$i$ 的响应向量为$Y_i={\left(Y_{i1},Y_{i2},\cdots ,Y_{im}\right)}^{\text{T}}$ ，协变量矩阵为$X_i={\left(X_{i1},X_{i2},\cdots ,X_{im}\right)}^{\text{T}}$ 。GEE具体表示为

$\begin{array}{c}S\left(\beta \right)=\frac{1}{n}{\displaystyle \sum _{i=1}^n{D}_i^{\text{T}}\left(\beta \right)V_i^{-1}\left(\beta \right)\left(Y_i-{\mu }_i\left(\beta \right)\right)}\\ =\frac{1}{n}{\displaystyle \sum _{i=1}^n{X}_i^{\text{T}}}{A}_i^{1/2}\left(\beta \right){\widehat{R}}^{-1}\left(\alpha \right)A_i{\left(\beta \right)}^{-1/2}\left(Y_i-{\mu }_i\left(\beta \right)\right)\\ =0,\end{array}$ (1)

其中$D_i\left(\beta \right)$ 表示个体均值函数对参数的导数矩阵，$V_i\left(\beta \right)$ 为个体内的工作协方差矩阵。在指数族分布的假设下$D_i\left(\beta \right)=A_i\left(\beta \right)X_i$ ，$V_i\left(\beta \right)=A_i^{1/2}\left(\beta \right)R\left(\alpha \right)A_i^{1/2}\left(\beta \right)$ ，$A_i\left(\beta \right)=\text{diag}\left\{{\sigma }_{i1}^2,{\sigma }_{i2}^2,\cdots ,{\sigma }_{im}^2\right\}$ 是一个由边际均值${\mu }_i\left(\beta \right)$ 决定的对角方差函数矩阵，$R\left(\alpha \right)$ 是由参数$\alpha$ 控制的工作相关矩阵，用于近似同一簇内的相关结构。

在高维情形下，协变量个数可能随着样本量的增加而同步增长，因此对估计方程的求解与数值实现提出更高要求。为此，Wang等[8]提出了PGEE，通过在估计方程中引入稀疏惩罚同时实现参数估计与变量选择。其统一形式为

$U\left(\beta \right)=S\left(\beta \right)-q_{\lambda }\left(\beta \right)=0$ ，

其中$q_{\lambda }(·)$ 为惩罚函数的一阶导数。

基于上述建模框架，本文进一步将SCAD-L2组合惩罚引入GEE框架，其估计方程形式为

$U\left(\beta \right)=S\left(\beta \right)-q_{{\lambda }_1,{\lambda }_2}\left(\beta \right)=0$ ， (2)

其中$S\left(\beta \right)$ 来自公式，惩罚函数导数部分为SCAD-L2惩罚函数的导数，形如

$q_{{\lambda }_1,{\lambda }_2}\left(\beta \right)={f^{\prime }}_{{\lambda }_1,SCAD}\left(\beta \right)+{f^{\prime }}_{{\lambda }_2,l_2}\left(\beta \right)$ ，

其中，

${f^{\prime }}_{{\lambda }_1,SCAD}\left(\theta \right)={\lambda }_1\left\{I\left(\theta \le {\lambda }_1\right)+\frac{{\left(a{\lambda }_1-\theta \right)}_{+}}{\left(a-1\right){\lambda }_1}I\left(\theta >{\lambda }_1\right)\right\}$ ，

${f^{\prime }}_{{\lambda }_2,l_2}\left(\theta \right)={\lambda }_2\theta$ ，

此处$\alpha$ ，${\lambda }_1$ ，${\lambda }_2$ 为调节参数，本文遵循Fan [13]的推荐取$\alpha =3.7$ ，${\lambda }_1$ ，${\lambda }_2$ 的选择通过网格搜索与交叉验证来确定。

2.1.2. 参数估计与迭代算法

为实现SCAD-L2惩罚下GEE的求解，本文遵循PGEE的算法框架，将牛顿–拉夫逊方法(Newton-Raphson method)与极小化–极大化(Majorization-Minimization)算法相结合[20]，形成一种稳定高效的混合求解策略。具体的参数更新方程为

${\widehat{\beta }}^k={\widehat{\beta }}^{k-1}+{\left[H\left({\widehat{\beta }}^{k-1}\right)+E\left({\widehat{\beta }}^{k-1}\right)\right]}^{-1}\left[S\left({\widehat{\beta }}^{k-1}\right)-E\left({\widehat{\beta }}^{k-1}\right)\right]$ ，

其中$H\left({\widehat{\beta }}^{k-1}\right)=\frac{1}{n}{\displaystyle \sum _{i=1}^n{X}_i^{\text{T}}}{A}_i^{1/2}\left({\widehat{\beta }}^{k-1}\right){\widehat{R}}^{-1}{A}_i^{1/2}\left({\widehat{\beta }}^{k-1}\right)X_i$为未加惩罚的GEE在当前参数估计值处的局部二阶近似，而对角矩阵$E\left({\widehat{\beta }}^{k-1}\right)=\mathrm{diag}\left\{\frac{q_{{\lambda }_1,{\lambda }_2}\left(\left|{\widehat{\beta }}_1\right|+\right)}{ϵ+\left|{\widehat{\beta }}_1\right|},\cdots ,\frac{q_{{\lambda }_1,{\lambda }_2}\left(\left|{\widehat{\beta }}_p\right|+\right)}{ϵ+\left|{\widehat{\beta }}_p\right|}\right\}$ 用于控制上界，此处$\left|{\widehat{\beta }}_j\right|+$ 表示略大于$\left|{\widehat{\beta }}_j\right|$ 的值，$ϵ>0$ 为任意给定的小正数常数，通常取10⁻⁶，该设计用于规避除以零的情形以保障数值稳定性。

2.1.3. 稳健性初始值选择

在迭代估计过程中，常用的策略是使用独立工作相关结构假设下的GEE估计值作为迭代初始值，但在多重共线性情形下信息矩阵接近奇异，导致该初始值方差膨胀甚至无法求解。为增强初始估计的数值稳定性并促进算法收敛，本文采用L₂正则化方法构建参数初始估计，该方法能够有效缓解共线性导致的不适定性，并为后续的非凸惩罚估计提供稳健的初始点。初始值估计方程形如

${\tilde{U}}_n\left({\beta }_n\right)=\frac{1}{n}{\displaystyle \sum _{i=1}^n{X}_i^{\text{T}}}\left(Y_i-\mu \left({\beta }_n\right)\right)-{\lambda }_{n2}{\beta }_n=0$ (3)

2.1.4. 计算开销与双参数调参复杂度

与仅含SCAD罚项的PGEE相比，SCAD-L2惩罚下GEE中叠加了L2正则化对应的Ridge项，该项等价于在更新方程中加入对角型的稳定化项或对参数向量施加线性收缩，其边际计算成本远小于构造与求解涉及$p$ 维参数的线性系统及工作相关结构更新等核心步骤，因此，在固定调参下模型拟合中，额外计算负担一般较为有限。

SCAD-L2惩罚下GEE的额外计算开销主要源于调参空间由一维扩展至二维所导致的模型重复拟合次数增加。对PGEE而言，惩罚选择通常集中在一维的${\lambda }_1$ 网格上，采用$B$ 折交叉验证的总体模型拟合次数约为$B\cdot L_{{\lambda }_1}$ 。但在同时包含稀疏化SCAD与稳定化L2两类惩罚强度时，需要在二维网格$\left({\lambda }_1,{\lambda }_2\right)$ 上进行搜索，使得总体拟合次数增加为$B\cdot L_{{\lambda }_1}\cdot L_{{\lambda }_2}$ ，该“乘法型”增长使得在$B$ 较大或网格较密时，总耗时会较PGEE明显上升。

综合来看，尽管SCAD-L2惩罚下GEE因采用基于L2惩罚GEE的初始值而增加了一次前置拟合开销，但在强多重共线性情形下，该初始值可显著提升数值稳定性，减少非凸惩罚迭代中的失败重启、震荡与额外迭代，从而以小幅前置代价换取更稳健的收敛。

2.2. 大样本渐近性质

本节在“大n，发散p”框架下研究SCAD-L2惩罚下GEE方法的变量选择与参数估计的渐近性质。记真实的系数向量为${\beta }_{n0}={\left({\beta }_{n10}^{\text{T}},{\beta }_{n20}^{\text{T}}\right)}^{\text{T}}$ ，其中${\beta }_{n20}=0$ ，${\beta }_{n10}$ 对应维度为$s_n$ 的非零系数。以下给出建立渐近性质所需的正则性条件：

(A1) 存在有限常数$C_X>0$ ，使得${\sup }_{i,j}\Vert X_{ij}\Vert \le C_X$ 。

(A2) 真实参数${\beta }_{n0}$ 严格位于紧集$ℬ\subset {ℝ}^{p_n}$ 内；估计量${\beta }_n$ 也被限制在$ℬ$ 中。

(A3) 存在固定常数$0<{\kappa }_1\le {\kappa }_2<\infty$ ，使得

${\kappa }_1\le {\lambda }_{\min }\left(n^{-1}{\displaystyle \sum _{i=1}^n{X}_i^{\text{T}}}{R}_0^{-1}{X}_i\right)+{\lambda }_{n2}\le {\lambda }_{\max }\left(n^{-1}{\displaystyle \sum _{i=1}^n{X}_i^{\text{T}}}{R}_0^{-1}{X}_i\right)+{\lambda }_{n2}\le {\kappa }_2$ 。

(A4) 真实的簇内相关矩阵$R_0$ 与固定的参考矩阵$\overline{R}$ 均为正定矩阵，且其特征值均被正数和有限常数所界定，即远离0和$+\infty$ 。此外，$\Vert {\widehat{R}}^{-1}-{\overline{R}}^{-1}\Vert =O_p\left(\sqrt{s_n/n}\right)$ ，该式表明估计矩阵${\widehat{R}}^{-1}$ 与参考矩阵${\overline{R}}^{-1}$ 的差距在概率意义下以$\sqrt{s_n/n}$ 的速率收敛。

(A5) 存在常数$c_1>0$ 和$\delta >0$ ，对所有$i$ ，有$E\left({\Vert {ϵ}_i\left({\beta }_{n0}\right)\Vert }^2\right)\le c_1$ ，其中${ϵ}_i\left(\beta \right)=A_i{\left(\beta \right)}^{-1/2}\left(Y_i-{\mu }_i\left(\beta \right)\right)$ 。

(A6) $s_n{\left(\log n\right)}^4=o\left(n{\lambda }_{n1}^2\right)$ ，$\log \left(p_n\right)=o\left(n{\lambda }_{n1}^2/{\left(\log n\right)}^2\right)$ 且$p_n{s}_n^4{\left(\log n\right)}^6=o\left(n^2{\lambda }_{n1}^2\right)$ 。

正则性条件(A1)~(A6)在很大程度上继承了PGEE的设定并做出了必要修正，这些条件共同保证了高维纵向数据下估计量选择一致性、估计一致性与渐近正态性。特别地，(A3)是为应对多重共线性引入的关键修正：与PGEE要求信息矩阵最小特征值有正下界不同，当协变量高度相关时，信息矩阵可能接近奇异，故而通过引入L2惩罚项${\lambda }_{n2}$ 作为扰动来保证理论分析中的稳定下界。此外，(A5)对标准化残差施加了矩条件，用于在高维设定下保证中心极限定理的适用性并控制随机误差项的尾部行为。

性质1：

假设满足正则化条件(A1)~(A3)，当$n\to \infty$ 时，满足$p_n^2/n\to 0$ 且${\lambda }_{n2}=O\left(n^{-1/2}\right)$ ，则初始值估计方程${\tilde{U}}_n\left({\beta }_n\right)=0$ (3)存在解${\tilde{\beta }}_n$ ，且该解满足：

$\Vert {\tilde{\beta }}_n-{\beta }_{n0}\Vert =O_p\left(\sqrt{p_n/n}+{\lambda }_{n2}{\beta }_{n0}\right)$ (4)

定理1：

假设满足正则化条件(A1)~(A5)，当$n\to \infty$ 时，满足$p_n^2/n\to 0$ ，${\lambda }_{n1}\to 0$ 且${\lambda }_{n2}=O\left(n^{-1}\right)$ ，则方程$U_n\left(\beta \right)=0$ (2)存在解${\widehat{\beta }}_n$ ，且该解满足：

$n^{-1}{p}_n^3=o\left(1\right)$ (5)

定理2：

假设满足正则化条件(A1)~(A6)，若$n^{-1}{p}_n^3=o\left(1\right)$ ，${\lambda }_{n1}\to 0$ 且${\lambda }_{n2}=O\left(n^{-1}\right)$ ，记${\widehat{\beta }}_n={\left({\widehat{\beta }}_{n1}^{\text{T}},{\widehat{\beta }}_{n2}^{\text{T}}\right)}^{\text{T}}$为方程$U_n\left({\beta }_n\right)=0$ 的解，以下结论成立：

(1) $P\left(\left|U_{nj}\left({\widehat{\beta }}_n\right)\right|=0,j=1,\cdots ,s_n\right)\to 1$ ，$P\left(\left|U_{nj}\left({\widehat{\beta }}_n\right)\right|\le \frac{{\lambda }_{n1}}{\log n},j=s_n+1,\cdots ,p_n\right)\to 1$ ；

(2) $P\left({\widehat{\beta }}_{n2}=0\right)\to 1$ ；

(3) 对于任意满足$\Vert {\alpha }_n\Vert =1$ 的${\alpha }_n\in {ℝ}^{s_n}$ ，有${\alpha }_n^{\text{T}}{\overline{M}}_{n1}^{-1/2}\left({\beta }_{n0}\right)\left({\overline{H}}_{n1}\left({\beta }_{n0}\right)+{\lambda }_{n2}{I}_{n1}\right)\left({\widehat{\beta }}_{n1}-{\beta }_{n10}\right)~N\left(0,1\right)$，

其中，

${\overline{M}}_{n1}\left({\beta }_{n0}\right)=\frac{1}{n^2}{\displaystyle \sum _{i=1}^n{X}_{i1}^{\text{T}}}{A}_i^{1/2}\left({\beta }_{n0}\right){\overline{R}}^{-1}{R}_0{\overline{R}}^{-1}{A}_i^{1/2}\left({\beta }_{n0}\right)X_{i1}+{\lambda }_{n2}{I}_{n1},$

${\overline{H}}_{n1}\left({\beta }_{n0}\right)=\frac{1}{n}{\displaystyle \sum _{i=1}^n{X}_{i1}^{\text{T}}}{A}_i^{1/2}\left({\beta }_{n0}\right)R_0^{-1}{A}_i^{1/2}\left({\beta }_{n0}\right)X_{i1}.$

注：定理2确立了估计量的Oracle性质，即能够一致地识别零系数、精确估计非零系数，并实现非零系数的渐近正态性。${\overline{M}}_{n1}$ 表示估计方程的渐近协方差矩阵，用于刻画得分函数的变异性。在理论推导过程中施加了若干条件，包括对惩罚参数${\lambda }_{n2}$ 较为严格的速率条件，这一设定旨在有效控制L₂惩罚引入的偏差，并为渐近正态性的成立提供支持。在实际应用中，限制条件可以适当放宽。

性质1的证明：

证明式(4)等价于证明：对于任何$ϵ>0$ ，存在一个$n$ ，使得对于所有的$\Delta >0$ 均有

$Pr\left\{\underset{\Vert {\beta }_n-{\beta }_{n0}\Vert =\Delta \left(\sqrt{\frac{p_n}{n}}+\left|{\lambda }_{n2}{\beta }_{n0}\right|\right)}{\sup }{\left({\widehat{\beta }}_n-{\beta }_{n0}\right)}^{\text{T}}{\tilde{U}}_n\left({\beta }_n\right)<0\right\}>1-ϵ$，

通过应用泰勒展开，对于介于${\beta }_{n0}$ 和${\beta }_n$ 之间的${\beta }_n^{\ast }$ ，得到：

$\begin{array}{c}{\left({\tilde{\beta }}_n-{\beta }_{n0}\right)}^{\text{T}}{\tilde{U}}_n\left({\beta }_n\right)={\left({\beta }_n-{\beta }_{n0}\right)}^{\text{T}}{\tilde{U}}_n\left({\beta }_{n0}\right)+{\left({\beta }_n-{\beta }_{n0}\right)}^{\text{T}}\frac{\partial }{\partial {\beta }_n}{U}_n\left({\beta }_n^{*}\right)\left({\beta }_n-{\beta }_{n0}\right)\\ =I_{n1}+I_{n2}.\end{array}$ (6)

对于$I_{n1}$ 项，有：

$\begin{array}{l}{I}_{n1}={\left({\beta }_n-{\beta }_{n0}\right)}^{\text{T}}{\tilde{U}}_n\left({\beta }_{n0}\right)\\ \le \Vert {\beta }_n-{\beta }_{n0}\Vert \cdot \Vert \frac{1}{n}{\displaystyle \sum _{i=1}^n{X}_i^{\text{T}}}\left(Y_i-{\mu }_i\left({\beta }_{n0}\right)\right)-{\lambda }_{n2}{\beta }_{n0}\Vert \\ \le \Delta \left(\sqrt{\frac{p_n}{n}}+\Vert {\lambda }_{n2}{\beta }_{n0}\Vert \right)\cdot \left(O_p\left(\sqrt{\frac{p_n}{n}}\right)+\Vert {\lambda }_{n2}{\beta }_{n0}\Vert \right)\\ =\Delta O_p\left({\gamma }_n^2\right),\end{array}$

其中${\gamma }_n\triangleq \sqrt{\frac{p_n}{n}}+\Vert {\lambda }_{n2}{\beta }_{n0}\Vert$ 。第二个不等式成立源于(A1)下能保证$E\left({\Vert \frac{1}{n}{\displaystyle \sum _{i=1}^n{X}_i^{\text{T}}}\left(Y_i-{\mu }_i\left({\beta }_{n0}\right)\right)\Vert }^2\right)\le O_p\left(\frac{p_n}{n}\right)$ 。

对于$I_{n2}$ 项，进一步拆分为$I_{n21}$ 和$I_{n22}$ ：

$\begin{array}{c}{I}_{n2}={\left({\beta }_n-{\beta }_{n0}\right)}^{\text{T}}\frac{\partial }{\partial {\beta }_n}{\tilde{U}}_n\left({\beta }_n^{*}\right)\left({\beta }_n-{\beta }_{n0}\right)\\ ={\left({\beta }_n-{\beta }_{n0}\right)}^{\text{T}}\left(-\frac{1}{n}{\displaystyle \sum _{i=1}^n{X}_i^{\text{T}}}{A}_i\left({\beta }_n^{*}\right)X_i-{\lambda }_{n2}{I}_n\right)\left({\beta }_n-{\beta }_{n0}\right)\\ ={\left({\beta }_n-{\beta }_{n0}\right)}^{\text{T}}\left(-\frac{1}{n}{\displaystyle \sum _{i=1}^n{X}_i^{\text{T}}}{A}_i\left({\beta }_n\right)X_i-{\lambda }_{n2}{I}_n\right)\left({\beta }_n-{\beta }_{n0}\right)\\ +{\left({\beta }_n-{\beta }_{n0}\right)}^{\text{T}}\left[-\frac{1}{n}{\displaystyle \sum _{i=1}^n{X}_i^{\text{T}}}\left(A_i^{*}\left({\beta }_n\right)-A_i\left({\beta }_n\right)\right)X_i\right]\left({\beta }_n-{\beta }_{n0}\right)\\ =I_{n21}+I_{n22}.\end{array}$

对于$I_{n22}$ 项，由$p^2{n}^{-1}=o\left(1\right)$ ，${\lambda }_{n2}^2{p}_n=o\left(1\right)$ 以及(A3)有：

$\begin{array}{l}{I}_{n22}\le {\Vert {\beta }_n-{\beta }_{n0}\Vert }^2\frac{1}{n}{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^m\left(A_i\left({\beta }_n\right)-A_i^{*}\left({\beta }_n\right)\right)X_{ij}{X}_{ij}^{\text{T}}}}\\ \le {\Vert {\beta }_n-{\beta }_{n0}\Vert }^2\frac{1}{n}\underset{i,j}{\sup }\Vert X_{ij}\Vert \cdot \Vert {\beta }_n-{\beta }_n^{*}\Vert {\lambda }_{max}\left({\displaystyle \sum X_i^{\text{T}}}{X}_i\right)\\ \le {\Delta }^2{\left(\sqrt{p_n/n}+\Vert {\lambda }_{n2}{\beta }_{n0}\Vert \right)}^2\frac{1}{n}{O}_p\left(p_n\right)O_p\left(\sqrt{p_n/n}+\Vert {\lambda }_{n2}{\beta }_{n0}\Vert \right)O_p\left(n\right)\\ ={\Delta }^2{O}_p\left(p_n^2/n^{3/2}\right)\\ ={\Delta }^2{o}_p\left(p_n/n\right).\end{array}$

接下来分析$I_{n21}$ 项：

$\begin{array}{l}{I}_{n21}=-{\left({\beta }_n-{\beta }_{n0}\right)}^{\text{T}}\left(\frac{1}{n}{\displaystyle \sum _{i=1}^n{X}_i^{\text{T}}}{A}_i\left({\beta }_n\right)X_i+{\lambda }_{n2}{I}_n\right)\left({\beta }_n-{\beta }_{n0}\right)\\ \le -{\Vert {\beta }_n-{\beta }_{n0}\Vert }^2\left({\lambda }_{min}\left(A_i\left({\beta }_n\right)\right){\lambda }_{min}\left(\frac{1}{n}{\displaystyle \sum _{i=1}^n{X}_i^{\text{T}}{X}_i}\right)+{\lambda }_{n2}\right)\\ \le -{\Vert {\beta }_n-{\beta }_{n0}\Vert }^2\left({\lambda }_{min}\left(A_i\left({\beta }_n\right)\right){\lambda }_{min}\left(R_0\right){\lambda }_{min}\left(\frac{1}{n}{\displaystyle \sum _{i=1}^n{X}_i^{\text{T}}{R}_0^{-1}{X}_i}\right)+{\lambda }_{n2}\right)\\ \le -{\Delta }^2{\gamma }_n^2\cdot \left(\left(a_A{a}_R{\lambda }_{min}\left(\frac{1}{n}{\displaystyle \sum _{i=1}^n{X}_i^{\text{T}}{R}_0^{-1}{X}_i}\right)+{\lambda }_{n2}\right)+{\lambda }_{n2}\left(1-a_A{a}_R\right)\right)\\ \le -C\cdot {\Delta }^2{\gamma }_n^2,\end{array}$

其中$a_A$ 与$a_R$ 分别表示${\lambda }_{min}\left(A_i\right)\left({\beta }_n\right)$ 与${\lambda }_{min}\left(R_0\right)$ ，且$0<a_A{a}_R<1$ 由标准化数据的协方差矩阵保证。

最后，泰勒展开式(6)由$I_{n21}$ 项支配，即对于所有充分大的$\Delta$ ，$I_{n21}$ 项具有较大的负值，由此性质1得证。

定理1的证明：

证明式(5)等价于证明：对于任何$ϵ>0$ ，存在一个$n$ ，使得对于所有的$\Delta >0$ ，均有

$Pr\left\{\underset{\Vert {\beta }_n-{\beta }_{n0}\Vert =\Delta \left(\sqrt{\frac{p_n}{n}}+\left|{\lambda }_{n2}{\beta }_{n0}\right|\right)}{\sup }{\left({\widehat{\beta }}_n-{\beta }_{n0}\right)}^{\text{T}}{\widehat{U}}_n\left({\beta }_n\right)<0\right\}>1-ϵ$，

通过应用泰勒展开，对于介于${\beta }_{n0}$ 和${\beta }_n$ 之间的${\beta }_n^{*}$ ，得到：

$\begin{array}{c}{\left({\beta }_n-{\beta }_{n0}\right)}^{\text{T}}{U}_n\left({\beta }_n\right)={\left({\beta }_n-{\beta }_{n0}\right)}^{\text{T}}{U}_n\left({\beta }_{n0}\right)+{\left({\beta }_n-{\beta }_{n0}\right)}^{\text{T}}\frac{\partial }{\partial {\beta }_n}{U}_n\left({\beta }_n^{*}\right)\left({\beta }_n-{\beta }_{n0}\right)\\ =I_{n1}+I_{n2},\end{array}$ (7)

其中，$I_{n1}={\left({\beta }_n-{\beta }_n\right)}^{\text{T}}{\overline{U}}_n\left({\beta }_n\right)+{\left({\beta }_n-{\beta }_{n0}\right)}^{\text{T}}\left(U_n\left({\beta }_{n0}\right)-{\overline{U}}_n\left({\beta }_n\right)\right)=I_{n11}+I_{n12}$。

对于$I_{n11}$ 项，有：

$\begin{array}{l}{I}_{n11}={\left({\beta }_n-{\beta }_{n0}\right)}^{\text{T}}\left({\overline{S}}_n\left({\beta }_n\right)-q_{scad,{\lambda }_{n1}}\left({\beta }_{n0}\right)-{\lambda }_{n2}{\beta }_{n0}\right)\\ \le \Vert {\beta }_{n0}-{\beta }_n\Vert \cdot \left(\Vert \overline{S_n}{\left(\beta \right)}_{n0}\Vert +\Vert {\lambda }_{n2}{\beta }_{n0}\Vert \right)\\ \le \Delta \left(\sqrt{\frac{p_n}{n}}+\Vert {\lambda }_{n2}{\beta }_{n0}\Vert \right)\cdot O_p\left(\sqrt{\frac{p_n}{n}}+\left|{\lambda }_{n2}{\beta }_{n0}\right|\right)\\ =\Delta {\gamma }_n{O}_p\left({\gamma }_n\right),\end{array}$

其中${\gamma }_n=\sqrt{p_n/n}+\Vert {\lambda }_{n2}{\beta }_{n10}\Vert$ 且$E\left({\Vert \overline{S}\left({\beta }_{n0}\right)\Vert }^2\right)\le O_p\left(\frac{p_n}{n}\right)$ 。由SCAD惩罚的性质可知，由于在$I_{n12}$ 时有${\lambda }_{n1}\to 0$ ，对于足够大的$I_{n12}$ ，当$\theta >a{\lambda }_{n1}$ 时，$I_{n12}$ 与$I_{n12}$ 成立，因此惩罚项对真实非零参数的影响可以忽略。

对于$I_{n12}$ 项，由Wang的引理3.1 [3]可知$\Vert U_n\left({\beta }_{n0}\right)-{\overline{U}}_n\left({\beta }_n\right)\Vert =O\left(\frac{p_n}{n}\right)$ ，因此$I_{n12}\le \Delta {\gamma }_{n}O\left(\frac{p_n}{n}\right)$ 。

由于$-\frac{\partial }{\partial {\beta }_n}{U}_n\left({\beta }_n^{*}\right)\triangleq D_n\left({\beta }_n^{*}\right)+{q^{\prime }}_{scad,{\lambda }_{n1}}\left({\beta }_n^{*}\right)+{\lambda }_{n2}I$ ，其中$D_n\left({\beta }_n\right)=-\frac{\partial }{\partial {\beta }_n}{S}_n\left({\beta }_n\right)$ ，因此，对于$I_{n2}$ 项有：

$\begin{array}{l}{I}_{n2}=-{\left({\beta }_n-{\beta }_{n0}\right)}^{\text{T}}\left({\overline{D}}_n\left({\beta }_n^{*}\right)+{q^{\prime }}_{scad,{\lambda }_{n1}}\left({\beta }_n^{*}\right)+{\lambda }_{n2}{I}_n\right)\left({\beta }_n-{\beta }_{n0}\right)\\ -{\left({\beta }_n-{\beta }_{n0}\right)}^{\text{T}}\left(D_n\left({\beta }_n^{*}{\beta }_{n0}\right)-{\overline{D}}_n\left({\beta }_n^{*}\right)\right)\left({\beta }_n-{\beta }_{n0}\right)\\ \triangleq I_{n21}+I_{n22},\end{array}$

由前述可知${q^{\prime }}_{scad,{\lambda }_{n1}}\left({\beta }_n^{*}\right)=0$ ，因此，

$\begin{array}{c}{I}_{n21}=-{\left({\beta }_n-{\beta }_{n0}\right)}^T\left({\overline{H}}_n\left({\beta }_{n0}\right)+{\lambda }_{n2}{I}_n\right)\left({\beta }_n-{\beta }_{n0}\right)-{\left({\beta }_n-{\beta }_{n0}\right)}^{\text{T}}\left({\overline{H}}_n\left({\beta }_n^{*}\right)-{\overline{H}}_n\left({\beta }_{n0}\right)\right)\left({\beta }_n-{\beta }_{n0}\right)\\ -{\left({\beta }_n-{\beta }_{n0}\right)}^T\left({\overline{D}}_n\left({\beta }_n^{*}\right)-{\overline{H}}_n\left({\beta }_n^{*}\right)\right)\left({\beta }_n-{\beta }_{n0}\right)\\ \triangleq I_{n21}^a+I_{n21}^b+I_{n21}^c,\end{array}$

易证$I_{n21}^b={\Delta }^2{o}_p\left(p_n/n\right)$ 且$I_{n21}^c={\Delta }^2{o}_p\left(p_n/n\right)$ 。结合${\overline{H}}_n\left({\beta }_{n0}\right)=\frac{1}{n}{\displaystyle \sum _{i=1}^n{X}_i^{\text{T}}}{A}_i^{1/2}\left({\beta }_{n0}\right){\overline{R}}^{-1}{A}_i^{1/2}\left({\beta }_{n0}\right)X_i$与正则化条件(A1)~(A4)，对于$I_{n21}^a$ ，可以得到

$\begin{array}{c}{I}_{n21}^a=-{\left({\beta }_n-{\beta }_{n0}\right)}^{\text{T}}\left(\frac{1}{n}{\displaystyle \sum _{i=1}^n{X}_i^{\text{T}}}{A}_i^{1/2}\left({\beta }_{n0}\right)R_0^{-1}{A}_i^{1/2}{X}_i+{\lambda }_{n2}{I}_n\right)\left({\beta }_n-{\beta }_{n0}\right)\\ -{\left({\beta }_n-{\beta }_{n0}\right)}^{\text{T}}\left(\frac{1}{n}{\displaystyle \sum _{i=1}^n{X}_i^{\text{T}}}{A}_i^{1/2}\left({\beta }_{n0}\right)\left[{\overline{R}}^{-1}-R_0^{-1}\right]A_i^{1/2}{X}_i\right)\left({\beta }_n-{\beta }_{n0}\right)\\ \le -{\Vert {\beta }_n-{\beta }_{n0}\Vert }^2\left[{\lambda }_{min}\left(A_i\left({\beta }_{n0}\right)\right)\cdot \frac{1}{n}{\displaystyle \sum _{i=1}^n{X}_i^{\text{T}}}{R}_0^{-1}{X}_i+{\lambda }_{n2}\right]\\ -{\Vert {\beta }_n-{\beta }_{n0}\Vert }^2\left[{\lambda }_{min}\left(\left[{\overline{R}}^{-1}-R_0^{-1}\right]\right)\cdot {\lambda }_{min}\left(A_i\left({\beta }_{n0}\right)\right)\cdot {\lambda }_{min}\left(\frac{1}{n}{\displaystyle \sum _{i=1}^n{X}_i^{\text{T}}}{X}_i\right)\right]\\ \le -{\Vert {\beta }_n-{\beta }_{n0}\Vert }^2\left[a_A{\lambda }_{min}\left(\frac{1}{n}{\displaystyle \sum _{i=1}^n{X}_i^{\text{T}}}{R}_0^{-1}{X}_i\right)+{\lambda }_{n2}\right]-{\Vert {\beta }_n-{\beta }_{n0}\Vert }^2{C}_2\\ \le -{\Vert {\beta }_n-{\beta }_{n0}\Vert }^2\left(C_1+C_2\right)\\ \le -C{\Delta }^2{\gamma }_n^2,\end{array}$

其中$a_A={\lambda }_{min}\left(A_i\left({\beta }_n\right)\right)$ ，$C_1>0$ 且$C_2\ge 0$ ，$C_1=\min \left(a_A,{\lambda }_{n2}\right)$ 。

最后，泰勒展开式(7)由$I_{n21}^a$ 支配，对于所有充分大的$\Delta$ ，$I_{n21}^a$ 是大且负的，因此定理1得证。

定理2的证明：

对于定理2中的(1)与(2)，Wang等[8]已在SCAD惩罚下的GEE框架中给出了完整的理论证明，其所依赖的正则条件与本文研究内容无冲突，为避免赘述，此处不再重复证明。下文将重点证明性质(3)。

首先证明${\alpha }_{n1}^{\text{T}}{\overline{M}}_{n1}^{-1/2}\left({\beta }_{n10}\right){\overline{U}}_{n1}\left({\beta }_{n10}\right)~N\left(0,1\right)$，现有：

$\begin{array}{l}{\alpha }_{n1}^{\text{T}}{\overline{M}}_{n1}^{-1/2}\left({\beta }_{n10}\right){\overline{U}}_{n1}\left({\beta }_{n10}\right)\\ ={\alpha }_{n1}^{\text{T}}{\overline{M}}_{n1}^{-1/2}\left\{\frac{1}{n}{\displaystyle \sum _{i=1}^n{X}_i^{\text{T}}}{A}_i^{1/2}\left({\beta }_{n10}\right){\overline{R}}^{-1}{A}_i^{-1/2}\left({\beta }_{n10}\right)\left(Y_i-{\mu }_i\left({\beta }_{n10}\right)\right)+{\lambda }_{n2}{\beta }_{n10}\right\}\\ ={\displaystyle \sum _{i=1}^n{Z}_{ni}},\end{array}$

定义$Z_{ni}=\frac{1}{n}{\alpha }_{n1}^{\text{T}}{\overline{M}}_{n1}^{-1/2}\left({\beta }_{n10}\right)X_i^{\text{T}}{A}_i^{1/2}\left({\beta }_{n10}\right){\overline{R}}^{-1}{\epsilon }_i\left({\beta }_{n10}\right)+{\lambda }_{n2}{\beta }_{n10}\triangleq I_{n1}+I_{n2}$，再根据柯西–施瓦茨不等式，有：

$Z_{ni}^2={\Vert I_{n1}+I_{n2}\Vert }^2\le {\Vert I_{n1}\Vert }^2+{\Vert I_{n2}\Vert }^2+2\left|I_{n1}\cdot I_{n2}\right|\le 2\left({\Vert I_{n1}\Vert }^2+{\Vert I_{n2}\Vert }^2\right),$

由于${\Vert I_{n2}\Vert }^2={\lambda }_{n2}^2{\Vert {\beta }_{n10}\Vert }^2=O\left(1/n^2\right)$ ，且

$\begin{array}{c}{\Vert I_{n1}\Vert }^2\le \frac{1}{n^2}{\alpha }_{n1}^{\text{T}}{\overline{M}}_{n1}^{-1/2}\left({\beta }_{n10}\right)X_i^{\text{T}}{A}_i^{1/2}\left({\beta }_{n10}\right){\overline{R}}^{-1}{\Vert {\epsilon }_i\left({\beta }_{n10}\right)\Vert }^2{\overline{R}}^{-1}{A}_i^{1/2}\left({\beta }_{n10}\right)X_i{\overline{M}}_{n1}^{-1/2}\left({\beta }_{n10}\right){\alpha }_{n1}\\ \le {\Vert {\epsilon }_i\left({\beta }_{n10}\right)\Vert }^2{\lambda }_{max}\left(A_i\left({\beta }_{n10}\right)\right){\lambda }_{max}\left({\overline{R}}^{-1}\right)\cdot \frac{1}{n^2}{\alpha }_{n1}^{\text{T}}{\overline{M}}_{n1}^{-1/2}\left({\beta }_{n10}\right)X_i^{\text{T}}{\overline{R}}^{-1}{X}_i{\overline{M}}_{n1}^{-1}\left({\beta }_{n10}\right){\alpha }_{n1}\\ =C\cdot {\psi }_{ni}\cdot {\Vert {\epsilon }_i\left({\beta }_{n10}\right)\Vert }^2.\end{array}$

定义其中${\psi }_{ni}\triangleq \frac{1}{n^2}{\alpha }_{n1}^{\text{T}}{\overline{M}}_{n1}^{-1/2}\left({\beta }_{n10}\right)X_i^{\text{T}}{\overline{R}}^{-1}{X}_i{\overline{M}}_{n1}^{-1/2}\left({\beta }_{n10}\right){\alpha }_{n1}$且满足${\psi }_{ni}\le \frac{1}{n^2}{\lambda }_{max}\left(X_i^{\text{T}}{\overline{R}}^{-1}{X}_i\right)\cdot {\lambda }_{max}\left({\overline{M}}_{n1}^{-1}\left({\beta }_{n10}\right)\right)=\frac{1}{n^2}\frac{{\lambda }_{max}\left(X_i^{\text{T}}{\overline{R}}^{-1}{X}_i\right)}{{\lambda }_{min}\left({\overline{M}}_{n1}\left({\beta }_{n10}\right)\right)}$。

进一步，矩阵${\overline{M}}_{n1}\left({\beta }_{n10}\right)$ 的最小特征值可以估计为

$\begin{array}{c}{b}_{n1}^{\text{T}}{\overline{M}}_{n1}\left({\beta }_{n10}\right)b_{n1}=b_{n1}^{\text{T}}\left(\frac{1}{n^2}{\displaystyle \sum _{i=1}^n{X}_i^{\text{T}}}{A}_i^{-1/2}\left({\beta }_{n10}\right){\overline{R}}^{-1}{R}_0{\overline{R}}^{-1}{A}_i^{1/2}\left({\beta }_{n10}\right)X_i+{\lambda }_{n2}{I}_{n1}\right)b_{n1}\\ \ge b_{n1}^{\text{T}}\left({\min }_i{\lambda }_{min}\left(A_i\left({\beta }_{n10}\right)\right)\frac{1}{n^2}{\displaystyle \sum _{i=1}^n{X}_i^{\text{T}}}{R}_0^{-1}{X}_i+{\lambda }_{n2}{I}_{n1}\right)b_{n1}\\ \ge \frac{c}{n}{\lambda }_{min}\left(\frac{1}{n}{\displaystyle \sum _{i=1}^n{X}_i^{\text{T}}}{R}_0^{-1}{X}_i\right)+{\lambda }_{n2}\end{array}$

由正则化假设(A4)，可得$\Vert {\overline{R}}^{-1}{R}_0{\overline{R}}^{-1}-R_0^{-1}\Vert =o\left(1\right)$ 。因此，

${\psi }_{ni}=\frac{1}{n^2}\frac{{\lambda }_{max}\left(X_i^{\text{T}}{\overline{R}}^{-1}{X}_i\right)}{\frac{c}{n}{\lambda }_{min}\left(\frac{1}{n}{\displaystyle \sum _{i=1}^n{X}_i^{\text{T}}}{R}_0^{-1}{X}_i\right)+{\lambda }_{n2}}=O\left(p_n/n\right),$

最终可以得到

$\begin{array}{c}{Z}_{ni}^2\le 2\left({\Vert I_{n1}\Vert }^2+{\Vert I_{n2}\Vert }^2\right)\\ =\text{2}\left(C\cdot {\psi }_{ni}\cdot {\Vert {\epsilon }_i\left({\beta }_{n10}\right)\Vert }^2+o\left(1\right)\right)\\ =C\cdot {\psi }_{ni}\cdot {\Vert {\epsilon }_i\left({\beta }_{n10}\right)\Vert }^2.\end{array}$

由(A5)中$E\left({\Vert {\epsilon }_{ni}\Vert }^2\right)=c_1$ ，根据Lindeberg条件，有

${\displaystyle \sum _{i=1}^{n}E}\left[Z_{ni}^{2}I\left(\left|Z_{ni}\right|>\epsilon \right)\right]\le {\displaystyle \sum _{i=1}^{n}C{\psi }_{ni}E}\left[{\Vert {\epsilon }_i\left({\beta }_{n0}\right)\Vert }^{2}I\left({\Vert {\epsilon }_i\left({\beta }_{n0}\right)\Vert }^2>\frac{{\epsilon }^2}{C{\psi }_{ni}}\right)\right].$

因此${\alpha }_{n1}^{\text{T}}{\overline{M}}_{n1}^{-1/2}\left({\beta }_{n10}\right){\overline{U}}_{n1}\left({\beta }_{n10}\right)~N\left(0,1\right)$得证。最后，重新整理原式有

$\begin{array}{l}{\alpha }_{n1}^{\text{T}}{\overline{M}}_{n1}^{-1/2}\left({\beta }_{n10}\right){\overline{U}}_{n1}\left({\beta }_{n10}\right)\\ ={\alpha }_{n1}^{\text{T}}{\overline{M}}_{n1}^{-1/2}\left({\beta }_{n10}\right)U_{n1}\left({\beta }_{n10}\right)+{\alpha }_{n1}^{\text{T}}{\overline{M}}_{n1}^{-1/2}\left({\beta }_{n10}\right)\left[{\overline{U}}_{n1}\left({\beta }_{n10}\right)-U_{n1}\left({\beta }_{n10}\right)\right]\\ ={\alpha }_{n1}^{\text{T}}{\overline{M}}_{n1}^{-1/2}\left({\beta }_{n10}\right)\left[D_{n1}\left({\beta }_{n1}^{*}\right)+f_{{{\lambda }^{\prime }}_{n1}}\left({\beta }_{n10}\right)+{\lambda }_{n2}{I}_{n1}\right]\left({\widehat{\beta }}_{n1}-{\beta }_{n10}\right)+{\alpha }_{n1}^{\text{T}}{\overline{M}}_{n1}^{-1/2}\left({\beta }_{n10}\right)\left[{\overline{U}}_{n1}\left({\beta }_{n10}\right)-U_{n1}\left({\beta }_{n10}\right)\right]\\ ={\alpha }_{n1}^{\text{T}}{\overline{M}}_{n1}^{-1/2}\left({\beta }_{n10}\right)\left[{\overline{H}}_{n1}\left({\beta }_{n10}\right)+f_{{{\lambda }^{\prime }}_{n1}}\left({\beta }_{n10}\right)+{\lambda }_{n2}{I}_{n1}\right]\left({\widehat{\beta }}_{n1}-{\beta }_{n10}\right)\\ +{\alpha }_{n1}^{\text{T}}{\overline{M}}_{n1}^{-1/2}\left({\beta }_{n10}\right)\left[D_{n1}\left({\beta }_{n1}^{*}\right)-{\overline{H}}_{n1}\left({\beta }_{n10}\right)\right]\left({\widehat{\beta }}_{n1}-{\beta }_{n10}\right)+{\alpha }_{n1}^{\text{T}}{\overline{M}}_{n1}^{-1/2}\left({\beta }_{n10}\right)\left[{\overline{U}}_{n1}\left({\beta }_{n10}\right)-U_{n1}\left({\beta }_{n10}\right)\right].\end{array}$

文献[3]指出，

$\underset{\Vert {\beta }_n-{\beta }_{n10}\Vert \le \Delta {\gamma }_n}{\sup }\left|{\alpha }_{n1}^T{\overline{M}}_{n1}^{-1/2}\left({\beta }_{n10}\right)\left(D_{n1}\left({\beta }_{n1}^{*}\right)-{\overline{H}}_{n1}\left({\beta }_{n10}\right)\right)\left({\widehat{\beta }}_n-{\beta }_{n10}\right)\right|=o_p\left(1\right),$

$\underset{\Vert {\beta }_n-{\beta }_{n10}\Vert \le \Delta {\gamma }_n}{\sup }\left|{\alpha }_{n1}^{\text{T}}{\overline{M}}_{n1}^{-1/2}\left({\beta }_{n10}\right)\left({\overline{U}}_{n1}\left({\beta }_{n10}\right)-U_{n1}\left({\beta }_{n10}\right)\right)\right|=o_p\left(1\right),$

最终可以得到：

$n=200,\rho =0.5$

至此定理2证毕。

3. 模拟实验

本节通过广泛的模拟研究评估并比较广义估计方程(GEE)、惩罚广义估计方程(PGEE)及SCAD-L2惩罚GEE (简称SLGEE)的性能表现，惩罚系数${\lambda }_{n1}$ 与${\lambda }_{n2}$ 的选取通过网格搜索结合交叉验证实现。在模型拟合中，为考察工作相关设定对估计与变量选择的影响，对GEE、PGEE和SLGEE分别采用三种工作相关结构进行计算：独立结构(indep)、可交换结构(exch)以及一阶自回归结构(AR-1)，每次模拟可得到“方法 × 工作相关结构”的多组结果，用于衡量不同工作相关矩阵假设下模型估计与变量选择的稳定性。由于MM算法固有的数值扰动，所估系数不会精确为零，因此将绝对值小于0.001的系数判定为非活跃变量，并将其从最终模型中剔除。

实验设计涵盖两类建模情形：其一为高斯型响应，采用恒等连接函数；其二为二分类响应，采用Logistic连接函数。对模拟实验结果采用多种指标评估模型性能：均方误差(MSE)，计算公式为$\text{MSE}={\Vert {\widehat{\beta }}_n-{\beta }_{n0}\Vert }^2/p$；过拟合数，即真实零系数被错误选中的数量；欠拟合数，即真实非零系数被错误排除的数量；识别准确率，即正确识别系数的比例；选择稳定性，通过Jaccard指数量化，即重复模拟中非零系数集交集与并集的比值。

3.1. 实验1

相关正态响应数据的生成模型为

$Y_{ij}=X_{ij}^{\text{T}}\beta +{\epsilon }_{ij},i=1,\cdots ,n;j=1,\cdots ,m$ ，

令$m=4$ ，则$X_{ij}={\left(x_{ij,1},\cdots ,x_{ij,200}\right)}^{\text{T}}$ 为包含200个协变量的向量，真实参数设为

${\beta }_{n0}={\left(\underset{15}{\underbrace{3,\cdots ,3}},\underset{15}{\underbrace{1.5,\cdots ,1.5}},\underset{15}{\underbrace{2,\cdots ,2}},\underset{155}{\underbrace{0,\cdots ,0}}\right)}^{\text{T}},$

即前45个系数非零，其余155个为零。随机误差${\epsilon }_i$ 来自多元正态分布，均值为0，协方差矩阵为$EX\left(\rho \right)$ ，即对角元为1、非对角元均为$\rho$ 的对称矩阵。实验设置$\rho =0.5$ 和$\rho =0.8$ 两种情形，以反映不同程度的组内相关性。对任意$i$ ，协变量的生成方式为

$\begin{array}{l}{X}_{jk}=Z_{1,j}+{\epsilon }_{jk}^x,Z_{1,j}~N\left(0,EX\left(\rho \right)\right),j=1,\cdots ,4;k=1,\cdots ,15,\\ X_{jk}=Z_{2,j}+{\epsilon }_{jk}^x,Z_{2,j}~N\left(0,EX\left(\rho \right)\right),j=1,\cdots ,4;k=16,\cdots ,30,\\ X_{jk}=Z_{3,j}+{\epsilon }_{jk}^x,Z_{3,j}~N\left(0,EX\left(\rho \right)\right),j=1,\cdots ,4;k=31,\cdots ,45,\\ X_{jk}~N\left(0,EX\left(\rho \right)\right),j=1,\cdots ,4;k=46,\cdots ,200,\end{array}$

其中${\epsilon }_{jk}^x{~}_{i.i.d.}N\left(0,0.001\right),k=1,\cdots ,100$ 。

为了更直观地看出不同参数组合对模型的影响，图1展示了实验1中通过二维网格搜索与五折交叉验证选取惩罚参数${\lambda }_1,{\lambda }_2$ 的热力图。对每一个网格点均计算平均交叉验证误差，并选择使其最小的参数组合作为最优惩罚参数。结果显示，交叉验证误差关于${\lambda }_1$ 呈现明显的“先下降后上升”的趋势：当${\lambda }_1$ 过小时惩罚不足易导致过拟合，CV误差较大；当${\lambda }_1$ 过大时，真实信号被过度收缩从而产生欠拟合，CV误差亦上升。对${\lambda }_2$ 而言，小幅L2惩罚有助于提高估计稳定性并改善预测误差，但当${\lambda }_2$ 增大至0.1或1时会引发过度收缩，导致CV误差显著增大。最终在所设网格内选得$\left({\lambda }_1^{\ast },{\lambda }_2^{\ast }\right)=\left(0.3,0.01\right)$ ，最小CV误差为966.29。

Figure 1. Heat map of five-fold cross-validation error for the SCAD-L2 penalized GEE over a two-dimensional parameter grid

图1. SCAD-L2化惩罚GEE在二维参数网格下的五折交叉验证误差热力图

表1和表2汇报了在相关正态连续型响应数据中，GEE、PGEE与SLGEE三种方法在样本量$n=200$ 与$n=600$ 下50次重复实验的模拟结果。结果显示，相较于不具备变量选择能力的传统GEE以及仅依赖单一惩罚机制的PGEE方法，SLGEE在所有实验设定下均表现出更优的统计性能。具体而言，SLGEE的均方误差在各方法中始终最低，略小于0.1，且在不同工作相关结构下结果保持一致，表明其估计过程具有较好的稳健性。在变量选择方面，SLGEE能够在有效控制过拟合的同时显著降低真实信号变量的遗漏率，使变量识别准确率稳定提升至97%~99%，Jaccard相似系数保持在0.94以上，显示出其在高维纵向数据情形下兼顾选择准确性与稳定性的优势。随着样本量的增加，PGEE和GEE的均方误差均呈下降趋势，而SLGEE即使在样本量较小的情况下也保持较低且稳定的均方误差，表现出更优的估计性能。

Table 1. Simulation results for GEE, PGEE, and SLGEE in Simulation Study 1 with n = 200

表1. 模拟实验1中n = 200时GEE、PGEE和SLGEE的模拟结果

方法	均方误差		过拟合数		欠拟合数		识别准确率		Jaccard
	均值	方差	均值	方差	均值	方差	均值	方差	均值
$n=200,\rho =0.5$
GEE.indep	0.303	0.073	152.88	1.35	0	0.00	23.56%	0.67%	0.9791
GEE.exch	0.22	0.049	151.98	1.67	0	0.00	24.01%	0.84%	0.9704
GEE.AR-1	0.249	0.053	151.8	1.94	0.02	0.14	24.09%	0.96%	0.9684
PGEE.indep	0.388	0.077	16.2	5.31	10.36	1.89	86.72%	2.98%	0.4075
PGEE.exch	0.344	0.064	7.72	3.92	10.22	1.73	91.03%	2.21%	0.4996
PGEE.AR-1	0.368	0.059	9.2	4.62	10.58	1.54	90.11%	2.60%	0.4729
SLGEE.indep	0.098	0.021	17.66	5.87	2.94	1.10	89.70%	2.96%	0.5672
SLGEE.exch	0.092	0.022	8.72	4.61	2.34	1.19	94.47%	2.22%	0.7008
SLGEE.AR-1	0.095	0.022	10.56	5.29	2.38	1.26	93.53%	2.63%	0.6656
$n=200,\rho =0.8$
GEE.indep	0.221	0.052	153.16	1.11	0.02	0.14	23.41%	0.56%	0.9816
GEE.exch	0.097	0.023	152.1	1.71	0	0.00	23.95%	0.85%	0.9717
GEE.AR-1	0.117	0.029	152	1.68	0	0.00	24.00%	0.84%	0.9705
PGEE.indep	0.375	0.057	24.24	5.37	10.18	1.30	82.79%	2.81%	0.3693
PGEE.exch	0.244	0.044	7.38	4.30	8.58	1.30	92.02%	2.29%	0.5454
PGEE.AR-1	0.267	0.045	9.9	5.44	8.9	1.34	90.60%	2.87%	0.5022
SLGEE.indep	0.1	0.021	26.1	6.34	3.34	0.72	85.28%	3.26%	0.5035
SLGEE.exch	0.086	0.021	7.52	4.78	1.28	1.05	95.60%	2.42%	0.7373
SLGEE.AR-1	0.087	0.020	10.6	5.65	1.68	1.13	93.86%	2.90%	0.6718

Table 2. Simulation results for GEE, PGEE, and SLGEE in Simulation Study 1 with n = 600

表2. 模拟实验1中n = 600时中GEE、PGEE和SLGEE的模拟结果

Methods	MSE		Overfitting		Underfitting		ER		Jaccard
	Mean	Var	Mean	Var	Mean	Var	Mean	Var	Mean
$n=600,\rho =0.5$
GEE.indep	0.093	0.021	150.84	2.23	0	0.00	24.58%	1.12%	0.9596
GEE.exch	0.063	0.015	148.58	2.42	0.02	0.14	25.70%	1.20%	0.9383
GEE.AR-1	0.071	0.018	149.96	2.83	0	0.00	25.02%	1.41%	0.9512
PGEE.indep	0.156	0.042	3.4	1.23	5.38	1.63	95.61%	1.09%	0.7043
PGEE.exch	0.139	0.038	0.48	0.58	5.36	1.59	97.08%	0.87%	0.7935
PGEE.AR-1	0.144	0.040	0.78	0.71	5.4	1.62	96.91%	0.82%	0.7818
SLGEE.indep	0.103	0.019	1.82	1.38	3.12	1.22	97.53%	0.97%	0.8924
SLGEE.exch	0.098	0.020	0.2	0.45	2.34	1.44	98.73%	0.71%	0.9476
SLGEE.AR-1	0.099	0.018	0.34	0.59	2.66	1.27	98.50%	0.65%	0.9463
$n=600,\rho =0.8$
GEE.indep	0.088	0.019	151.98	1.76	0	0.00	24.01%	0.88%	0.9704
GEE.exch	0.027	0.005	149.16	2.46	0	0.00	25.42%	1.23%	0.9439
GEE.AR-1	0.034	0.006	150.04	2.40	0	0.00	24.98%	1.20%	0.9519
PGEE.indep	0.164	0.042	9.16	2.22	5.76	1.42	92.54%	1.26%	0.5806
PGEE.exch	0.114	0.034	0.7	0.86	5.24	1.44	97.03%	0.81%	0.7953
PGEE.AR-1	0.12	0.033	1.1	1.02	5.3	1.40	96.80%	0.81%	0.7796
SLGEE.indep	0.097	0.024	8	2.77	2.76	1.19	94.62%	1.52%	0.7145
SLGEE.exch	0.079	0.021	0.12	0.33	1.42	1.28	99.23%	0.67%	0.9542
SLGEE.AR-1	0.08	0.022	0.32	0.68	1.16	1.20	99.26%	0.64%	0.9515

3.2. 实验2

相关二元响应数据的边际均值${\pi }_{ij}$ 满足

$\mathrm{Pr}\left(Y_{ij}=1|X_{ij}\right)=\frac{\exp X_{ij}^{\text{T}}\beta }{1+\exp X_{ij}^{\text{T}}\beta },i=1,\cdots ,n;j=1,\cdots ,m,$

令$n=200,m=10$ ，$X_{ij}={\left(x_{ij,1},\cdots ,x_{ij,100}\right)}^{\text{T}}$ 为包含100个协变量的向量，真实参数设为

$\beta ={\left(\underset{10}{\underbrace{0.7,\cdots ,0.7}},\underset{10}{\underbrace{-0.7,\cdots ,-0.7}},\underset{10}{\underbrace{-0.4,\cdots ,-0.4}},\underset{70}{\underbrace{0,\cdots ,0}}\right)}^{\text{T}},$

协变量的生成方式为

$\begin{array}{l}{X}_{jk}=Z_1+{\epsilon }_{jk}^x,Z_{1,j}~\text{Uniform}\left(0,0.5\right),j=1,\cdots ,10;k=1,\cdots ,10,\\ X_{jk}=Z_{2,j}+{\epsilon }_{jk}^x,Z_{2,j}~\text{Uniform}\left(0,0.5\right),j=1,\cdots ,10;k=11,\cdots ,20,\\ X_{jk}=Z_{3,j}+{\epsilon }_{jk}^x,Z_{3,j}~\text{Uniform}\left(0,0.5\right),j=1,\cdots ,10;k=21,\cdots ,30,\\ X_{jk}~\text{Uniform}\left(0,0.5\right),j=1,\cdots ,10;k=31,\cdots ,100,\end{array}$

其中${\epsilon }_{jk}^x{~}_{i.i.d.}N\left(0,0.001\right)$ 。

表3与表4给出了相关二元响应数据情形下的模拟结果，该情形相较于连续型响应具有更强的非线性特征和信息稀疏性。数值结果显示，SLGEE在该建模场景下仍能够保持良好的估计与选择性能，其回归系数均方误差在不同样本量设定下均处于较低水平，仅有0.03左右，表明该方法对二元响应模型具有较好的适应性。在变量选择方面，随着样本量的增加，SLGEE能够保持对过拟合与欠拟合的有效控制，当样本量达到600时，其过拟合变量数接近于0，欠拟合变量数稳定在1~2个，使变量识别准确率提升至95%~98%，Jaccard相似系数超过0.85，显示出较高的选择准确性与稳定性。并且在不同工作相关结构及相关性的情形下，SLGEE的估计精度和变量选择表现变化较小，体现了其在二元响应纵向数据分析中的稳健性。

Table 3. Simulation results for GEE, PGEE, and SLGEE in Simulation Study 2 with n = 200

表3. 模拟实验2中n = 200时GEE、PGEE和SLGEE的模拟结果

Methods	MSE		Overfitting		Underfitting		ER		Jaccard
	Mean	Var	Mean	Var	Mean	Var	Mean	Var	Mean
$n=200,\rho =0.5$
GEE.indep	0.290	0.019	69.88	0.14	0.13	0.14	30.00%	0.00%	0.9950
GEE.exch	0.281	0.019	69.88	0.14	0.00	0.00	30.13%	0.00%	0.9975
GEE.AR-1	0.287	0.019	69.88	0.14	0.00	0.00	30.13%	0.00%	0.9975
PGEE.indep	0.279	0.019	35.75	0.77	7.63	0.67	56.63%	1.97%	0.4269
PGEE.exch	0.269	0.019	34.25	0.72	7.50	0.60	58.25%	2.11%	0.4154
PGEE.AR-1	0.276	0.019	35.38	0.75	7.63	0.67	57.00%	2.05%	0.4235
SLGEE.indep	0.026	0.002	6.75	0.82	2.13	0.50	91.13%	0.84%	0.6364
SLGEE.exch	0.026	0.002	4.00	0.86	1.63	0.56	94.38%	1.25%	0.7270
SLGEE.AR-1	0.026	0.002	6.25	0.76	2.00	0.52	91.75%	0.76%	0.6543
$n=200,\rho =0.8$
GEE.indep	0.444	0.037	69.75	0.19	0.00	0.00	30.25%	0.19%	0.9950
GEE.exch	0.418	0.034	69.75	0.28	0.00	0.00	30.25%	0.28%	0.9950
GEE.AR-1	0.443	0.036	70.00	0.00	0.00	0.00	30.00%	0.00%	1.0000
PGEE.indep	0.424	0.031	39.88	1.88	8.00	0.68	52.13%	1.89%	0.4469
PGEE.exch	0.399	0.028	39.00	1.97	7.88	0.71	53.13%	1.95%	0.4381
PGEE.AR-1	0.419	0.029	39.75	2.12	8.25	0.68	52.00%	1.89%	0.4423
SLGEE.indep	0.029	0.002	9.50	0.60	2.13	0.26	88.38%	0.71%	0.5809
SLGEE.exch	0.034	0.001	2.00	0.43	2.25	0.28	95.75%	0.67%	0.7832
SLGEE.AR-1	0.029	0.002	7.13	0.84	2.00	0.30	90.88%	1.08%	0.6263

Table 4. Simulation results for GEE, PGEE, and SLGEE in Simulation Study 2 with n = 600

表4. 模拟实验2中n = 600时GEE、PGEE和SLGEE的模拟结果

Methods	MSE		Overfitting		Underfitting		ER		Jaccard
	Mean	Var	Mean	Var	Mean	Var	Mean	Var	Mean
$n=600,\rho =0.5$
GEE.indep	0.081	0.003	69.75	0.19	0.00	0.00	30.25%	0.19%	0.9950
GEE.exch	0.078	0.004	69.88	0.14	0.00	0.00	30.13%	0.14%	0.9975
GEE.AR-1	0.080	0.003	69.63	0.21	0.00	0.00	30.38%	0.21%	0.9925
PGEE.indep	0.071	0.003	13.88	1.49	7.50	0.64	78.63%	1.58%	0.3880
PGEE.exch	0.071	0.004	12.25	1.73	8.13	0.50	79.63%	1.76%	0.3691
PGEE.AR-1	0.071	0.003	13.63	1.57	7.63	0.64	78.75%	1.68%	0.3826
SLGEE.indep	0.022	0.001	0.38	0.21	1.50	0.21	98.13%	0.26%	0.8873
SLGEE.exch	0.025	0.001	0.13	0.14	1.63	0.30	98.25%	0.28%	0.8980
SLGEE.AR-1	0.022	0.001	0.25	0.18	1.75	0.28	98.00%	0.37%	0.8872
$n=600,\rho =0.8$
GEE.indep	0.128	0.006	69.88	0.14	0.25	0.18	29.88%	0.26%	0.9925
GEE.exch	0.114	0.006	70.00	0.00	0.00	0.00	30.00%	0.00%	1.0000
GEE.AR-1	0.124	0.006	69.75	0.19	0.13	0.14	30.13%	0.14%	0.9925
PGEE.indep	0.120	0.007	23.38	1.46	8.75	0.90	67.88%	1.29%	0.3496
PGEE.exch	0.107	0.006	20.50	1.19	8.88	0.78	70.63%	1.42%	0.3379
PGEE.AR-1	0.117	0.007	22.75	1.40	8.88	0.89	68.38%	1.50%	0.3445
SLGEE.indep	0.023	0.001	2.13	0.85	1.63	0.65	96.25%	1.18%	0.8007
SLGEE.exch	0.032	0.001	0.00	0.00	1.88	0.50	98.13%	0.50%	0.8952
SLGEE.AR-1	0.025	0.001	1.25	0.46	1.75	0.47	97.00%	0.74%	0.8313

4. 结论

本文提出了一种将SCAD-L2惩罚融入广义估计方程框架的高维变量选择方法，以在纵向数据的协变量存在强相关性时同时实现参数估计的稀疏性与稳定性。理论分析表明，在“大n，发散p”的高维设定下，该方法的参数估计值满足估计一致性与Oracle性质，从而确保其在高维情境下的理论可靠性。模拟实验进一步从多角度验证了该方法在参数估计精度与变量选择稳定性方面的显著优势，并显示其在不同相关结构和响应类型下均具有稳健的统计性能，体现了其在复杂纵向数据分析中的有效性与实用价值。

NOTES

^*通讯作者。