基于秩回归的多阈值变平面模型研究

Research on Multi-Threshold Variable Plane Model Based on Rank Regression

摘要: 为了更好地应用个性化医疗技术，识别出导致治疗效果出现异质性的亚组人群，进而为这些特定的患者群体提供更加精准有效的治疗方法，文章介绍了一种多阈值的变平面模型，该模型将研究对象划分为具有不同协变量效应的亚组。针对此模型，文章提出了一种新的基于秩回归的两阶段估计方法来确定亚组数量、阈值的位置以及其他的回归参数。在第一阶段，采用分组选择原则，一致地识别亚组的数量；在第二阶段，采用惩罚诱导平滑技术来精确估计变点位置和模型的其他参数。通过将该方法应用到艾滋病临床试验组研究175所得数据中，并对比其他的模型方法，可以发现基于秩回归的估计方法结果更好，同时该方法具有很好的鲁棒性。

Abstract: In order to better apply personalized medical technology, identify subgroups that cause heterogeneity in treatment outcomes, and provide more accurate and effective treatment methods for these specific patient groups, this paper introduces a multi-threshold variable plane model, which divides the research subjects into subgroups with different covariate effects. This article proposes a new two-stage estimation method based on rank regression to determine the number of subgroups, the position of thresholds, and other regression parameters for this model. In the first stage, the principle of grouping selection is adopted to consistently identify the number of subgroups. In the second stage, the penalty-induced smoothing technique is used to accurately estimate the position of the inflection point and other parameters of the model. By applying the method proposed in this paper to the data obtained from the AIDS clinical trial group study 175, compared with other model methods, we can find that the estimation method based on rank regression has better results, and this method has good robustness.

1. 引言

随着大数据时代的到来，对每一个患者的医疗健康记录也越来越详细。通常情况下，不同患者对不同疾病的治疗方案会产生不同的反应，患者之间治疗效果的异质性使医疗决策复杂化，从而促使医疗研究者们从传统的“一刀切”方法过渡到量身定制疗法，即个性化医疗。

个性化医疗的一个新兴研究方向是使用适当的算法技术对患者进行分类，得到不同的亚组，并分析这些亚组的治疗效果。在现有的文献中，已经开发出了多种基于数据驱动的亚组识别方法。Zhao等人[1]引入了结果加权学习(O-Learning)框架，从分类的角度直接寻找最佳的二元处理方法。Loh [2]提出了广义无偏交互检测与估计方法，旨在识别那些在接受治疗后效果得到显著提升的患者亚组。此外，Foster等人[3]设计了一种虚拟双胞胎方法，用于发现治疗效果增强的特定亚组。同时，Cai等人[4]和Zhao等人[5]则采用参数评分体系来估算个体患者的治疗差异，并据此筛选出从新治疗中可能获得更大益处的人群。

事实上，上述大多数亚组识别工作使用了与变点检测相似的技术[6]，并且可以使用传统的变点理论来严格证明。因此，可以考虑基于变点模型提出亚组识别方法，变点模型也称为阈值回归或分段回归模型，它的一个自然扩展是变平面模型。变平面模型相对于变点模型的优点是，它允许潜在的分组变量是协变量的线性组合，而不是局限于单个协变量，可以处理的变量更多。

到目前为止，由于难以处理由变平面定义的多维断点和不连续的亚组指标，对多阈值变平面模型只进行了少量的研究。Li和Jin [7]提出了一种基于惩罚的加速失效时间回归模型框架。在Jin等人[8]的研究中，他们将阈值问题定义为群体模型选择问题，并应用了快速计算工具。Wei和Kosorok [9]使用切片逆回归技术提出了协变量空间中具有变平面的不规则Cox模型，Fan等人[10]使用双鲁棒性评分统计量考虑了使用变平面方法来检验亚组是否存在，但是该方法只允许一个阈值(即只有两个子组)，并且在一个单位球上搜索平方分数测试统计数据的最大值相当具有挑战性，尤其是在针对多个组时。为了解决上述问题，本文考虑了一个具有未知阈值的变平面模型，该模型允许存在多个变平面，自动生成具有不同协变量效应的亚组，从而促进个性化医疗在医学领域的应用[11]。

本文提出的方法其优势主要体现在以下三个方面。首先，本文研究的变平面的概念不只是使用预先指定的索引变量，而且允许协变量的线性组合，可以得到更有意义的亚组定义，进而为精准医疗提供更灵活的工具。其次，本文对普通的基于最小二乘的多阈值变平面模型估计进行了重要改进，考虑基于秩回归的多阈值变平面模型估计，将阈值识别问题转化为模型选择问题，所得结果在鲁棒性方面具有更显著的优势。此外，在实践中，亚组可能只是在少数选定的协变量效应上有所不同，在其他协变量上具有相同的效应。因此，本文允许某些效应为零，从而获得稀疏解，这是通过惩罚诱导平滑估计方法实现的。

2. 多阈值变平面

2.1. 多阈值变平面模型基本理论

假设$\left\{{\left(Y_i,X_i^{\text{T}},Z_i^{\text{T}}\right)}^{\text{T}},i=1,\cdots ,n\right\}$ 遵循以下具有s个阈值的变平面模型：

$Y_i=X_i^{\text{T}}\left[\beta +{\displaystyle \sum _{j=1}^s{\delta }_j}I\left(Z_i^{\text{T}}\theta >a_j\right)\right]+{\epsilon }_i,i=1,\cdots ,n$ (1)

其中$X={\left(X_{i1},\cdots ,X_{ip}\right)}^{\text{T}}$ 是p维向量，$Y_i\in ℝ$ 是与第i个研究对象$X_i$ 对应的响应变量。$\beta$ 为基线组的协变量效应向量，${\delta }_j$ 为第j个亚组相对于基线组的增强效应向量，其子组由变平面$I\left(Z_i^{\text{T}}\theta >a_j\right)$ 定义。当$\delta =0$ 时，参数$\theta$ 不被识别。$Z_i$ 是不包含截距项的d维分组变量，$\theta$ 为变平面参数，$a_1,\cdots ,a_s$ 是阈值位置，且满足$a_1<a_2<\cdots <a_s$ 。${\epsilon }_i$ 是由${\epsilon }_{ij}$ 组成的独立的误差项，$i=1,2,\cdots ,n$ ，$1\le j\le p$ ，且${\epsilon }_i$ 的概率密度函数均为$f(\cdot )$ 。${\epsilon }_{ij}$ 是连续的随机变量，且对于任何成对的差异值${\epsilon }_{ik}-{\epsilon }_{jl}$ ，$i,j=1,\cdots ,n$ ，$j,l=1,\cdots ,p$ ，其中位数均为0。在这个模型中，可以通过将不同的治疗组编码为$X_j$ 中的虚拟协变量来合并多个治疗效果，相应的系数$\beta$ 和${\delta }_j$ 表示不同处理方案的效果。为了更精确地识别模型，本文考虑对$\theta$ 做出假设：$\theta \in \Theta =\left\{\theta \in {ℝ}^d:\Vert \theta \Vert =1,{\theta }_r>0,1\le r\le d\right\}$ ，且固定的第r个元素为正。

