纵向数据下变系数线性模型的惩罚秩回归

Penalized Rank Regression for Varying-Coefficient Linear Models with Longitudinal Data

摘要: 针对纵向数据变系数线性模型中数据具有组内相关性和组间独立性，以及组内协方差矩阵估计效率低的特点，本文基于修正的Cholesky分解和B样条函数处理非参数部分，对组内协方差矩阵提出更有效的估计，改善了不平衡纵向数据的估计效率。同时，文章结合了两个非凸惩罚函数SCAD和MCP实现稳健的变量选择，模拟研究和实证分析表明新方法能够获得更有效的估计，并且相同参数准则下，对高相关性的样本数据，MCP惩罚函数的显著性更强。

Abstract: For longitudinal data in varying coefficient linear models characterized by within-group correlation and between-group independence, as well as the low estimation efficiency of within-group covariance matrices, this paper proposes a more efficient estimation method for the within-group covariance matrix based on modified Cholesky decomposition and B-spline functions to handle the nonparametric components, thereby improving the estimation efficiency for unbalanced longitudinal data. Additionally, the study incorporates two non-convex penalty functions, SCAD and MCP, to achieve robust variable selection. Simulation studies and empirical analyses demonstrate that the proposed method yields more efficient estimates. Under the same parameter criteria, the MCP penalty exhibits stronger significance for highly correlated sample data.

1. 引言

纵向数据是指通过对观测个体在不同时间节点上进行多次重复测量而获得的数据。因此，纵向数据分析可以对同一个体内部的变化规律以及协变量对响应变量的影响进行动态分析。在独立数据的回归分析中引入数据的组内相关性，提高统计推断的效率的同时，纵向数据分析也对统计方法提出了巨大的挑战[1]-[3]。常见的适用于纵向数据的模型，如边际回归模型、转移回归模型和混合效应模型。传统的参数回归模型只考虑参数部分，在处理实际问题，尤其是分析长期或动态环境下的数据时，模型解释能力不足。因此，考虑系数参数随时间变化是提升模型适应性的重要方向。变系数线性模型的首次提出是由Shumway [4]，Cleveland [5]和Hastie以及Tibshirani [6]给出的。其他基于该模型的讨论和扩展可从Fan和Zhang [7]以及Chiang等人[8]的论文中得到。

传统的变系数线性模型通常利用最小二乘方法估计系数参数，然而其对异常值敏感且依赖误差分布假设。但秩回归作为一种能有效降低异常值影响的非参数回归方法，它的引入恰好弥补了这一缺陷。该方法也在纵向数据下得到了广泛研究，如将基于秩回归的方法运用到纵向数据中[9]-[11]，并扩展到变量选择研究中。Fu等[12]基于独立性假设提出了基于秩的变量选择方法；Fu等[13]基于加权秩提出了带有SCAD (平滑裁剪绝对偏差)惩罚的Wilcoxon离散函数。但是在估计纵向数据的协方差矩阵时，大部分文献选择提前指定具体结构，如Jung等人[9]在独立工作结构下构造了广义的Wilcoxon-Mann-Whitney秩统计量，但得到的估计不是渐进有效的。Fu等人[14]在可交换结构下考虑了基于秩的估计函数，但是当组内相关结构不是可交换时估计效率会降低。因此，建立协方差模型，在更一般的相关结构下研究变系数模型是值得关注的话题。吕晶等[15]利用修正的Cholesky分解处理了纵向秩回归的组内相关性，提出了更加有效的无偏性秩回归得分函数。

以上提到的文献模型虽然在变系数估计方面取得了显著进展，但大多未系统解决变量选择问题。往往在实际建模中，研究者倾向于初始纳入过多变量以避免遗漏重要预测因子，但这一做法不可避免地导致模型过拟合、预测精度下降及计算效率损失。对此，Wang等人[16]结合了基差函数扩展和SCAD惩罚函数提出了正则化变量选择方法，成功识别并选出具有时变效应的重要变量。Zhang [17]提出了MCP (极小极大凹惩罚)惩罚函数，在保留SCAD惩罚优点的同时对回归系数进行有差别惩罚，从而得到更加精确的估计，同时证明了基于MCP惩罚的回归模型具有Oracle性质。但是相较SCAD惩罚函数，基于MCP惩罚的函数形式更为简洁、参数估计偏差更小、模型结果也更加精简。

综合上述讨论，本文结合B样条函数近似和修正的Cholesky分解提出一种纵向数据下变系数线性模型的秩回归方法。利用B样条光滑函数近似变系数参数，并在修正的Cholesky分解下估计更具一般性结v构的协方差工作矩阵。然而，基于秩回归的估计函数非凸，不连续且不可微。本文将利用Brown等人[18]提出的诱导光滑的方法获得更精确的估计。在光滑估计函数的基础下，文章分别添加了非凸惩罚函数SCAD和MCP进行对比分析，以针对不同相关性数据，实现更有效的变量选择。

2. 基于修正的Cholesky分解的秩回归

2.1. B样条近似

假定纵向数据$\left\{\left(t_{ik},X_{ik},Y_{ik}\right),1\le i\le n,1\le k\le n_i\right\}$ 满足如下形式的变系数模型：

$Y_{ik}=X_{ik}^{\text{T}}\beta \left(t_{ik}\right)+{\epsilon }_{ik},$ (1)

其中，$Y_{ik}$ 为响应变量，$X_{ik}={\left(X_{ik}^1,X_{ik}^2,\cdots ,X_{ik}^p\right)}^{\text{T}}$ 是一个$p$ 维协变量，$\beta \left(t_{ik}\right)={\left({\beta }_1\left(t_{ik}\right),{\beta }_2\left(t_{ik}\right),\cdots ,{\beta }_p\left(t_{ik}\right)\right)}^{\text{T}}$ ，${\beta }_r\left(t\right)$ $\left(r=1,2,\cdots ,p\right)$ 是未知的非参数函数，并且${\epsilon }_{ik}$ 是第$i$ 个个体在$t_{ik}$ 时刻测得的零均值随机过程。