2.2. 基于秩回归的多阈值变平面模型求解

Hettmansperger [12]指出基于秩回归的估计方法是不依赖于总体分布的、具有稳健性的、高度有效的。此外，相对于一般的最小二乘估计法，秩回归方法能够有效地处理异方差和异常值问题。

记$\eta ={\left({\beta }^{\text{T}},{\delta }^{\text{T}},a^{\text{T}},{\theta }^{\text{T}}\right)}^{\text{T}}$ 。根据Jaeckel [13]的思想，结合秩回归的理论知识可知，秩估计的目的是最小化残差的离散度。因此，当$a_j,j=1,\cdots ,s$ 已知时，未知参数$\beta ,\delta ,\theta$ 可以通过最小化，如下目标函数来估计：

$L_n\left(\eta \right)=\frac{1}{n}{\displaystyle \sum _{i<j}\left|e_i-e_j\right|}$ (2)

其中$e_i=Y_i-X_i^{\text{T}}\beta -X_i^{\text{T}}{\displaystyle \sum _{j=1}^s{\delta }_j}I\left(Z_i^{\text{T}}\theta >a_j\right)$ ，并约束$\theta \in \Theta$ 。

然而，在一般情况下，变平面的数量和位置都是未知的，而且在回归模型中，需要在选择出显著协变量的同时及时更新参数估计值。显然，公式(2)并不能同时进行显著变量选择和参数估计，因此考虑使用惩罚估计。在实践中，本文提出了一种迭代的两阶段多阈值变平面估计方法。在第一阶段中，对于任意给定的一致估计量$\widehat{\theta }$ ，可通过带有惩罚的变点检测算法得到s的一致估计量$\widehat{s}$ 。得到$\widehat{s}$ 之后，在第二阶段使用诱导平滑方法来估计参数$\left(\beta ,\delta ,\theta ,a\right)$ 。具体过程如下。

2.2.1. 初步分割阶段

对于给定的$\widehat{\theta }$ ，记${\widehat{W}}_i=Z_i^{\text{T}}\widehat{\theta },i=1,\cdots ,n$，生成秩映射$\left\{l_{\left(i\right)}:1\le i\le n\right\}$ ，使得${\widehat{W}}_{l_{\left(i\right)}}$ 是$\left\{{\widehat{W}}_i:1\le i\le n\right\}$ 中的第i个最小值，为了避免位置估计不一致，将其按升序排列，即${\widehat{W}}_{l_{\left(1\right)}}\le {\widehat{W}}_{l_{\left(2\right)}}\le \cdots \le {\widehat{W}}_{l_{\left(n\right)}}$ 。

首先，本文根据${\widehat{W}}_{l_{\left(i\right)}}$ 将数据序列分为$q_n+1$ 段，其中当$n\to \infty$ 时，$q_n$ 趋于无穷大。即拆分数据序列使得第一段${ℱ}_1=\left\{i:{\widehat{W}}_i\le {\widehat{W}}_{l_{\left(n-q_{n}m\right)}}\right\}$ ，包含第$n-q_{n}m$ 个观测值，其他$q_n$ 段${ℱ}_j=\left\{i:{\widehat{W}}_{l_{\left(n-\left(q_n-j+2\right)m\right)}}\le {\widehat{W}}_i\le {\widehat{W}}_{l_{\left(n-\left(q_n-j+1\right)m\right)}}\right\}$ ，$j=2,\cdots ,q_n+1$ ，分别包含m个观测值，其中$m=\lceil n/q_n\rceil$ 。定义$b_j=\left|{ℱ}_j\right|,j=1,\cdots ,q_n+1$ ，表示集合${ℱ}_j$ 中所含元素个数，即集合${ℱ}_j$ 的大小。令$Y_{\left(j\right)}={\left(Y_i,i\in {ℱ}_j\right)}^{\text{T}}$ ，$X_{\left(j\right)}={\left(X_i,i\in {ℱ}_j\right)}^{\text{T}}$ ，$\tilde{Y}={\left(Y_{\left(1\right)}^{\text{T}},\cdots ,Y_{\left(q_n+1\right)}^{\text{T}}\right)}^{\text{T}}$，$\tilde{X}=\left(X^{\left(1\right)},\cdots ,X^{\left(q_n+1\right)}\right)$ ，$X^{\left(1\right)}={\left(X_{\left(1\right)}^{\text{T}},\cdots ,X_{\left(q_n+1\right)}^{\text{T}}\right)}^{\text{T}}$ ，$X^{\left(j\right)}={\left(0_{p\times {\displaystyle {\sum }_{i=1}^{j-1}{b}_j}},X_{\left(j\right)}^{\text{T}},\cdots ,X_{\left(q_n+1\right)}^{\text{T}}\right)}^{\text{T}}$ ，$j=2,\cdots ,q_n+1$ ，则估计值${\tilde{\gamma }}^{\prime }={\left({\tilde{\beta }}^{\text{T}},{{\tilde{\delta }}^{\prime }}_1^{\text{T}},\cdots ,{{\tilde{\delta }}^{\prime }}_{q_n}^{\text{T}}\right)}^{\text{T}}$可以通过最小化以下带有惩罚的目标函数得到：

${L^{\prime }}_n\left({\tilde{\gamma }}^{\prime }\right)=\frac{1}{n}{\displaystyle \sum _{i<j}\left|{e^{\prime }}_i-{e^{\prime }}_j\right|}+{\displaystyle \sum _{j=1}^{q_n}{p}_{{\lambda }_n}}\left(\Vert {{\delta }^{\prime }}_j\Vert \right)$ (3)

其中${e^{\prime }}_i=Y_{\left(i\right)}-X^{\left(i\right)}{\gamma }_i$ ，$p_{{\lambda }_n}(\cdot )$ 为惩罚函数。