对于模型的非参数函数，利用B样条函数近似未知变系数${\beta }_r\left(t\right)$ ，可得到${\beta }_r\left(t\right)\approx B_r\left(t\right){\gamma }_r={\displaystyle \sum _{h=1}^H{\gamma }_{rh}{B}_{rh}\left(t\right)}$ 。其中，$B_{rh}\left(t\right)\ge 0$ ，${\displaystyle \sum _{h=1}^H{B}_{rh}\left(t\right)}=1$ ，且控制点总数$H=h+m+1$ ，$m$ 为样条曲线的阶数，$h$ 为内节点个数。经过上述样条近似，可得到如下模型：

$Y_{ik}\approx {\displaystyle \sum _{r=1}^p{\displaystyle \sum _{h=1}^H{X}_{ik}^{\left(r\right)}}}{\gamma }_{rh}{B}_{rh}\left(t_{ik}\right)+{\epsilon }_{ik},$ (2)

现定义${\pi }^{\left(r\right)}(·)={\left(B_{r1}(·),\cdots ,B_{rH}(·)\right)}^{\text{T}}$ 为标准化B样条基函数，定义$\prod \left(T,X\right)={\left(X^{\left(1\right)}{\pi }^{\left(1\right)}{\left(T\right)}^{\text{T}},\cdots ,X^{\left(p\right)}{\pi }^{\left(p\right)}{\left(T\right)}^{\text{T}}\right)}^{\text{T}}$ ，${\prod }_{ik}=\prod \left(T_{ik},X_{ik}\right)$ ，${\pi }_{ik}^{\left(r\right)}={\pi }^{\left(r\right)}\left(t_{ik}\right)$ 。因此，(2)可以写为$Y_{ik}\approx {\prod }_{ik}^T\gamma +{\epsilon }_{ik}$ ，其中$\gamma ={\left({\gamma }_1^{\text{T}},\cdots ,{\gamma }_p^{\text{T}}\right)}^{\text{T}}$ ，${\gamma }_r={\left({\gamma }_{r1},\cdots ,{\gamma }_{rH}\right)}^{\text{T}}$ 。

基于上面的近似，我们可以得到在$t_{ik}$ 处的残差${\widehat{\epsilon }}_{ik}\left(\widehat{\gamma }\right)=Y_{ik}-{\prod }_{ik}^{\text{T}}\widehat{\gamma }$。对于模型(2)，定义$S_{ik}$ 为${\epsilon }_{ik}\left(\gamma \right)$ 的平均秩，$S_{ik}=M^{-1}{\displaystyle \sum _{j=1}^n{\displaystyle \sum _{l=1}^{n_j}\left[I\left\{{\epsilon }_{jl}\left(\gamma \right)<{\epsilon }_{ik}\left(\gamma \right)\right\}-\frac{1}{2}\right]}}$ ，其中$M={\displaystyle \sum _{i=1}^n{n}_i}$ ，记$S_i\left(\gamma \right)={\left(S_{i1},\cdots ,S_{i{n}_i}\right)}^{\text{T}}$ ， $S\left(\gamma \right)={\left\{S_1{\left(\gamma \right)}^{\text{T}},\cdots ,S_n{\left(\gamma \right)}^{\text{T}}\right\}}^{\text{T}}$ 。

对$S\left(\gamma \right)$ 有秩估计函数：

$U_{OPT}\left(\gamma \right)={\overline{D}}^{\text{T}}{V}^{-1}S\left(\gamma \right),$ (3)

其中$V=Cov\left(S\left(\gamma \right)\right)$ ，$\overline{D}=\left\{{\overline{D}}_{ik},i=1,\cdots ,n;k=1,\cdots ,n_i\right\}$ 为$S\left(\gamma \right)$ 的期望关于$\gamma$ 的导数，且${\overline{D}}_{ik}\left({\gamma }_0\right)=M^{-1}{\displaystyle \sum _{j=1}^n{\displaystyle \sum _{l=1}^{n_j}\left({\Pi }_{jl}^{\text{T}}-{\Pi }_{ik}^{\text{T}}\right)f_{ikjl}\left(0\right)}}$，$f_{ikjl}(·)$ 为${\epsilon }_{ik}-{\epsilon }_{jl}$ 的密度函数，${\overline{D}}_{ik}$ 涉及未知的密度函数$f_{ikjl}(·)$ 。假设$f_{ikjl}\left(0\right)$ 为常数，并在模型(3)中用$D_{ik}^{\text{T}}=\overline{\Pi }-{\Pi }_{ik}$代替${\overline{D}}_{ik}^{\text{T}}$。

假定使用一个块对角矩阵$diag\left(V_1,\cdots ,V_n\right)$ 作为$V$ 的工作协方差矩阵，可得到$\gamma$ 的估计函数为$U_B\left(\gamma \right)={\displaystyle \sum _{i=1}^n{D}_i^{\text{T}}{V}_i^{-1}{S}_i\left(\gamma \right)}$ ，其中，$V_i=\mathrm{var}\left\{S_i\left(\gamma \right)\right\}$ 。

2.2. 一般工作结构的秩估计函数

在以往的秩估计方法中，Jung和Ying [9]假定误差具有独立工作结构用于估计$\beta$ ，但估计结果${\widehat{\beta }}_I$ 在其他工作结构情况下会严重损失估计效率。基于此，Fu和Wang [14]在可交换结构下提出了最佳秩估计函数，并利用矩方法估计协方差矩阵的逆中的参数${\sigma }_i^2$ 和${\sigma }_{ii}$ 。在该模型下，只有${\epsilon }_i$ 具有可交换结构时，$S_i$ 和${\epsilon }_i$ 才具有相同的相关结构。因此，当${\epsilon }_i$ 不具有可交换相关结构时，${\widehat{\beta }}_F$ 会损失估计效率。

已知上述两个模型只有在给定结构下的估计才是有效的，当误差具有其它结构时，估计效率大大降低。此时，修正的Cholesky分解可以保证协方差矩阵正定，同时提高数值计算的稳定性。故本文将基于修正的Cholesky分解建模，基于Yao和Li [19]，$V_i=Cov\left\{S_i\left({\gamma }_0\right)\right\}$ 可进行如下分解，

$Cov\left\{{\Phi }_i{S}_i\right\}={\Phi }_i^{\text{T}}{V}_i{\Phi }_i={\Lambda }_i$ (4)

则$V_i={\left({\Phi }_i^{-1}\right)}^{\text{T}}{\Lambda }_i{\Phi }_i^{-1}$ 。其中，${\Phi }_i$ 是主对角元素均为1的下三角矩阵，${\Lambda }_i=diag\left(d_{i1}^2,\cdots ,d_{i{n}_i}^2\right)$ 为$n_i\times n_i$ 的对角矩阵。记$e_i={\left(e_{i1},\cdots ,e_{i{n}_i}\right)}^{\text{T}}={\Phi }_i{S}_i$ ，${\Phi }_i$ 的第$\left(k,l\right)$ 个元素为$-{\varphi }_{k,l}^{\left(i\right)}\left(k=1,\cdots ,n_i;l=1,\cdots ,k-1\right)$ ，则有如下关系式：

$\{\begin{array}{l}{S}_{i1}=e_{i1}\\ S_{ik}={\varphi }_{k,1}^{\left(i\right)}{S}_{i1}+\cdots +{\varphi }_{k,k-1}^{\left(i\right)}{S}_{i,k-1}+e_{ik},i=1,2,\cdots ,n;k=2,\cdots ,n_i\end{array}$

已知${\Lambda }_i$ 为对角矩阵，故$e_{ik}$ 互不相关，且$\mathrm{var}\left(e_{ik}\right)=d_{ik}^2$ 。根据Liu和Li [20]的思想，在此构造线性模型和非参数模型，引入$s\times 1$ 维系数向量$\theta$ 和一维未知光滑函数估计自回归系数${\varphi }_{k,l}^{\left(i\right)}=z_{k,l}^{\left(i\right)\text{T}}\theta$ 和更新方差$d_{ik}^2=d_{ik}^2\left(t_{ik}\right)=g\left(t_{ik}\right)$ 。此处取$z_{k,l}^{\left(i\right)}={\left(1,t_{ik}-t_{il},\cdots ,{\left(t_{ik}-t_{il}\right)}^{s-1}\right)}^{\text{T}}$ 。

利用最小二乘法对参数向量$\theta$ 进行估计，求解下面的最小二乘问题即可得到估计值$\widehat{\theta }$ ：

$\begin{array}{c}\widehat{\theta }=\arg \min {\displaystyle \sum _{i=1}^n{\displaystyle \sum _{k=2}^{n_i}{\left({\widehat{S}}_{ik}-{\widehat{S}}_{i1}{z}_{k,1}^{\left(i\right)\text{T}}\theta -\cdots -{\widehat{S}}_{i,k-1}{z}_{k,k-1}^{\left(i\right)\text{T}}\theta \right)}^2}}\\ =\arg \min {\displaystyle \sum _{i=1}^n{\displaystyle \sum _{k=2}^{n_i}{\left({\widehat{S}}_{ik}-{\widehat{\upsilon }}_{ik}\theta \right)}^2}}\\ ={\left({\displaystyle \sum _{i=1}^n{\displaystyle \sum _{k=2}^{n_i}\left({\widehat{\upsilon }}_{ik}{\widehat{\upsilon }}_{ik}^{\text{T}}\right)}}\right)}^{-1}\left({\displaystyle \sum _{i=1}^n{\displaystyle \sum _{k=2}^{n_i}\left({\widehat{\upsilon }}_{ik}{\widehat{S}}_{ik}\right)}}\right)\end{array}$

其中，${\widehat{\upsilon }}_{ik}={\left({\displaystyle \sum _{l=1}^{k-1}{\widehat{S}}_{il}{z}_{k,l,1}^{\left(i\right)\text{T}}},\cdots ,{\displaystyle \sum _{l=1}^{k-1}{\widehat{S}}_{il}{z}_{k,l,s}^{\left(i\right)\text{T}}}\right)}^{\text{T}}$。

根据Fan and Yao [21]的思想，对光滑函数$g\left(t_{ik}\right)$ 的估计可通过求解如下最小二乘问题得到：

$\left({\widehat{\eta }}_1,{\widehat{\eta }}_2\right)=\arg \min {\displaystyle \sum _{i=1}^n{\displaystyle \sum _{k=2}^{n_i}{\left[{\widehat{e}}_{ik}^2-{\eta }_1-{\eta }_2\left(t_{ik}-t\right)\right]}^2{K}_h\left(t_{ik}-t\right)}}$

其中，$K_h\left(t_{ik}-t\right)=K\left(\frac{t_{ik}-t}{h}\right)/h$ ，h为带宽，$K(·)$ 为核函数，通过上述问题的求解可知$d_{ik}^2\left(t_{ik}\right)={\widehat{\eta }}_1$ 。

结合上述讨论，可得到估计函数$U_C\left(\gamma \right)={\displaystyle \sum _{i=1}^n{D}_i^{\text{T}}{\widehat{V}}_i^{-1}{S}_i\left(\gamma \right)}$。估计函数的解记为${\gamma }_C$ ，求出结果需要对估计函数求导，但$S_{ik}=M^{-1}{\displaystyle \sum _{j=1}^n{\displaystyle \sum _{l=1}^{n_j}\left[I\left\{{\epsilon }_{jl}\left(\gamma \right)<{\epsilon }_{ik}\left(\gamma \right)\right\}-\frac{1}{2}\right]}}$ 不可微，则采用Brown B和Wang [18]提出的诱导光滑方法得到$S_{ik}$ 的近似式，即

${\overline{S}}_{ik}\left(\gamma \right)=M^{-1}{\displaystyle \sum _{j=1}^n{\displaystyle \sum _{l=1}^{n_j}\left[\Phi \left({\epsilon }_{ikjl}/r_{ikjl}\right)-0.5\right]}},\left(i,k\right)\ne \left(j,l\right)$ (5)

得到诱导光滑秩估计函数为：

${\tilde{U}}_C\left(\gamma \right)={\displaystyle \sum _{i=1}^n{D}_i^{\text{T}}{\widehat{V}}_i^{-1}{\overline{S}}_i\left(\gamma \right)}$ (6)

其中，$\Phi (·)$ 表示标准正态的卷积分布函数，${\epsilon }_{ikjl}={\epsilon }_{ik}-{\epsilon }_{jl}$ ，$r_{ikjl}^2=d_{ikjl}^{\text{T}}\Omega d_{ikjl}$ ，$d_{ikjl}={\Pi }_{ik}-{\Pi }_{jl}$ ，$\Omega$ 是${\tilde{\gamma }}_C$ 的$p\times p$ 维协方差矩阵，${\tilde{\gamma }}_C$ 为诱导光滑秩估计函数的解，${\overline{S}}_i\left(\gamma \right)={\left({\overline{S}}_{i1}\left(\gamma \right),\cdots ,{\overline{S}}_{i{n}_i}\left(\gamma \right)\right)}^{\text{T}}$。此时${\tilde{U}}_C\left(\gamma \right)$ 函数的特性满足了我们的需求。

2.3. 迭代算法

已知上述求导结果，算法具体迭代步骤如下：

步骤1 给定初始值${\tilde{\gamma }}^{\left(0\right)}={\left(0,0,\cdots ,0\right)}^{\text{T}}$和${\Omega }^{\left(0\right)}=n_{-1}{I}_p$ ；

步骤2 利用第$m$ 次迭代得到的${\tilde{\gamma }}^{\left(m\right)}$ 和${\tilde{\Omega }}^{\left(m\right)}$ ，在下面公式的基础上更新${\tilde{\gamma }}^{\left(m+1\right)}$ 和${\tilde{\Omega }}^{\left(m+1\right)}$ ；

${\tilde{\gamma }}^{\left(m+1\right)}={\tilde{\gamma }}^{\left(m\right)}-{\left({\frac{\partial {\tilde{U}}_C\left(\gamma \right)}{\partial \gamma }|}_{{\tilde{\gamma }}^{\left(m\right)},{\tilde{\Omega }}^{\left(m\right)}}\right)}^{-1}{\tilde{U}}_C\left(\gamma \right)$ (7)

${\tilde{\Omega }}^{\left(m+1\right)}={\left({\frac{\partial {\tilde{U}}_C\left(\gamma \right)}{\partial \gamma }|}_{{\tilde{\gamma }}^{\left(m\right)},{\tilde{\Omega }}^{\left(m\right)}}\right)}^{-1}Cov\left({\tilde{U}}_C\left({\tilde{\gamma }}^{\left(m\right)},{\tilde{\Omega }}^{\left(m\right)}\right)\right){\left({\left({\frac{\partial {\tilde{U}}_C\left(\gamma \right)}{\partial \gamma }|}_{{\tilde{\gamma }}^{\left(m\right)},{\tilde{\Omega }}^{\left(m\right)}}\right)}^{-1}\right)}^{\text{T}}$ (8)

其中，$Cov\left({\tilde{U}}_C\left({\tilde{\gamma }}^{\left(m\right)},{\tilde{\Omega }}^{\left(m\right)}\right)\right)={\displaystyle \sum _{i=1}^n{D}_i^{\text{T}}{\widehat{V}}_i^{-1}{\tilde{S}}_i\left(\gamma \right){\tilde{S}}_i^{\text{T}}\left(\gamma \right){\widehat{V}}_i^{-1}{D}_i}$。

步骤3 重复上述步骤2直至收敛。

3. 变量选择与惩罚秩估计函数

在实际情况中，真实的模型通常是未知的，并且我们不知道哪些系数为0，但提前假定系数为0往往会造成损失，因此为避免错误假设带来的损失，考虑到参数的稀疏性问题，我们在上述诱导光滑秩估计函数的基础下引入满足无偏性的SCAD惩罚函数和MCP惩罚函数改进模型，通过构建惩罚秩估计函数得到系数参数更一致的估计。

3.1. 非凸惩罚函数

SCAD惩罚函数$p_{\lambda }^{SCAD}\left(\left|{\gamma }_k\right|\right)$ 形式如下：

$p_{\lambda }^{SCAD}\left(\left|\gamma \right|\right)=\{\begin{array}{ll}\lambda \left|\gamma \right|,\hfill & \left|\gamma \right|\le \lambda \hfill \\ \left(a^2-1\right){\lambda }^2-{\left(\left|\gamma \right|-a\lambda \right)}^2/2\left(a-1\right),\hfill & \lambda <\left|\gamma \right|\le a\lambda \hfill \\ \left(a+1\right){\lambda }^2/2,\hfill & \left|\gamma \right|>a\lambda \hfill \end{array}$ (9)

其中SCAD惩罚函数的一阶导为

$p_{\lambda }^{SCAD}{}^{\prime }\left(\left|\gamma \right|\right)=\{\begin{array}{ll}\lambda sign\left(\left|\gamma \right|\right),\hfill & \left|\gamma \right|\le \lambda \hfill \\ \left(a\lambda sign\left(\left|\gamma \right|\right)-\left|\gamma \right|\right)/\left(a-1\right),\hfill & \lambda <\left|\gamma \right|\le a\lambda \hfill \\ 0,\hfill & \left|\gamma \right|>a\lambda \hfill \end{array}$

已知SCAD惩罚函数没有连续的二阶导数，但是对于该惩罚函数可以通过一个二次函数近似。现假定有初始值${\gamma }_0$ ，惩罚函数的近似函数如下所示：

${\displaystyle \sum _{k=1}^p{p}_{\lambda }^{SCAD}\left(\left|{\gamma }_k\right|\right)}\approx {\displaystyle \sum _{k=1}^p\frac{1}{2}}\frac{p_{\lambda }^{SCAD}{}^{\prime }\left(\left|{\tilde{\gamma }}_k\right|\right)}{\left|{\tilde{\gamma }}_k\right|}{\left|{\gamma }_k\right|}^2={\gamma }^{\text{T}}{M}_{SCAD}\gamma$ (10)

在近似式下，得到近似的秩估计函数：

${\tilde{U}}_{SCAD}\left(\gamma \right)={\displaystyle \sum _{i=1}^n{D}_i^{\text{T}}{\widehat{V}}_i^{-1}{\overline{S}}_i\left(\gamma \right)}+n{\gamma }^{\text{T}}{M}_{SCAD}$ (11)

其中，$\gamma ={\left({\gamma }_1^{\text{T}},\cdots ,{\gamma }_p^{\text{T}}\right)}^{\text{T}}$ ，$M_{SCAD}=diag\left(\frac{p_{\lambda }^{SCAD}{}^{\prime }\left(\left|{\tilde{\gamma }}_1\right|\right)}{2\left|{\tilde{\gamma }}_1\right|},\cdots ,\frac{p_{\lambda }^{SCAD}{}^{\prime }\left(\left|{\tilde{\gamma }}_p\right|\right)}{2\left|{\tilde{\gamma }}_p\right|}\right)$ 。

可得到近似的估计等式为：${\displaystyle \sum _{i=1}^n{\displaystyle \sum _{k=1}^{n_i}{D}_{ik}^{\text{T}}{\widehat{V}}_{ik}^{-1}{\Delta }_{ik}\left(\gamma \right)}}+n{\displaystyle \sum _{k=1}^p{p}_{\lambda }^{SCAD}{}^{\prime }\left(\left|{\gamma }_k^0\right|\right)\left|{\gamma }_k\right|}=0$，其中，${\gamma }_k^0$ 为一个初始估计值，即非惩罚加权Wilcoxon估计量。

计算得到算法迭代式如下：

${\tilde{\gamma }}^{\left(m+1\right)}={\tilde{\gamma }}^{\left(m\right)}-{\left(\frac{\partial {\tilde{U}}_{SCAD}\left(\gamma \right)}{\partial \gamma }\right)}^{-1}{\tilde{U}}_{SCAD}\left(\gamma \right)$ (12)

${\tilde{\Omega }}^{\left(m+1\right)}={\left({\frac{\partial {\tilde{U}}_{SCAD}\left(\gamma \right)}{\partial \gamma }|}_{{\tilde{\gamma }}^{\left(m\right)},{\tilde{\Omega }}^{\left(m\right)}}\right)}^{-1}Cov\left({\tilde{U}}_{SCAD}\left({\tilde{\gamma }}^{\left(m\right)},{\tilde{\Omega }}^{\left(m\right)}\right)\right){\left({\left({\frac{\partial {\tilde{U}}_{SCAD}\left(\gamma \right)}{\partial \gamma }|}_{{\tilde{\gamma }}^{\left(m\right)},{\tilde{\Omega }}^{\left(m\right)}}\right)}^{-1}\right)}^{\text{T}}$ (13)

其中，$\frac{\partial {\tilde{U}}_{SCAD}\left(\gamma \right)}{\partial \gamma }={\displaystyle \sum _{i=1}^n{\displaystyle \sum _{k=1}^{n_i}{D}_{ik}^{\text{T}}{\widehat{V}}_{ik}^{-1}{\Delta }_{ik}\left(\gamma \right)}}+n{M}_{SCAD}$，

$Cov\left({\tilde{U}}_{SCAD}\left({\tilde{\gamma }}^{\left(m\right)},{\tilde{\Omega }}^{\left(m\right)}\right)\right)={\displaystyle \sum _{i=1}^n{D}_i^{\text{T}}{\widehat{V}}_i^{-1}{\tilde{S}}_i\left({\tilde{\gamma }}^{\left(m\right)}\right){\tilde{S}}_i^{\text{T}}\left({\tilde{\gamma }}^{\left(m\right)}\right){\widehat{V}}_i^{-1}{D}_i}+n{M}_{SCAD}^{\text{T}}{\tilde{\gamma }}^{\left(m\right)}{\tilde{\gamma }}^{\left(m\right)}{}^{\text{T}}{M}_{SCAD}$。

MCP惩罚函数也同样满足Oracle性质，且在处理特征之间有很高相关性数据时，表现要比SCAD惩罚函数更好。MCP惩罚函数的惩罚项$p_{\lambda }^{MCP}\left(\left|{\gamma }_k\right|\right)$ 为：

$p_{\lambda }^{MCP}\left(\left|\gamma \right|\right)=\{\begin{array}{ll}\lambda \left|\gamma \right|-\frac{{\left|\gamma \right|}^2}{2a},\hfill & \left|\gamma \right|\le a\lambda \hfill \\ \frac{1}{2}a{\lambda }^2,\hfill & \left|\gamma \right|>a\lambda \hfill \end{array}$ (14)

MCP惩罚函数的一阶导为

$p_{\lambda }^{MCP}{}^{\prime }\left(\left|\gamma \right|\right)=\{\begin{array}{ll}\lambda sign\left(\left|\gamma \right|\right)-\frac{\left|\gamma \right|}{a},\hfill & \left|\gamma \right|\le a\lambda \hfill \\ 0,\hfill & \left|\gamma \right|>a\lambda \hfill \end{array}$

同样利用二次函数对MCP惩罚近似，得到MCP惩罚下的秩估计函数为

${\tilde{U}}_{MCP}\left(\gamma \right)={\displaystyle \sum _{i=1}^n{D}_i^{\text{T}}{\widehat{V}}_i^{-1}{\overline{S}}_i\left(\gamma \right)}+n{\gamma }^{\text{T}}{M}_{MCP},$ (15)

其中，$M_{MCP}=diag\left(\frac{p_{\lambda }^{MCP}{}^{\prime }\left(\left|{\tilde{\gamma }}_1\right|\right)}{2\left|{\tilde{\gamma }}_1\right|},\cdots ,\frac{p_{\lambda }^{MCP}{}^{\prime }\left(\left|{\tilde{\gamma }}_p\right|\right)}{2\left|{\tilde{\gamma }}_p\right|}\right)$ 。则有估计等式${\displaystyle \sum _{i=1}^n{\displaystyle \sum _{k=1}^{n_i}{D}_{ik}^{\text{T}}{\widehat{V}}_{ik}^{-1}{\Delta }_{ik}\left(\gamma \right)}}+n{\displaystyle \sum _{k=1}^p{p}_{\lambda }^{MCP}{}^{\prime }\left(\left|{\gamma }_k^0\right|\right)\left|{\gamma }_k\right|}=0$。

同样可得算法迭代式：

${\tilde{\gamma }}^{\left(m+1\right)}={\tilde{\gamma }}^{\left(m\right)}-{\left(\frac{\partial {\tilde{U}}_{MCP}\left(\gamma \right)}{\partial \gamma }\right)}^{-1}{\tilde{U}}_{MCP}\left(\gamma \right)$ (16)

${\tilde{\Omega }}^{\left(m+1\right)}={\left({\frac{\partial {\tilde{U}}_{MCP}\left(\gamma \right)}{\partial \gamma }|}_{{\tilde{\gamma }}^{\left(m\right)},{\tilde{\Omega }}^{\left(m\right)}}\right)}^{-1}Cov\left({\tilde{U}}_{MCP}\left({\tilde{\gamma }}^{\left(m\right)},{\tilde{\Omega }}^{\left(m\right)}\right)\right){\left({\left({\frac{\partial {\tilde{U}}_{MCP}\left(\gamma \right)}{\partial \gamma }|}_{{\tilde{\gamma }}^{\left(m\right)},{\tilde{\Omega }}^{\left(m\right)}}\right)}^{-1}\right)}^{\text{T}}$ (17)

其中，$\frac{\partial {\tilde{U}}_{MCP}\left(\gamma \right)}{\partial \gamma }={\displaystyle \sum _{i=1}^n{\displaystyle \sum _{k=1}^{n_i}{D}_{ik}^{\text{T}}{\widehat{V}}_{ik}^{-1}{\Delta }_{ik}\left(\gamma \right)}}+n{M}_{MCP}$，

$Cov\left({\tilde{U}}_{MCP}\left({\tilde{\gamma }}^{\left(m\right)},{\tilde{\Omega }}^{\left(m\right)}\right)\right)={\displaystyle \sum _{i=1}^n{D}_i^T{\widehat{V}}_i^{-1}{\tilde{S}}_i\left({\tilde{\gamma }}^{\left(m\right)}\right){\tilde{S}}_i^{\text{T}}\left({\tilde{\gamma }}^{\left(m\right)}\right){\widehat{V}}_i^{-1}{D}_i}+n{M}_{MCP}^{\text{T}}{\tilde{\gamma }}^{\left(m\right)}{\tilde{\gamma }}^{\left(m\right)}{}^{\text{T}}{M}_{MCP}$。

3.2. 最佳调节参数选择

在(9)和(14)式中，还有参数a需要估计，a的不同取值会直接影响估计性能。对于(9)式中的参数a，Fan和Li [22]通过蒙特卡洛模拟得出a的最优值约等于3.7，(14)式的a在实际使用中通常默认为3。对于参数$\lambda$ ，本文使用BIC (贝叶斯信息准则)进行参数选择。具体计算式如下：

${\text{BIC}}_{\lambda }=\log \left({\displaystyle \sum _{i=1}^n{D}_i^{\text{T}}{\widehat{V}}_i^{-1}{\overline{S}}_i\left(\gamma \right)}\right)+\frac{df\left(\lambda \right)\log \left(n\right)}{n}$ (18)

4. 模拟研究

数据来源于如下变系数模型

$Y_{ik}={\displaystyle \sum _{k=1}^{n_i}{X}_{ik}^{\text{T}}\beta \left(t_{ik}\right)}+{\epsilon }_{ik},i=1,\cdots ,n;j=1,\cdots ,n_i$ (19)

其中${\beta }_1\left(t\right)=\frac{3}{2}-\cos \left(\left(t-2\right)\pi \right)$ ，${\beta }_2\left(t\right)={\left(\frac{1}{4}-2t\right)}^2$ ，${\beta }_3\left(t\right)=-\frac{2}{3}+\frac{5.5-4t}{300}$ ，${\beta }_{p-3}\left(t\right)=0$ 。$t_{ik}~U\left(-1,1\right)$ ，$X_{ik}^1$ 服从$\left[t_{ik},3+t_{ik}\right]$ 上的均匀分布，在$X_{ik}^1$ 的条件下，$X_{ik}^2$ 产生于正态分布$N\left(0,{\Sigma }_2\right)$ 生成的随机数加1，其中${\Sigma }_2=\frac{X_{ik}^1+1}{X_{ik}^1+2}$ ；$X_{ik}^3$ 服从参数为0.8的伯努利分布，其余协变量$\left(X_{ik}^4,\cdots ,X_{ik}^{p-3}\right)$ 服从均值为0，协方差矩阵为$\Upsilon$ 的多元正态分布。

为对比不同相关性数据下惩罚函数变量选择的结果，考虑两种生成协变量的协方差矩阵，具体如下：

CASE Ⅰ协方差矩阵$\Upsilon =Cov\left(X_{ik}^r,X_{jl}^r\right)=\{\begin{array}{ll}15{\text{e}}^{-\left|t_{ik}-t_{il}\right|},\hfill & i=j\hfill \\ 0,\hfill & i\ne j\hfill \end{array}$ ；

CASE Ⅱ协方差矩阵$\Upsilon =Cov\left(X_{ik}^r,X_{jl}^r\right)=\{\begin{array}{ll}30{\text{e}}^{-0.3\left|t_{ik}-t_{il}\right|},\hfill & i=j\hfill \\ 0,\hfill & i\ne j\hfill \end{array}$ 。

随机误差${\epsilon }_i$ 服从多元正态分布$N\left(0,{\Xi }_i\right)$ 。这里协方差矩阵${\Xi }_i$ 采用定义如下的矩阵${\Xi }_i={\Delta }_i^{-1}{B}_i{\left({\Delta }_i^{\text{T}}\right)}^{-1}$ ，其中$B_i$ 是$m_i\times m_i$ 的对角矩阵且第$k$ 个对角元素为$\frac{\sin \left(\pi {\varsigma }_{ik}\right)}{3}+0.5$ ，${\varsigma }_{ik}~U\left(0,2\right)$ ，以及${\Delta }_i$ 是对角元素全为1的下三角矩阵且第$\left(k,l\right)$ 个元素为$-{\delta }_{k,l}^i\left(k\ne l\right)$ ，${\delta }_{k,l}^i=W_{k,l}^{\left(i\right)\text{T}}{\theta }^{*}$ ，其中$W_{k,l}^{\left(i\right)}={\left(1,{\omega }_{ik}-{\omega }_{il}\right)}^{\text{T}}$ ，${\theta }^{*}={\left(0.2,0.5\right)}^{\text{T}}$ 和${\omega }_{ik}~U\left(0,1\right)$ 。

对于重复观测次数$n_i$ ，此处我们考虑不平衡的情形，令$n_i$ 取2到10之间的随机数。取$n=30,50,100$ 为小中大样本，模拟次数为100次。

为了体现模型评估的准确性，我们采用均方根平均误差(RASE)来对模型结果进行评价，RASE的具体计算表达式如下：

$\text{RASE}={\left\{n^{-1}{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^{n_i}{\left({\widehat{\beta }}_j\left(U_i\right)-{\beta }_j\left(U_i\right)\right)}^2}}\right\}}^{\frac{1}{2}}$ (20)

将本文方法下的估计结果记为${\beta }_{PB}^{SCAD}$ 、${\beta }_{PB}^{MCP}$ ，与三种结构下的PGEE方法进行对比分析，估计结果记为${\beta }_{PGEE}^{indep}$ 、${\beta }_{PGEE}^{AR}$ 、${\beta }_{PGEE}^{CS}$ 。具体见表1。

Table 1. RASE results for CASE Ⅰ and CASE Ⅱ (${\beta }_{PB}^{SCAD}$ 、${\beta }_{PB}^{MCP}$ (×10⁻¹))

表1. CASE Ⅰ和CASE Ⅱ的RASE结果(其中${\beta }_{PB}^{SCAD}$ 、${\beta }_{PB}^{MCP}$ (×10⁻¹))

$n$	$p_n$		CASE Ⅰ	CASE Ⅱ
			RASE
30	12	${\beta }_{PGEE}^{indep}$	0.7663	0.2021
		${\beta }_{PGEE}^{AR}$	0.7637	0.1908
		${\beta }_{PGEE}^{CS}$	0.7025	0.2174
		${\beta }_{PB}^{SCAD}$	0.0946	0.0870
		${\beta }_{PB}^{MCP}$	0.0944	0.0870
	15	${\beta }_{PGEE}^{indep}$	0.6214	0.1936
		${\beta }_{PGEE}^{AR}$	0.7213	0.1523
		${\beta }_{PGEE}^{CS}$	0.6547	0.1237
		${\beta }_{PB}^{SCAD}$	0.0088	0.0073
		${\beta }_{PB}^{MCP}$	0.0077	0.0072
50	12	${\beta }_{PGEE}^{indep}$	0.6688	0.1939
		${\beta }_{PGEE}^{AR}$	0.7633	0.1417
		${\beta }_{PGEE}^{CS}$	0.6981	0.0871
		${\beta }_{PB}^{SCAD}$	0.0855	0.0850
		${\beta }_{PB}^{MCP}$	0.0795	0.0747
	15	${\beta }_{PGEE}^{indep}$	0.5652	0.1543
		${\beta }_{PGEE}^{AR}$	0.5012	0.1540
		${\beta }_{PGEE}^{CS}$	0.4980	0.1246
		${\beta }_{PB}^{SCAD}$	0.0654	0.0531
		${\beta }_{PB}^{MCP}$	0.0554	0.0526
100	12	${\beta }_{PGEE}^{indep}$	0.5165	0.2016
		${\beta }_{PGEE}^{AR}$	0.5275	0.2058
		${\beta }_{PGEE}^{CS}$	0.4932	0.2091
		${\beta }_{PB}^{SCAD}$	0.0393	0.0300
		${\beta }_{PB}^{MCP}$	0.0393	0.0290
	15	${\beta }_{PGEE}^{indep}$	0.3960	0.1180
		${\beta }_{PGEE}^{AR}$	0.3910	0.1324
		${\beta }_{PGEE}^{CS}$	0.2862	0.1461
		${\beta }_{PB}^{SCAD}$	0.0386	0.0319
		${\beta }_{PB}^{MCP}$	0.0376	0.0291

根据上面的表格，可以得到如下结论：

1) 不同维度下处理效果

为了体现本文提出的方法对不同维度样本数据的处理效果差异，逐步增加模型维数，分别取$p_n=12$ 和$p_n=15$ 。对比上述结果，在小中大三种样本类型下，模型维数越大，估计的均方根平均误差越小，说明模型的估计效率越高。

2) 不同协方差工作结构对估计结果的影响

基于文章研究，在相同惩罚函数和样本数据结构下，PGEE方法对应的模型平均误差略高于本文提出方法下的模型误差，说明利用修正的Cholesky分解处理变系数纵向秩回归模型的组内协方差矩阵是有效的，明显降低了估计效率。同时，在模拟研究中，会出现两种方法的模型误差非常接近的情况，这说明此时对协方差矩阵结构的两种估计都是正确的。总之，与现有方法相比，本文提出的方法具有更强的优势，主要原因在于它允许协方差结构具有更一般的形式，不需提前指定矩阵结构，从而可以达到更强的稳健性和有效性。

3) 不同协方差相关性对惩罚函数变量选择的影响

对比相同维度下CASE Ⅰ和CASE Ⅱ的模型结果，发现当样本协方差矩阵相关性变强时，CASE Ⅱ的模型平均误差更小。这是因为当变量高度相关时，MCP惩罚的硬阈值特性和无偏性使其能更有效地筛选关键变量并减少估计偏差，而SCAD的渐进惩罚可能导致冗余变量残留。结合两个惩罚的表达式分析，$\left|\gamma \right|>a\lambda$ 时，MCP惩罚完全为0，而SCAD仍然存在轻微惩罚，为非零常数。

5. 模拟研究

我们将本文提出的方法应用到纵向孕酮数据中，该数据的具体表述见文章[23]。这项研究收集了34名女性的尿样，每个妇女的重复观测次数从11到28变化且这组数据共有492个观测值，每位女性的观测次数覆盖了完整的月经周期。但由于不同女性的周期长度不同，研究采用时间标准化，将每位女性的月经周期按比例缩放到28天的参考周期。同时，考虑到激素数据通常呈右偏分布，研究对原始孕酮值进行对数转换，使其更接近正态分布，便于后续统计分析。已经有一些学者利用线性模型分析这组数据，前面的文献旨在研究血清孕酮水平的对数与年龄(Age)、身体质量指数(BMI)和时间($t_{ij}$ )的关系。Fan等人[24]利用Huber得分函数分析这组数据，尽管有界的Huber得分函数可以减弱响应变量中异常值的影响从而得到稳健估计，但是在非稳健估计情形下会损失一些效率。吕晶等人[15]利用修正的Cholesky分解对这组数据进行了有效的秩推断分析，但是文章并没有考虑变量选择时造成的模型损失。考虑到孕酮数据随时间变化的特殊性，本文建立了B样条近似下的变系数秩回归模型，叠加非凸惩罚函数进行变量选择，研究血清孕酮水平的对数与Age、BMI、$t_{ij}$ 以及它们交互项之间的关系，具体模型如下：

$Y_{ij}={\beta }_0+{\beta }_1{\text{Age}}_i+{\beta }_2{\text{BMI}}_i+{\beta }_3{t}_{ij}+{\beta }_4{\text{Age}}_i\ast t_{ij}+{\beta }_5{\text{BMI}}_i\ast t_{ij}+{\beta }_6{t}_{ij}^2+{\beta }_7{\text{Age}}_i\ast {\text{BMI}}_i+{\epsilon }_{ij}$ (21)

下面比较相同样本数据下，本文所提方法(${\widehat{\beta }}_{PB}^{SCAD}$ 、${\widehat{\beta }}_{PB}^{MCP}$ )、PGEE方法(${\widehat{\beta }}_{PGEE}^{CS}$ , ${\widehat{\beta }}_{PGEE}^{indep}$ , ${\widehat{\beta }}_{PGEE}^{AR}$ )以及吕等人[15]提出的方法(${\widehat{\beta }}_C$ )的结果。表2列举了200次Bootstrap抽样得到的估计系数(Estimate)和对应的标准差(SE)。

由表2可得如下结论：

1) 在上述所有方法中，本文提出的方法具有最小标准差，这说明本文提出的方法更有效；

2) 对比${\widehat{\beta }}_{PGEE}^{CS}$ 、${\widehat{\beta }}_{PGEE}^{indep}$ 、${\widehat{\beta }}_{PGEE}^{AR}$ 以及${\widehat{\beta }}_{PB}^{SCAD}$ 、${\widehat{\beta }}_{PB}^{MCP}$ 可知，利用修正的Cholesky分解处理纵向秩回归协方差矩阵的结构是有效的，明显降低了模型估计过程中损失的效率；

3) 在本文提出的方法下，${\widehat{\beta }}_{PB}^{MCP}$ 较${\widehat{\beta }}_{PB}^{SCAD}$ 的标准差更小，说明同样作为非凸惩罚函数，在修正的Cholesky分解下的变系数纵向秩回归模型中，MCP惩罚函数对模型有效性的提高更显著。图1是对200次Bootstrap抽样的一次随机样本数据绘制得到的相关性热图。可知由于交互项的存在，变量$X_3$ 、$X_4$ 、$X_5$ 、$X_6$ 之间出现了强相关性，比较${\widehat{\beta }}_{PB}^{MCP}$ 和${\widehat{\beta }}_{PB}^{SCAD}$ 对应的两种估计方法下得到的四个变量的标准差结果，可以看出，强相关下，MCP惩罚函数表现出更明显的优势；