本文选取Fan和Li [14]提出的具有较好性质的非凹惩罚函数，即平滑截断的绝对偏差(SCAD)惩罚，其定义为其一阶导数，并且是关于原点对称的。即对于$\theta >0$ ，

${p^{\prime }}_{{\lambda }_n}\left(\theta \right)=\lambda \left\{I\left(\theta \le {\lambda }_n\right)+\frac{{\left(a\lambda -\theta \right)}_{+}}{\left(a-1\right)\theta }I\left(\theta >{\lambda }_n\right)\right\}$

其中$a>2$ ，且${\lambda }_n$ 是一个非负惩罚参数，用来控制模型的变量选择或稀疏性。一般地，当p较大时，通常假设$\gamma$ 为稀疏结构。现考虑使用迭代局部二次近似算法来找到式(3)的最小值。Leng [15]指出式(2)中的目标函数可以视为Jaeckel [13]提出的Wilcoxon分数的秩离散函数，所以式(3)中的第一项可以近似为：

$\frac{1}{n}{\displaystyle \sum _{i<j}\left|{e^{\prime }}_i-{e^{\prime }}_j\right|}\approx {\displaystyle \sum _{i=1}^n{\omega }_i{\left({e^{\prime }}_i-\xi \right)}^2}$

其中$\xi$ 是${\left\{{e^{\prime }}_i\right\}}_{i=1}^n$ 的中位数，且