Table 2. Estimated coefficients (×10⁻²) and standard deviations (×10⁻¹) derived from the progesterone data

表2. 孕酮数据中所得的估计系数(×10⁻²)以及标准差(×10⁻¹)

		Age	BMI	Time	Age*Time	BMI* Time	Time* Time	Age* BMI
${\widehat{\beta }}_C$	Estimate	−5.0348	−5.3870	−15.481	0.2718	0.2670	−0.0201	0.1434
	SE	5.6838	8.2926	1.8319	0.0344	0.0442	0.0243	0.2211
${\widehat{\beta }}_{PGEE}^{CS}$	Estimate	7.3184	11.5133	−4.2705	0.1122	0.0557	−0.0121	−0.2473
	SE	1.4396	2.2622	1.1004	0.0244	0.0313	0.0173	0.0658
${\widehat{\beta }}_{PGEE}^{indep}$	Estimate	−7.2348	−12.870	−3.2181	0.1218	0.3670	−0.0211	0.1134
	SE	5.1661	7.7001	1.4750	0.0321	0.0362	0.0339	0.1058
${\widehat{\beta }}_{PGEE}^{AR}$	Estimate	7.1611	11.6469	−3.8152	0.1128	0.0395	−0.0138	−0.2435
	SE	1.4663	2.2829	1.1037	0.0244	0.0327	0.0176	0.0669