${\omega }_i=\{\begin{array}{l}\frac{\frac{R\left({e^{\prime }}_i\right)}{n+1}-\frac{1}{2}}{{e^{\prime }}_i-\xi },{e^{\prime }}_i\ne \xi \\ 0,其他\end{array}$

使用简单的泰勒展开，给定目标函数(3)的初始估计${\gamma }^{\prime }$ (等效于给定了$\beta$ 和${{\delta }^{\prime }}_j$ )，可以相应地获得权重${\omega }_i$ 和残差的中位数$\xi$ 。对于式(3)的第二项，本文考虑近似正则项为：

$p_{{\lambda }_n}\left(\Vert {{\delta }^{\prime }}_j\Vert \right)\approx p_{{\lambda }_n}\left(\Vert {\widetilde{{\delta }^{\prime }}}_j\Vert \right)+\frac{1}{2}\frac{{p^{\prime }}_{{\lambda }_n}\left(\Vert {\widetilde{{\delta }^{\prime }}}_j\Vert \right)}{\Vert {\widetilde{{\delta }^{\prime }}}_j\Vert }\left\{{\Vert {{\delta }^{\prime }}_j\Vert }^2-{\Vert {\widetilde{{\delta }^{\prime }}}_j\Vert }^2\right\}$

因此，式(3)可以近似为：

${L^{″}}_n\left({\tilde{\gamma }}^{\prime }\right)\approx {\left(\tilde{Y}-\tilde{X}{\tilde{\gamma }}^{\prime }\right)}^{\text{T}}\tilde{W}\left(\tilde{Y}-\tilde{X}{\tilde{\gamma }}^{\prime }\right)+\frac{n^2}{2}{{\tilde{\gamma }}^{\prime }}^{\text{T}}\Omega {\tilde{\gamma }}^{\prime }$ (4)

其中$\tilde{Y}={\left(Y_{\left(1\right)}^{\text{T}},\cdots ,Y_{\left(q_n+1\right)}^{\text{T}}\right)}^{\text{T}}$，且

$\Omega =diag\left(\frac{{p^{\prime }}_{{\lambda }_n}\left(\Vert {\widetilde{{\delta }^{\prime }}}_1\Vert \right)}{\Vert {\widetilde{{\delta }^{\prime }}}_1\Vert },\cdots ,\frac{{p^{\prime }}_{{\lambda }_n}\left(\Vert {\widetilde{{\delta }^{\prime }}}_{q_n}\Vert \right)}{\Vert {\widetilde{{\delta }^{\prime }}}_{q_n}\Vert }\right)$ ，$\tilde{W}=diag\left({\tilde{\omega }}_1,\cdots ,{\tilde{\omega }}_n\right)$

为方便起见，记估计量${\tilde{\gamma }}^{\prime }=\left({{\tilde{\gamma }}^{\prime }}_1^{\text{T}},\cdots ,{{\tilde{\gamma }}^{\prime }}_{q_n+1}^{\text{T}}\right)$，$\widehat{\mathcal{A}}=\left\{j:{{\tilde{\gamma }}^{\prime }}_j\ne 0,j=1,\cdots ,q_n+1\right\}$ ，

${\widehat{\mathcal{A}}}^{\ast }=\left\{j:j\in \widehat{\mathcal{A}},j-1\notin \widehat{\mathcal{A}},j=2,\cdots ,q_n+1\right\}=\left\{{\widehat{\kappa }}_1,\cdots ,{\widehat{\kappa }}_{\widehat{s}}\right\},{\widehat{\kappa }}_1<\cdots <{\widehat{\kappa }}_{\widehat{s}}$ (5)

显然，由上述定义可知${\widehat{\mathcal{A}}}^{\ast }$ 是$\widehat{\mathcal{A}}$ 的子集。若$j-1\notin \widehat{\mathcal{A}},j\in \widehat{\mathcal{A}}$ 和$j+1\in \widehat{\mathcal{A}}$ ，则$j\in {\widehat{\mathcal{A}}}^{\ast }$ 且$j+1\notin {\widehat{\mathcal{A}}}^{\ast }$ 。因此，当给定每个估计量$\widehat{\theta }$ 的值，就可以获得变平面的阈值数量估计值$\widehat{s}=\left|{\widehat{\mathcal{A}}}^{\ast }\right|$ 。如果$\widehat{s}=0$ ，则表明没有子组。如果$\widehat{s}>0$ ，则真阈值$a_j$ 极有可能位于区间$\left({\widehat{W}}_{l_{\left(n-\left(q_n-{\widehat{\kappa }}_j+3\right)m\right)}},{\widehat{W}}_{l_{\left(n-\left(q_n-{\widehat{\kappa }}_j+1\right)m\right)}}\right),j=1,\cdots ,\widehat{s}$ 。若给定的估计量$\widehat{\theta }$ 是一致的，则初步分割阶段得到的估计值$\widehat{s}$ 也将以很高的概率收敛。

2.2.2. 平滑优化阶段

根据初步分割阶段得到的变平面阈值数量的估计值$\widehat{s}$ ，可通过最小化以下平滑目标函数来估计模型中的参数$a,\theta$ 和$\gamma ={\left({\beta }^{\text{T}},{\delta }_1^{\text{T}},\cdots ,{\delta }_{\widehat{s}}^{\text{T}}\right)}^{\text{T}}$：

$\frac{1}{n}{\displaystyle \sum _{i<j}\left|e_i-e_j\right|}$ (6)

其中$e_i=Y_i-X_i^{\text{T}}\beta -X_i^{\text{T}}{\displaystyle \sum _{k=1}^{\widehat{s}}{\delta }_{k}I\left(\frac{Z_i^{\text{T}}\theta -a_k}{h}\right)}$。

针对以上目标函数，本文考虑使用一个迭代估计过程，它可以产生相对稳定的解。对于给定的$\theta$ ，目标函数(6)可以简单地视为分段函数，可以通过秩回归方法来估计基线系数$\beta$ 和增强效应$\delta$ 。然而，对于给定的$\gamma$ ，目标函数(6)不是连续的，找到它的最小值比较困难。因此，本文考虑使用光滑的分布函数$\Phi \left(\cdot /h\right)$ 作为不连续的示性函数的平滑逼近，其中带宽h随着样本量n的增加逐渐趋于0。如果$Z_i^{\text{T}}\theta >a_k$ ，则当$h\to 0$ 时，$\Phi \left(\left(Z_i^{\text{T}}\theta -a_k\right)/h\right)\to 1$ 。因此，目标函数(6)可以用如下函数来近似：

$Q_n^{\ast }\left(\eta \right)=\frac{1}{n}{\displaystyle \sum _{i<j}\left|e_i^{\ast }-e_j^{\ast }\right|}$ (7)

其中$e_i^{\ast }=Y_i-X_i^{\text{T}}\beta -X_i^{\text{T}}{\displaystyle \sum _{k=1}^{\widehat{s}}{\delta }_k\Phi \left(\frac{Z_i^{\text{T}}\theta -a_k}{h}\right)}$。

记${\tilde{\eta }}^{\ast }=\arg {\min }_{\theta \in \Theta }\left\{Q_n^{\ast }\left(\eta \right)\right\}$ 。对于非稀疏问题，可以直接通过Newton-Type算法来最小化目标函数(7)。对于稀疏问题，考虑在式(7)中增加惩罚函数，即最小化以下惩罚目标函数来估计$\eta$ ：

$Q_n^{\ast \ast }\left(\eta \right)=\frac{1}{n}{\displaystyle \sum _{i<j}\left|e_i^{\ast }-e_j^{\ast }\right|}+{\displaystyle \sum _{j=1}^{q_n+1}{p}_{{\lambda }_n}\left(\Vert {\gamma }_j\Vert \right)}$ (8)

记${\tilde{\eta }}^{\ast \ast }=\arg {\min }_{\theta \in \Theta }\left\{Q_n^{\ast \ast }\left(\eta \right)\right\}$ ，可通过迭代诱导平滑方法来处理上述带有惩罚的目标函数，进而得到参数的估计值。具体的计算类似初步分割阶段的处理方法，对式(8)的第一项和第二项进行近似即可。

重复进行上述初步分割阶段和平滑优化阶段，直到满足收敛条件。尤其当变平面的估计数量保持不变时，考虑终止迭代。现将本文提出的求解基于秩回归的多阈值变平面模型的算法记为TSMCPLD，其具体算法过程如下。

算法. TSMCPLD

步骤0：给定$\theta$ 的初始估计，记为${\tilde{\theta }}_{int}^{\ast }$ ，并设置${\tilde{\theta }}^{\ast }={\tilde{\theta }}_{int}^{\ast }/\Vert {\tilde{\theta }}_{int}^{\ast }\Vert$ ；

步骤1：实现初步分割阶段。最小化式(4)得到估计值${\tilde{\gamma }}^{\ast }$ ，然后计算式(5)中定义的索引集${\widehat{\mathcal{A}}}^{\ast }$ ，通过$\widehat{s}=\left|{\widehat{\mathcal{A}}}^{\ast }\right|$ 获得阈值数量；

步骤2：对于给定的$\widehat{s}$ ，在平滑优化阶段中，通过最小化惩罚目标函数(8)来更新参数$\left({\widehat{\theta }}^{\ast },{\widehat{a}}^{\ast },{\widehat{\gamma }}^{\ast }\right)$ ；

步骤3：迭代步骤1和2直至收敛。

在上述算法中，初步分割阶段的性能依赖于片段长度m，并且最佳m的选择可以遵循[7]中的建议，将数据序列进行L次初步分割阶段，事件(不包括第一段)的公共长度为$m_{\ell },\ell =1,\cdots ,L$ ，设$m_{\ell }=\lfloor {\kappa }_{\ell }\sqrt{n}\rfloor ,\ell =1,\cdots ,L$ ，其中${\kappa }_{\ell }$ 取区间内$\ell$ 个网格点的值，贝叶斯信息准则(BIC)可用于选择最佳索引$\ell$ 。

在平滑优化阶段的具体计算中，可以使用R包BB中的BBoptim函数来优化高维非线性目标函数。此外，参数的数量可能相当大，尤其是有大量的异构子组时，包含惩罚函数会导致稀疏解。在中等或高维环境下，调谐参数${\lambda }_n$ 也可以通过贝叶斯信息准则(BIC)来选择。具体的计算公式为：

$\text{BIC}\left({\lambda }_n\right)=n\log \left(\text{RSS}/n\right)+\log \left(n\right)\left(p+1\right)\left(\widehat{s}+1\right)$

其中RSS为残差平方和。

3. 渐进理论性质

为了更好地建立渐近理论性质，本文考虑增加以下必要条件：

条件1：a) $E\left(X_i{X}_i^{\text{T}}\right)={\Sigma }_0$ 是有限且正定的。$E\left(Z_i{Z}_i^{\text{T}}\right)$ 是正定的。$Z_i$ 和${\epsilon }_i$ 是独立的，$i=1,\cdots ,n$ 。$E\left({\epsilon }_i|X_i\right)=0$ 几乎必然成立。

b) 对于某些$\xi >1$ ，设$0<E{\Vert {\left(X_i^{\text{T}},Z_i^{\text{T}}\right)}^{\text{T}}\left(X_i^{\text{T}},Z_i^{\text{T}}\right)\Vert }^{\xi }<\infty$ ，且$E{\Vert {\left(X_i^{\text{T}},Z_i^{\text{T}}\right)}^{\text{T}}{\epsilon }_i\Vert }^{\xi }<\infty$ ，$E\left(X_i{X}_i^{\text{T}}|Z_i\right)>0$ 。