${\widehat{\beta }}_{PB}^{SCAD}$	Estimate	24.9950	39.2499	−3.5169	0.1956	−0.0744	−0.0484	−0.9459
	SE	1.9106	3.0075	6.9158	0.1284	0.1664	0.0711	0.1113
${\widehat{\beta }}_{PB}^{MCP}$	Estimate	0.4673	0.7336	−2.6240	0.0725	0.0848	0.0049	0.1542
	SE	0.6155	1.0663	2.9786	0.0671	0.0650	0.0382	0.0498

Figure 1. Bootstrap resampling-based correlation heatmap

图1. Bootstrap随机抽样样本相关性热图

Figure 2. Boxplot of estimation results from 200 bootstrap resampling trials

图2. 200次Bootstrap随机抽样估计结果的箱线图

4) 绘制上述五种方法的参数估计结果对应箱线图如下，结合图2的特点分析，MCP惩罚函数下的箱线图相对集中，说明估计结果的稳定性和一致性更强，而SCAD惩罚下的箱体较宽，这是因为SCAD惩罚对高相关性变量的处理不稳定，会导致系数被过度压缩或方向不一致。总之，根据箱线图分析可知，对于维度相同的样本数据，MCP惩罚函数对变量的选择更果断，能有效降低估计方差。

6. 结论

本文针对纵向数据中存在的时变效应与异质性结构，提出了一种基于变系数线性模型的稳健估计方法。该方法让每个变量的影响大小可以随时间变化，比如在研究药物疗效时，药物对病人的作用效果可能会随着治疗时间的推移逐渐增强或减弱。相比之下，传统方法只能假设药物的效果在整个治疗过程中保持不变，这往往不符合实际情况。同时，为进一步提升估计效率，本文设计了基于修正Cholesky分解的协方差矩阵估计方法，避免了预设参数结构的局限性，并通过惩罚秩回归技术实现了对离群值与厚尾分布的高鲁棒性。理论证明与数值模拟均表明，所提方法在参数估计的一致性、效率性及稳健性方面均优于现有方法。除了纵向数据，本文提出的方法还可扩展至非平稳时间序列分析与时变比例风险模型等场景，为复杂动态系统的统计建模提供了通用工具。