条件2：参数$\eta$ 的参数空间是紧空间，其中${\min }_{1\le l\le k\le s}\left\{\left|a_l-a_k\right|\right\}$ 且$\min \left\{\Vert \beta /\sqrt{p}\Vert ,\Vert {\delta }_1/\sqrt{p}\Vert ,\cdots ,\Vert {\delta }_s/\sqrt{p}\Vert \right\}$ 以零为界。设$\rho \left(t\right)={\lambda }_n^{-1}{p}_{{\lambda }_n}\left(t\right)$ 且$\overline{\rho }\left(t\right)={\rho }^{\prime }\left(\left|t\right|\right)\mathrm{sgn}\left(t\right)$ 。本文假设罚函数$p_{{\lambda }_n}(\cdot )$ 满足以下条件：

条件3：$p_{{\lambda }_n}(\cdot )$ 是一个对称函数，它在$\left[0,\infty \right)$ 上是不递减的凹函数。存在一个常数$v>0$ ，使得$\rho \left(t\right)$ 对于所有$\left|t\right|\ge v{\lambda }_n$ 都是常数，且$\rho \left(0\right)=0$ 。除了有限的t，${\rho }^{\prime }\left(t\right)$ 存在且连续，而且${\rho }^{\prime }\left(0+\right)=1$ 。本文设$a_0=-\infty$ 且$a_{s+1}=\infty$ 。记${ℱ}_j^0=\left\{i:a_{j-1}<Z_i^{\text{T}}\theta <a_j\right\},j=1,\cdots ,s$ ，其包含真实向量阈值位置$a$ ，变平面参数$\theta$ 。与本文中$\tilde{X}$ 和$\tilde{Y}$ 的定义类似，可以通过将${ℱ}_j^0$ 替换为${ℱ}_j$ 来定义${\tilde{X}}_a$ 和${\tilde{Y}}_a$ 。根据条件1可知，${\tilde{X}}_a^{\text{T}}{\tilde{X}}_a/n\stackrel{a.s.}{\to }\Upsilon$，其中$\Upsilon$ 是一个正定矩阵。为了获得式(8)中${\widehat{\gamma }}^{\ast }$ 的渐近性质，考虑以下假设条件：

条件4：${\max }_{u\ge 0}\left\{{p^{″}}_{{\lambda }_n}\left(u\right)\right\}+{\Lambda }_{\left(s+1\right)p}\left(\Upsilon \right)>0$ ，其中${\Lambda }_{\left(s+1\right)p}\left(\Upsilon \right)$ 是$\Upsilon$ 的最小特征值。

条件5：设置$W_i=Z_i^{\text{T}}\theta$ 且$V_i=\left(X_i,Z_i\right)$ ，$f_{W|V}(\cdot )$ 表示给定$V_i=V$ 时$W_i$ 的条件密度，$f_W(\cdot )$ 表示$W_i$ 的密度，其中$f_{W|V}(\cdot )$ 有紧支撑且二阶导数有界。$P\left(W_i\le a_j\right)={\tau }_j,0<{\tau }_1<\cdots <{\tau }_s<1$ 。此外，对于某些$M<\infty$ ，$E\left({\epsilon }_i^4|V_i\right)<M$ 几乎必然成立。

条件6：$h\to 0$ 和$n{h}^2\to 0$ ，当$n\to \infty$ 时。

设计矩阵的条件1允许特定制度的异方差性。误差假设可以被放宽为${\epsilon }_i=\sigma \left(X_i^{\text{T}}\beta \right)e_i$ ，其中$e_i$ 独立于$X_i$ ，且$e_1,\cdots ,e_n$ 是i.i.d.的，均值为零和方差为${\sigma }^2$ 。条件2是关于参数空间的，它通过要求$a_{j-1}<a_j,j=1,\cdots ,s$ 来排除小于$s+1$ 个亚组的简化模型的可能性。条件3和4通常用于高维数据环境的收缩回归中。SCAD这类凹形惩罚满足条件3，对于SCAD惩罚，条件4等价于${\Lambda }_{\left(s+1\right)p}\left(\Upsilon \right)>1/\left(v-1\right)$ ，这确保了目标函数(8)是全局凸的。条件5是一个标准的平滑条件，其表明存在s个不同的跳跃，否则模型无法识别。条件6可以来确定h的速率。

根据大数定律和条件5，可以得到${\displaystyle \sum _{i=1}^{n}I\left(a_{j-1}<W_i\le a_j\right)}/n\stackrel{p}{\to }{\tau }_j-{\tau }_{j-1}>0$ ，且当$m\to \infty$ 和$m/n\to 0$ 时，有${\displaystyle \sum _{i=1}^{n}I\left({\widehat{W}}_{l_{\left(n-\left(q_n-j+2\right)m\right)}}<{\widehat{W}}_i\le {\widehat{W}}_{l_{\left(n-\left(q_n-j+1\right)m\right)}}\right)}/n=m/n\to 0$ 。因此，当$\theta$ 是已知的或是一致估计值，即${\widehat{W}}_i=W_i+o_p\left(1\right)$ 时，对于足够大的n，在每个分段${ℱ}_j=\left\{{\widehat{W}}_{l_{\left(n-\left(q_n-j+2\right)m\right)}}<{\widehat{W}}_i\le {\widehat{W}}_{l_{\left(n-\left(q_n-j+1\right)m\right)}}\right\}$ 中，最多有一个阈值，其中，${\widehat{W}}_{l_{\left(n-\left(q_n-j+2\right)m\right)}}$ 和${\widehat{W}}_{l_{\left(n-\left(q_n-j+1\right)m\right)}},j=1,\cdots ,q_n+1$ 的定义在第2.2.1节中。由以下定理1即可保证初步分割阶段得到变平面数量估计值$\widehat{s}$ 是s的一致估计。

定理1. 假设$m\to \infty$ 且$m=O\left(n^r\right)$ ，其中$0<r\le 1/2$ 为常数。当$n\to \infty$ 时，${\lambda }_n\to 0$ 和${\lambda }_n\sqrt{n}/\log n\to \infty$ 。如果条件1~5成立，那么有$\underset{n\to \infty }{\lim }P\left(\widehat{s}=s\right)=1$ 。

设$\gamma ={\left({\gamma }_1,\cdots ,{\gamma }_{\left(s+1\right)p}\right)}^{\text{T}}={\left({\beta }^{\text{T}},{\delta }_1^{\text{T}},\cdots ,{\delta }_s^{\text{T}}\right)}^{\text{T}}$ 是(1)中的回归参数，$\mathcal{G}=\left\{j:{\gamma }_j\ne 0,j=1,\cdots ,\left(s+1\right)p\right\}$ 是模型中

的重要变量集。对于给定的一致估计量$\widehat{s}$ ，由定理1可以获得最小化非正则目标函数(7)的平滑秩回归估计量${\tilde{\eta }}^{\ast }$ 的一致性，其中$s=1$ 。本文考虑通过最小化惩罚平滑目标函数(8)来获得估计量${\widehat{\eta }}^{\ast }$ ，以下定理保证了我们估计量的一致性。

定理2. 在条件1~6下，当$n\to \infty$ 时，$\widehat{s}=s$ 且${\lambda }_n\to 0$ ，$L_n^{\ast }\left(\eta \right)$ 有局部最小值${\widehat{\eta }}^{\ast }$ ，使得$\Vert {\widehat{\gamma }}^{\ast }-\gamma \Vert =O_p\left(\sqrt{1/n}\right)$ ，$\Vert {\widehat{a}}^{\ast }-a\Vert =O_p\left(\sqrt{h/n}\right)$ ，且$\Vert {\widehat{\theta }}^{\ast }-\theta \Vert =O_p\left(\sqrt{h/n}\right)$ ，其中$\Vert {\widehat{\theta }}^{\ast }\Vert =\Vert \theta \Vert =1$ 。

现重新记$\mathcal{G}=\left\{g_1,\cdots ,g_{s+1}\right\}$ ，其中$g_{j+1}=\left\{j_1,\cdots ,j_{p_j}\right\}$ 是第j个子群中$p_j$ 个非零协变量集的索引集，$j=0,1,\cdots ,s$ 。不失一般性，考虑将$\gamma$ 记为${\gamma }_p={\left({\gamma }_{\left(1\right)}^{\text{T}},{\gamma }_{\left(2\right)}^{\text{T}}\right)}^{\text{T}}$ ，其中${\gamma }_{\left(1\right)}=\left({\gamma }_{g_1}^{\text{T}},\cdots ,{\gamma }_{g_{s+1}}^{\text{T}}\right)$ ，${\gamma }_{g_{j+1}}={\left({\gamma }_{j_1},\cdots ,{\gamma }_{j_{p_j}}\right)}^{\text{T}}$ 且${\gamma }_{\left(2\right)}=0$ 。记$X_{i,g_{j+1}}={\left(X_{i,j_1},\cdots ,X_{i,j_{p_j}}\right)}^{\text{T}},j=0,1,\cdots ,s$ 。记${\Sigma }_1={\left({\sigma }_{1,jk}\right)}_{0\le j,k\le s}$ 表示${\displaystyle \sum _{j=0}^s{p}_j}\times {\displaystyle \sum _{j=0}^s{p}_j}$ 个分块矩阵，其中每一块矩阵分别为${\sigma }_{1,jk}=4{\sigma }^{2}E{X}_{i,g_{j+1}}{X}_{i,g_{k+1}}^{\text{T}}I\left(Z_i^{\text{T}}\theta >a_j\vee a_k\right)$ ，${\Sigma }_2=\frac{4}{nh}diag\left\{\frac{{\sigma }^2}{2\sqrt{\pi }}{A}_j+\Pi \cdot B_j,j=2,\cdots ,s\right\}$ 是$s\times s$ 维对角矩阵，其中$A_j=E\left\{{\left({\delta }_j^{\text{T}}{X}_i\right)}^2|Z_i^{\text{T}}\theta =a_j\right\}{f}_W\left(a_j\right)$ ，$B_j=E\left\{{\left({\delta }_j^{\text{T}}{X}_i\right)}^4|Z_i^{\text{T}}\theta =a_j\right\}{f}_W\left(a_j\right)$ ，且$\Pi ={\displaystyle {\int }_{-\infty }^{\infty }\varphi {\left(s\right)}^2{\left(I\left(s>0\right)-\Phi \left(s\right)\right)}^2\text{d}s}$ ，而且${\Sigma }_3=\frac{4}{nh}{\displaystyle \sum _{j=1}^s\left\{\frac{{\sigma }^2}{2\sqrt{\pi }}{G}_j+\Pi \cdot H_j\right\}}$ 是$d\times d$ 维矩阵，其中$G_j=E\left\{Z_i{Z}_i^{\text{T}}{\left({\delta }_j^{\text{T}}{X}_i\right)}^2|Z_i^{\text{T}}\theta =a_j\right\}{f}_W\left(a_j\right)$ ，$H_j=E\left\{Z_i{Z}_i^{\text{T}}{\left({\delta }_j^{\text{T}}{X}_i\right)}^4|Z_i^{\text{T}}\theta =a_j\right\}{f}_W\left(a_j\right)$ 。

设$V_{11}={\left(v_{1k,l}\right)}_{0\le k,l\le s}$ ，其中$v_{1k,l}=2E{X}_{i,g_{k+1}}{X}_{i,g_{k+1}}^{\text{T}}I\left(Z_i^{\text{T}}\theta >a_k\vee a_l\right)$ ，$V_{22}=diag\left(\frac{A_k}{\sqrt{\pi }},k=1,\cdots ,s\right)$ ， $V_{23}={\left(v_{23vk}\right)}_{1\le v\le d,1\le k\le s}$ 且$v_{23vk}=-\frac{1}{\sqrt{\pi }}E\left[Z_{iv}{\left({\delta }_{g_{k+1}}^{\text{T}}{X}_{i,g_{k+1}}\right)}^2|Z_i^{\text{T}}\theta =a_k\right]f_W\left(a_k\right)$ ，$V_{33}=\frac{1}{\sqrt{\pi }}{\displaystyle \sum _{j=1}^s{G}_j}$ 。记：

${\Gamma }_{{\lambda }_n}=diag\left\{{p^{″}}_{{\lambda }_n}\left(\left|{\gamma }_{0_1}\right|\right),\cdots ,{p^{″}}_{{\lambda }_n}\left(\left|{\gamma }_{0_{p_0}}\right|\right),\cdots ,{p^{″}}_{{\lambda }_n}\left(\left|{\gamma }_{s_1}\right|\right),\cdots ,{p^{″}}_{{\lambda }_n}\left(\left|{\gamma }_{s_{p_s}}\right|\right)\right\}$

以下定理中可以证明估计量的极限分布是成立的。

定理3. 在条件1~6下，当$n\to \infty$ 时，${\lambda }_n\to 0$ 和${\lambda }_n\sqrt{n}/\log n\to \infty$ ，在概率趋近于1的情况下，定理2中的惩罚平滑估计量${\widehat{\eta }}^{\ast }=\left({\widehat{\gamma }}_{\left(1\right)}^{\ast },{\widehat{\gamma }}_{\left(2\right)}^{\ast },{\widehat{a}}^{\ast },{\widehat{\theta }}^{\ast }\right)$ 满足：

a) 稀疏性：${\widehat{\gamma }}_{\left(2\right)}^{\ast }\to 0_{}a.s.$

b) 渐近正态性：

$\sqrt{n}\left({\widehat{\gamma }}_{\left(1\right)}^{\ast }-{\gamma }_{\left(1\right)}\right)\stackrel{D}{\to }N\left(0,{\left(V_{11}+{\Gamma }_{{\lambda }_N}\right)}^{-1}{\Sigma }_1{\left(V_{11}+{\Gamma }_{{\lambda }_n}\right)}^{-1}\right)$

$\sqrt{n/h}\tilde{V}\left(\begin{array}{c}{\widehat{a}}^{\ast }-a\\ {\widehat{\theta }}^{\ast }-\theta \end{array}\right)\stackrel{D}{\to }N\left(0,\Omega \right)$

其中$\tilde{V}=\left(\begin{array}{cc}{V}_{22}& V_{23}\\ P_{\theta }{V}_{23}^{\text{T}}& P_{\theta }{V}_{33}\end{array}\right)$，$\Omega =diag\left({\Sigma }_2,P_{\theta }{\Sigma }_3\right)$ 且$p_{\theta }=1-\theta {\theta }^{\text{T}}$ 。此外，$\sqrt{n/h}\tilde{V}\left(\begin{array}{c}{\widehat{a}}^{\ast }-a\\ {\widehat{\theta }}^{\ast }-\theta \end{array}\right)$ 和$\sqrt{n}\left({\widehat{\gamma }}_{\left(1\right)}^{\ast }-{\gamma }_{\left(1\right)}\right)$ 都是渐进独立的。

定理3可以为许多比我们更简单但文献中尚未研究的模型提供推理工具。例如，考虑一维阈值变量(即$d=1$ )的情况很有趣，其中$\theta =1$ ，$Z_i^{\text{T}}\theta =Z_i$ 。然后，可以通过本文提出的估计方法估计${\widehat{\eta }}^{\ast }$ ，并在以下推论中得到所得估计量的分布理论：

推论1. 假设条件1~6成立，那么有$\underset{n\to \infty }{\lim }P\left(\widehat{s}=s\right)=1$ ，此外，$\sqrt{n}\left({\gamma }^{\ast }-\gamma \right)$ 和$\sqrt{n/h}\left({\widehat{a}}^{\ast }-a\right)$ 是渐近独立的，且有：

$\sqrt{n}\left({\widehat{\gamma }}_{\left(1\right)}^{\ast }-{\gamma }_{\left(1\right)}\right)\stackrel{D}{\to }N\left(0,{\left(V_{11}-{\Gamma }_{{\lambda }_N}\right)}^{-1}{\Sigma }_1{\left(V_{11}-{\Gamma }_{{\lambda }_n}\right)}^{-1}\right)$

$\sqrt{n/h}\left({\widehat{a}}^{\ast }-a\right)\stackrel{D}{\to }N\left(0,V_{22}^{-1}{\Sigma }_2{V}_{22}^{-1}\right)$

4. 实证分析

现考虑将本文所提出的方法应用于艾滋病临床试验组研究175所得数据。该数据是R软件的内置数据集(ACTG175)，可直接在R中调用。这项随机临床试验比较了不同药物治疗对感染人类免疫缺陷病毒I型成人的治疗效果[16]。本文的研究目标是对获取到的患者数据进行亚组分析，以便在(20 ± 5)周时得到更精确的每个组的CD4计数(细胞/mm³)预测值。现结合实际病例与本文的研究目标，考虑选择以下9个协变量进行分析，具体见下表1。

Table 1. Various covariates and their meanings

表1. 各协变量及其含义

协变量	含义	单位
$X_1$	血友病(0 = 否，1 = 是)	—
$X_2$	性别(0 = 女性，1 = 男性)	—
$X_3$	基线CD4计数	—
$X_4$	直到首次出现以下情况的天数： 1) CD4T细胞计数至少下降50； 2) 表明进展为艾滋病的事件；3) 死亡	天
$X_5$	年龄	岁
$X_6$	重量	kg
$X_7$	Karnofsky分数	—
$X_8$	基线CD8计数	—
$X_9$	之前接受抗逆转录病毒治疗的天数	天

本文首先将$X_i={\left(1,X_{i1},\cdots ,X_{i9}\right)}^{\text{T}}$ 作为自变量，$Y=20\pm 5$ 周时的CD4计数作为因变量，结合相应数据拟合得到线性回归模型，即认为其没有亚组，并用${\widehat{\beta }}^{ols}$ 表示对数据进行普通的最小二乘估计得到的各个参数的估计值。选择$Z_i={\left(X_{i1},X_{i2},X_{i3},X_{i4},X_{i5}\right)}^{\text{T}}$ 作为阈值变量。

现考虑在艾滋病临床试验组研究所得数据中随机选取500名患者对应的变量数据，将其应用于本文所提出的基于秩回归估计的多阈值变平面模型中，可以得到系数估计值$\widehat{\beta }$ 和$\widehat{\delta }$ ，以及它们的标准误差(SE)，具体见下表2。

由表2分别可以得到两个亚组的各种解释变量的推论。例如，第一组基线CD4计数变量$X_3$ 的系数为0.452，前两组间的系数差为0.011，说明第二亚组的系数为0.452 + 0.011 = 0.463，第二和第三亚组之间的系数差为0.202，即说明第三亚组的系数为0.463 + 0.202 = 0.665。然而，在没有进行的亚组分析的简单线性回归模型中，对于基线CD4计数变量，普通最小二乘法只返回一个恒定的系数0.538，没有考虑到亚组内的基线CD4计数差异。其他系数也可以类似地解释。

Table 2. Estimated values of various coefficients and corresponding standard errors

表2. 各系数估计值以及相应的标准误差

	$\widehat{\beta }$		${\widehat{\delta }}_1$		${\widehat{\delta }}_2$		${\widehat{\beta }}^{ols}$
	系数	标准误差	系数	标准误差	系数	标准误差	系数	标准误差
$X_0$	−0.073	−0.054	0.123	0.105	0.241	0.082	−0.003	0.126
$X_1$	−0.165	0.015	−0.002	0.025	−0.011	0.104	−0.086	0.024
$X_2$	0.001	0.020	0	—	0	—	0.023	0.066
$X_3$	0.452	0.052	0.011	0.035	0.202	0.032	0.538	0.035
$X_4$	0.222	0.048	0	—	0.053	0.024	0.244	0.043
$X_5$	0	—	−0.036	0.045	0	—	−0.025	0.014
$X_6$	−0.006	0.014	0	—	0	—	−0.006	0.031
$X_7$	0.075	0.031	0.018	0.021	0.163	0.012	0.032	0.019
$X_8$	0.079	0.019	−0.023	0.029	−0.105	0.084	−0.054	0.034
$X_9$	−0.175	0.042	−0.234	0.037	−0.047	0.033	−0.087	0.051

现考虑将选取的艾滋病临床试验数据以7:3的比例分为训练集和测试集两部分，选择前70%作为训练集，后30%作为测试集。将得到的训练集对应的数据带入基于秩回归的多阈值变平面模型中，通过使用本文提出的两阶段计算方法，可以得到各参数估计结果，再将各参数估计结果带入测试集中，结合测试集中自变量对应的数据，即可获得相应的因变量，即CD4计数预测值。类似地，可以得到基于普通最小二乘估计的多阈值变平面模型的各参数估计结果，将其带入测试集中也可以得到预测的CD4计数。

根据得到的两种不同方法的CD4计数的预测值，借助R软件，可以画出它们预测的CD4计数对比的散点图，具体见下图1。

Figure 1. Comparison of prediction results by different methods scatter plot

图1. 不同方法预测结果对比散点图

由上图1得到的不同方法预测结果的散点图对比，可以发现秩回归估计方法和最小二乘估计方法的预测值都与实际观测值有微小差距，因此考虑绘制实际观测的CD4计数的数据分布直方图和Q-Q图，进一步说明本文选取的秩回归估计方法的合理性，具体见下图2。

(a) (b)

Figure 2. Histogram of CD4 count data distribution (left) and Q-Q chart (right)

图2. CD4计数的数据分布直方图(左)和Q-Q图(右)

根据图2展示的CD4计数数据分布直方图与Q-Q图可知，该数据并未完全遵循正态分布，特别是在分布的两端尾部，部分数据点显著偏离了正态分布的参照线，这表明数据集中存在异常值或离群数据，例如第190个与第464个数据点。因此，相较于一般的适用于独立同分布数据的最小二乘估计法，本文所提出的基于秩回归的多阈值变平面模型估计方法在此类数据上表现出更高的适用性。该方法展现出良好的鲁棒性，其估计效果更为优越。

此外，考虑将本文所提出的基于秩回归估计的多阈值变平面模型(MCPL)的预测性能与基于秩回归估计的单阈值变化平面(SCPL)模型、相关文献[7]提出的单阈值变点模型(MCPT)以及等加权多变量$Z_i=\left(X_1+X_2+X_3+X_4+X_5\right)/5$ 的变平面模型(E-MCPL)的预测性能进行了比较。由于本例选取的数据中，$X_1$ 、$X_2$ 为不连续变量，所以不能应用于单阈值变点模型。本文总结了上述提出的所有方法的预测误差，具体见下表3。

Table 3. Prediction errors of model estimation results using different methods

表3. 不同方法的模型估计结果的预测误差

方法	预测误差	$\widehat{s}$	阈值$\widehat{a}$	分组结果
MCPL	0.565	2	(−0.436, 0.315)	390:66:44
SCPL	0.581	1	0.213	105:395
MCPT-X3	0.578	1	−0.102	198:302
MCPT-X4	0.585	0	—	—
MCPT-X5	0.575	2	(−1.312, 1.513)	148:305:47
E-MCPL	0.587	0	—	—
OLS	0.597	0	—	—

注：MCPL表示本文提出的基于秩回归的多阈值变平面模型估计方法；SCPL表示基于秩回归的单阈值变平面模型估计方法；MCPT-X₃表示阈值变量为X₃的单阈值变点模型估计方法；MCPT-X₄表示阈值变量为X₄的单阈值变点模型估计方法；MCPT-X₅表示阈值变量为X₅的单阈值变点模型估计方法；E-MCPL表示等加权多变量$Z_i=\left(X_1+X_2+X_3+X_4+X_5\right)/5$ 的变平面模型估计方法；OLS表示普通的最小二乘模型估计方法。

由上表3结果可知，对比其他几种方法的估计结果，本文提出的基于秩回归的多阈值变平面模型估计方法的预测误差是最小的。即表明本文所提出的方法估计更精确、更有效。此外，为了研究亚组，本文总结了几种模型估计方法的分组结果，见上述表3的最后一列。现考虑通过所有协变量对应的不同亚组的平均值来展示其分组结果，具体见下图3。

(a) (b)

(c) (d)

Figure 3. Average values of covariates corresponding to subgroups for different model estimation methods

图3. 不同模型估计方法的协变量对应亚组的平均值

由上述图3中不同模型估计方法得到的每个协变量分组后组内平均值的条形分布图可知，本文提出的基于秩回归的多阈值变平面模型的分组效果更显著，能够更好地根据病人的不同特征(即自变量)将其分到相应的组别，采取适合该病人的治疗方法，即更好地体现了个性化医疗。

5. 结论

本文的创新点是，相对于普通的最小二乘估计法，本文所提出的基于秩回归的多阈值变平面模型估计方法具有很好的鲁棒性，而且此方法还添加了惩罚函数，可以扩展到高维的数据分析。此外，当特征空间的维数超高且为样本大小的指数阶时，可以通过基于秩回归的多阈值变平面模型估计方法来估计单个或多个标记的影响，进而考虑单个或多个阈值的变化，以便得到更精确的结果。

本文得出的结论可以应用于医学中的个性化医疗技术，以便为更多病情复杂的患者找到适合他们的治疗方案，从而更好、更有效地对其进行精准治疗，帮助患者早日康复，不再忍受病痛的折磨。