基于成对融合惩罚的左删失数据子群分析

Subgroup analysis for left-censored data based on pairwise fusion penalty

摘要

基于Tobit回归模型，使用成对融合惩罚的正则化方法，对具有异质性的左删失数据进行子群分析，实现了回归参数估计与子群识别的同步进行。通过引入一组新变量，将原优化问题转化为可以用交替方向乘子法求解的仅含等式约束的多变量优化问题。并且，将每一步迭代的目标函数中与损失相关的多变量函数，利用广义坐标下降算法转化为一组二次优化单变量函数。证明所提算法的收敛性，并建立所得参数估计量的大样本性质。模拟研究和实际数据分析表明所提方法具有良好性能。

Abstract

We use pairwise fusion penalty regularization method， based on Tobit regression model， to perform subgroup analysis on left-censored data with heterogeneity， simultaneously estimating regression parameters and identifying subgroups. By introducing a set of new parameters， the original optimization problem is transformed into a multivariate optimization problem with equality constraints only that can be solved by alternating direction method of multipliers. Moreover， the multivariate function related to the loss in each iteration is transformed into a group of quadratic surrogate functions of single variable by generalized coordinate descent algorithm. We prove that the proposed algorithm is convergent， and establish the large sample properties of the obtained parameter estimators. Simulation studies and real data analysis show that the proposed method has good performance.

关键词

左删失数据 / Tobit模型 / 子群识别 / 成对融合惩罚

Key words

left-censored data / Tobit model / subgroup identification / pairwise fusion penalty

删失数据是经济学、医学、农业、生物和工程等领域研究中经常需要处理的一种数据类型，其观测值仅部分已知。Tobin^［1］对删失因子进行多元线性回归分析，得到Tobit模型。对 Tobit模型的现有研究主要集中于放宽模型假设和发展不同参数估计方法两方面。Tobit模型依赖于误差项的正态假设。有学者放宽此分布假设，仅使用误差项的均值和方差^［2-3］。Miller^［3］使用未删失观测值的加权线性组合作为斜率和截距项的估计量。Wei等^［4］指出，不指定误差项的分布阻碍了有现实用处和理论支持的统计方法的发展。针对此问题做出贡献的有文献［2-3，5-6］等。尤其是Prentice^［5］的工作，在指定不同的误差项分布下，基于广义秩向量的边际概率提出一类线性秩统计量。此类统计量被证明是渐近正态的，并且无需假设误差项分布^［7］。秩估计量也得到进一步研究，并发展出R包^［4，8］。在发展Tobit模型参数估计方法方面，有最小二乘估计^［9-11］和基于经验似然方法的研究^［12-13］。然而，这些模型考虑的删失数据并不存在异质性。

在精准医学、经济学和社会学中，人们需要区分异质性个体，并对属于不同子群个体的响应变量和预测因子之间的关系进行建模。一个关键问题是如何刻画个体之间的差异。现有研究对异质性的刻画主要从分布的参数和解释变量着手。在分布的参数方面，有限混合效应模型^［14］将所有样本视为来自多个不同总体之混合，不同总体用不同的均值、协方差阵区别。这样的模型框架简单，得到许多发展和扩充^［15-16］。Logistic-Cox模型^［17］把有限混合模型推广到删失数据。然而，混合效应模型要求事先指定子群数目以及每个子群的潜在分布类型。Ma和Huang^［18］使用线性回归模型中个体各自的截距项解释异质性。这样的模型设定相对简单，无需分布假设和指定子群数目，且需要估计的参数数量相对混合效应模型减少许多。还有研究允许异质性体现在解释变量之中^［19-20］，用多个协变量的线性组合解释异质性，模型中不同个体的协变量有各自的系数向量。以上文献展示了对潜在异质性的刻画。然而，这些研究所考虑的异质性数据鲜有删失背景。

正则化方法已被推广应用到子群识别问题中。Ma和Huang^［18］使用一些具有优良性质（如无偏性）的凹惩罚函数，对成对截距项之差进行压缩。然而，此研究针对连续数据使用的最小二乘损失无法直接应用于删失数据。关于删失数据的子群分析，文献［21］基于加速失效模型，使用平均插补法处理右删失数据之后，对最小二乘损失压缩成对截距项之差来识别子群。

本文基于Tobit模型，对异质左删失数据提出一种新方法，结合成对融合凹惩罚进行子群分析。该方法不需要异质性子群数目的先验信息，而是通过压缩不同个体之间的截距之差，就可以在估计参数的同时，自动实现子群识别（相关理论结果的证明请见期刊官网http：//journal.ucas.ac.cn在线补充文档）。

1 模型与方法

1.1 模型设定

假设对个体i观测到一组协变量{\mathit{x}}_i=(x_{i\mathrm{1}},\dots ,x_{ip}{)}^{\mathrm{T}}以及一个观测响应变量y_i，i=\mathrm{1},\dots ,n，记\mathit{X}=({\mathit{x}}_{\mathrm{1}},\dots ,{\mathit{x}}_n{)}^{\mathrm{T}}。在Tobit回归中，由于响应变量的观测删失特点，假设存在一个潜在响应变量y_i^{\mathrm{*}}，使得y_i=\mathrm{m}\mathrm{a}\mathrm{x}\{y_i^{\mathrm{*}},c\}。其中，c是已知常数。为方便起见，在接下来的研究中假设c=\mathrm{0}。

考虑来自异质总体的观测响应变量\mathit{y}=(y_{\mathrm{1}},\cdots ,y_n{)}^{\mathrm{T}}，假设不同子群体在基准水平上不同，即使用个体具有各自截距的回归模型进行建模。因此，考虑如下Tobit线性模型

y_i^{\mathrm{*}}={\mu }_i+{\mathit{x}}_i^{\mathrm{T}}\mathit{\beta }+{ϵ}_i,i=\mathrm{1},\cdots ,n,（1）

其中， \mathit{\mu }=({\mu }_{\mathrm{1}},\cdots ,{\mu }_n{)}^{\mathrm{T}}是异质个体的截距， \mathit{\beta }=({\beta }_{\mathrm{1}},\cdots ,{\beta }_p{)}^{\mathrm{T}}是公共的系数向量，独立同分布的误差项\mathit{ϵ}=({ϵ}_{\mathrm{1}},\cdots ,{ϵ}_n{)}^{\mathrm{T}}满足{ϵ}_i~N(\mathrm{0},{\sigma }^{\mathrm{2}})。易得P(y_i^{\mathrm{*}}\le \mathrm{0})=\Phi ((-{\mu }_i-{\mathit{x}}_i^{\mathrm{T}}\mathit{\beta })/\sigma )，其中， \Phi (·)为标准正态分布函数。记d_i=I_{y_i>\mathrm{0}}，则删失数据y_i的Tobit似然如下

L_n\left(\mathit{\mu },\mathit{\beta },\sigma \right)={{\displaystyle \prod _{i=\mathrm{1}}^n}\left[\frac{\mathrm{1}}{\sqrt[]{\mathrm{2}\mathrm{\pi }}\sigma }\mathrm{e}\mathrm{x}\mathrm{p}\left\{-\frac{\mathrm{1}}{\mathrm{2}{\sigma }^{\mathrm{2}}}(y_i-{\mu }_i-{\mathit{x}}_i^{\mathrm{T}}{\mathit{\beta })}^{\mathrm{2}}\right\}\right]}^{d_i}

{\left[\Phi ((-{\mu }_i-{\mathit{x}}_i^{\mathrm{T}}\mathit{\beta })/\sigma )\right]}^{\mathrm{1}-d_i}.（2）

假定n个个体来自K\ge \mathrm{1}个不同的组\mathrm{𝒢}=\left\{{\mathrm{𝒢}}_{\mathrm{1}},\cdots ,{\mathrm{𝒢}}_K\right\}，其中，{\mathrm{𝒢}}_i表示第i个组的个体组(i=\mathrm{1},\cdots ,K)。在实际中，K是未知的，假定K\ll n是合理的。对同属一组的个体i\in {\mathrm{𝒢}}_k，有相同的未知潜在截距{\tau }_k，\mathit{\tau }=({\tau }_{\mathrm{1}},\cdots ,{\tau }_K{)}^{\mathrm{T}}。

Amemiya^［9］证明了式（2）的似然函数不是凹函数。然而， Olsen^［10］提出重新参数化的方法，即令\mathit{\alpha }=\mathit{\mu }/\sigma，\mathit{\delta }=\mathit{\beta }/\sigma，\gamma =\mathrm{1}/\sigma，使得

\begin{array}{l}\mathrm{l}\mathrm{n}{L}_n(\mathit{\alpha },\mathit{\delta },\gamma )={\sum }_{i=\mathrm{1}}^n{d}_i\left[\mathrm{l}\mathrm{n}\gamma -\frac{\mathrm{1}}{\mathrm{2}}(\gamma y_i-{\alpha }_i-{\mathit{x}}_i^{\mathrm{T}}{\mathit{\delta })}^{\mathrm{2}}\right]+\\ (\mathrm{1}-d_i)\mathrm{l}\mathrm{n}\left[\Phi (-{\alpha }_i-{\mathit{x}}_i^{\mathrm{T}}\mathit{\delta })\right]\end{array}

为凹函数；同时，他还证明了重新参数化后Tobit对数似然的全局凹性。

Ma和Huang^［18］对子群分析提出成对融合惩罚方法。该方法在凸损失之后加上对成对参数之差的凹惩罚函数，可以在不知道任何分组先验信息的情况下，达到自动分组的效果。据此，本节考虑以下目标函数

（3）

其中， {\mathit{𝓁}}_n(\mathit{\alpha },\mathit{\delta },\gamma )=-\mathrm{l}\mathrm{n}{L}_n(\mathit{\alpha },\mathit{\delta },\gamma )/n是凸损失函数，凹惩罚函数p(·,\lambda )的调节参数\lambda \ge \mathrm{0}。

惩罚函数的选择很重要，不合适的惩罚函数不能有效恢复子群结构。Ma和Huang^［18］指出， L₁惩罚p_{}(t,\lambda )=\lambda t存在识别不出子群效应或者过度拟合子群结构的问题。因此，本节倾向于使用具有良好性质（包括无偏性、稀疏性）的惩罚函数，比如平滑剪切绝对偏差惩罚函数（smoothly clipped absolute deviation，SCAD）^［22］和最小曲率惩罚函数（minimax concave penalty，MCP）^［23］。SCAD具有以下形式

p_a(t,\lambda )=\lambda {\int }_{\mathrm{0}}^{t}I(x\le \lambda )+\frac{{(a\lambda -x)}_{+}}{(a-\mathrm{1})\lambda }I(x>\lambda )\mathrm{d}x,a>\mathrm{2},

MCP具有以下形式

p_a(t,\lambda )={\int }_{\mathrm{0}}^t\frac{{(a\lambda -x)}_{+}}{a}\mathrm{d}x,a>\mathrm{1},

其中，参数a控制惩罚函数的凹性。根据Fan和Li^［22］及Zhang^［23］的研究，本文把a当作固定值，对SCAD惩罚取a=\mathrm{3.7}，对MCP惩罚取a=\mathrm{3}。

1.2 参数估计算法

通过引入一列新参数{\eta }_{ij}={\alpha }_i-{\alpha }_j，，最小化Q_n(\mathit{\alpha },\mathit{\delta },\gamma )的问题转化为如下问题

\mathrm{s}.\mathrm{t}.\hspace{1em}{\alpha }_i-{\alpha }_j-{\eta }_{ij}=\mathrm{0}.（4）

在交替方向乘子法（alternating direction multiplier method， ADMM）下，上述只含等式约束的优化问题等价于最小化如下增广拉格朗日函数

（5）

其中，对偶变量\mathit{\eta }是拉格朗日乘子，\vartheta是惩罚参数。利用ADMM的框架，优化式（4）的算法过程如下：

对k=\mathrm{0,1},\dots，重复如下步骤直至收敛

\begin{array}{l}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{p}\mathrm{1}:({\mathit{\alpha }}^{(k+\mathrm{1})},{\mathit{\delta }}^{(k+\mathrm{1})},{\gamma }^{(k+\mathrm{1})})=\\ \mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}{\mathrm{n}}_{\mathit{\alpha },\mathit{\delta },\gamma }\mathrm{ℒ}(\mathit{\alpha },\mathit{\delta },\gamma ,{\mathit{\eta }}^{(k)},{\mathit{\nu }}^{(k)}),\end{array}

\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{p}\mathrm{2}:{\mathit{\eta }}^{(k+\mathrm{1})}=\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}{\mathrm{n}}_{\mathit{\eta }}\mathrm{ℒ}({\mathit{\alpha }}^{(k+\mathrm{1})},{\mathit{\delta }}^{(k+\mathrm{1})},{\gamma }^{(k+\mathrm{1})},\mathit{\eta },{\mathit{\nu }}^{(k)}),

\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{p}\mathrm{3}:{\nu }_{ij}^{(k+\mathrm{1})}={\nu }_{ij}^{(k)}+\vartheta ({\alpha }_i^{(k+\mathrm{1})}-{\alpha }_j^{(k+\mathrm{1})}-{\eta }_{ij}^{(k+\mathrm{1})}).（6）

step 2中，当给定{\eta }_{ij}之外的其他参数时， \mathrm{ℒ}(\mathit{\alpha },\mathit{\delta },\gamma ,\mathit{\eta },\mathit{\nu })作为{\eta }_{ij}的单变量函数是凹函数（要求SCAD下a>{\vartheta }^{-\mathrm{1}}+\mathrm{1}，MCP下a>{\vartheta }^{-\mathrm{1}}）。此时， step 2等价于

{\eta }_{ij}^{(k+\mathrm{1})}=\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}{\mathrm{n}}_{{\eta }_{ij}}p(|{\eta }_{ij}|,\lambda )+\frac{\vartheta }{\mathrm{2}}({\eta }_{ij}-c_{ij}{)}^{\mathrm{2}},

其中， c_{ij}={\alpha }_i^{(k+\mathrm{1})}-{\alpha }_j^{(k+\mathrm{1})}+{\vartheta }^{-\mathrm{1}}{\nu }_{ij}^{(k)}。L₁惩罚下

{\widehat{\eta }}_{ij}=\mathrm{S}\mathrm{T}(c_{ij},\lambda /\vartheta ),（7）

其中，软阈值函数\mathrm{S}\mathrm{T}\left(t,\lambda \right)=\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left(t\right)\mathrm{m}\mathrm{a}\mathrm{x}\left\{|t|-\lambda ,\mathrm{0}\right\}。MCP惩罚下

{\widehat{\eta }}_{ij}=\left\{\begin{array}{ll}\frac{\mathrm{S}\mathrm{T}(c_{ij},\lambda /\vartheta )}{\mathrm{1}-{(a\vartheta )}^{-\mathrm{1}}},& |c_{ij}| \le a\lambda ,\\ c_{ij},& |c_{ij}| >a\lambda .\end{array}\right.（8）

SCAD惩罚下

（9）

step 1中，形式复杂的Tobit似然函数使得无法得到参数迭代解析式。因此，需要使用Yang和Zou^［24］提出的广义坐标下降法（generalized coordinate descent， GCD），GCD算法收敛后即得到({\mathit{\alpha }}^{(k+\mathrm{1})},{\mathit{\delta }}^{(k+\mathrm{1})},{\gamma }^{(k+\mathrm{1})})。GCD把满足二次优化条件的多变量优化问题转化为单变量优化问题，使得每一步坐标更新变得简单并且高效。

定义1.1 对于一个单变量函数f:\mathbb{R}\to \mathbb{R}，若存在M\in {\mathbb{R}}^{+}，使得对任意t,a\in \mathbb{R}，有

f(t+a)\le f(t)+f^{\text{'}}(t)a+\frac{M}{\mathrm{2}}{a}^{\mathrm{2}},

则称函数f满足二次优化条件。

值得注意的是，重新参数化后得到的Tobit损失满足二次优化条件。

引理1.1 当固定其他参数值，将{\mathit{𝓁}}_n看作{\alpha }_i的单变量函数时， {\mathit{𝓁}}_n满足二次优化条件，M_i=\mathrm{1}；当固定其他参数值，将{\mathit{𝓁}}_n看作{\delta }_j的单变量函数时， {\mathit{𝓁}}_n满足二次优化条件，M_j={\sum }_{i=\mathrm{1}}^n{x}_{ij}^{\mathrm{2}}/n，对协变量进行标准化处理后，M_j=\mathrm{1}。

为方便起见，本节假设协变量已进行标准化处理。

首先基于当前值(\tilde{\mathit{\alpha }},\tilde{\mathit{\delta }},\tilde{\gamma },\tilde{\mathit{\eta }},\tilde{\mathit{\nu }})更新{\mathit{\alpha }}^{(k,t+\mathrm{1})}。注意到，Tobit损失只有第i项和{\alpha }_i有关。固定其他参数，得到关于{\alpha }_i的单变量函数{\mathit{𝓁}}_n({\alpha }_i|\tilde{\mathit{\alpha }},\tilde{\mathit{\delta }},\tilde{\gamma })=\frac{\mathrm{1}}{n}\left[\frac{\mathrm{1}}{\mathrm{2}}{d}_i(\tilde{\gamma }{y}_i-{\alpha }_i-\right.\left.{\mathit{x}}_i^{\mathrm{T}}{\tilde{\mathit{\delta }})}^{\mathrm{2}}-(\mathrm{1}-d_i)\mathrm{l}\mathrm{n}\Phi (-{\alpha }_i-{\mathit{x}}_i^{\mathrm{T}}\tilde{\mathit{\delta }})\right]。根据引理1.1，构造二次优化函数Q_{\alpha }({\alpha }_i|\tilde{\mathit{\alpha }},\tilde{\mathit{\delta }},\tilde{\gamma })={\mathit{𝓁}}_n({\tilde{\alpha }}_i|\tilde{\mathit{\alpha }},\tilde{\mathit{\delta }},\tilde{\gamma })+{\mathit{𝓁}}_n^{\text{'}}({\tilde{\alpha }}_i|\tilde{\mathit{\alpha }},\tilde{\mathit{\delta }},\tilde{\gamma })\left({\alpha }_i-{\tilde{\alpha }}_i\right)+{\left({\alpha }_i-{\tilde{\alpha }}_i\right)}^{\mathrm{2}}/\mathrm{2}。根据优化最小化准则， 用Q_{\alpha }({\alpha }_i|\tilde{\mathit{\alpha }},\tilde{\mathit{\delta }},\tilde{\gamma })替代{\mathit{𝓁}}_n({\alpha }_i|\tilde{\mathit{\alpha }},\tilde{\mathit{\delta }},\tilde{\gamma })， 故{\alpha }_i^{(k,t+\mathrm{1})}可以通过最小化Q_{\alpha }({\alpha }_i|\tilde{\mathit{\alpha }},\tilde{\mathit{\delta }},\tilde{\gamma })+{\tilde{\mathit{\nu }}}^{\mathrm{T}}(\mathit{\Delta }\mathit{\alpha }-\tilde{\mathit{\eta }})+\frac{\vartheta }{\mathrm{2}}{‖\mathit{\Delta }\mathit{\alpha }-\tilde{\mathit{\eta }}‖}_{\mathrm{2}}^{\mathrm{2}}得到。为使用全部样本信息，将上式对i求和，有

{\mathit{\alpha }}^{(k,t+\mathrm{1})}=({\mathit{I}}_n+n\vartheta {\mathit{\Delta }}^{\mathrm{T}}{\mathit{\Delta })}^{-\mathrm{1}}(\tilde{\mathit{\alpha }}-{\tilde{\mathit{h}}}_{\alpha }+n\vartheta {\mathit{\Delta }}^{\mathrm{T}}\tilde{\mathit{\eta }}-n{\mathit{\Delta }}^{\mathrm{T}}\tilde{\mathit{\nu }}),（10）

其中， {\tilde{\mathit{h}}}_{\alpha }=\{{\mathit{𝓁}}_n^{\text{'}}({\tilde{\alpha }}_i|\tilde{\mathit{\alpha }},\tilde{\mathit{\delta }},\tilde{\gamma }),i=\mathrm{1},\dots ,n\}\in {\mathbb{R}}^n。此时，\tilde{\mathit{\alpha }}={\mathit{\alpha }}^{(k,t)},\tilde{\mathit{\delta }}={\mathit{\delta }}^{(k,t)},\tilde{\gamma }={\gamma }^{(k,t)},\tilde{\mathit{\eta }}={\mathit{\eta }}^{(k)},\tilde{\mathit{\nu }}={\mathit{\nu }}^{(k)}。

接下来基于当前值(\tilde{\mathit{\alpha }},\tilde{\mathit{\delta }},\tilde{\gamma },\tilde{\mathit{\eta }},\tilde{\mathit{\nu }})更新{\delta }_j^{(k,t+\mathrm{1})}。固定其他参数，将Tobit损失看作{\delta }_j的单变量函数，有

\begin{array}{l}{\mathit{𝓁}}_n({\delta }_j|\tilde{\mathit{\alpha }},\tilde{\mathit{\delta }},\tilde{\gamma })=\frac{\mathrm{1}}{n}{\sum }_{i=\mathrm{1}}^n\frac{\mathrm{1}}{\mathrm{2}}{d}_i(\tilde{\gamma }{y}_i-{\tilde{\alpha }}_i-x_{ij}{\delta }_j-{\mathit{x}}_{i,(-j)}^{\mathrm{T}}{\tilde{\mathit{\delta }}}_{(-j)}{)}^{\mathrm{2}}-\\ (\mathrm{1}-d_i)\mathrm{l}\mathrm{n}\Phi (-{\tilde{\alpha }}_i-x_{ij}{\delta }_j-{\mathit{x}}_{i,(-j)}^{\mathrm{T}}{\tilde{\mathit{\delta }}}_{(-j)}).\end{array}

根据引理1.1和优化最小化准则，用二次优化函数{\mathit{𝓁}}_n({\tilde{\delta }}_j|\tilde{\mathit{\alpha }},\tilde{\mathit{\delta }},\tilde{\gamma })+{\mathit{𝓁}}_n^{\text{'}}({\tilde{\delta }}_j|\tilde{\mathit{\alpha }},\tilde{\mathit{\delta }},\tilde{\gamma })({\delta }_j-{\tilde{\delta }}_j)+({\delta }_j-{\tilde{\delta }}_j{)}^{\mathrm{2}}/\mathrm{2}替代{\mathit{𝓁}}_n({\delta }_j|\tilde{\mathit{\alpha }},\tilde{\mathit{\delta }},\tilde{\gamma })，故

{\delta }_j^{(k,t+\mathrm{1})}={\tilde{\delta }}_j-{\tilde{h}}_{{\delta }_j},（11）

其中， {\tilde{h}}_{{\delta }_j}={\mathit{𝓁}}_n^{\text{'}}({\tilde{\delta }}_j|\tilde{\mathit{\alpha }},\tilde{\mathit{\delta }},\tilde{\gamma })\in \mathbb{R}。此时\tilde{\mathit{\alpha }}={\mathit{\alpha }}^{(k,t)},。

最后更新{\gamma }^{(k,t+\mathrm{1})}。给定({\mathit{\alpha }}^{(k,t+\mathrm{1})},{\mathit{\delta }}^{(k,t+\mathrm{1})})，{\gamma }^{(k,t+\mathrm{1})}=\begin{array}{l}\frac{\mathrm{1}}{\mathrm{2}{\sum }_{i=\mathrm{1}}^n{d}_i{y}_i^{\mathrm{2}}}{\displaystyle \sum _{i=\mathrm{1}}^n}({\alpha }_i^{(k,t+\mathrm{1})}+{\mathit{x}}_i^T{\mathit{\delta }}^{(k,t+\mathrm{1})})d_i{y}_i+\frac{\mathrm{1}}{\mathrm{2}{\sum }_{i=\mathrm{1}}^n{d}_i{y}_i^{\mathrm{2}}}\\ \sqrt[]{{\left(({\alpha }_i^{(k,t+\mathrm{1})}+{\mathit{x}}_i^{\mathrm{T}}{\mathit{\delta }}^{(k,t+\mathrm{1})})d_i{y}_i\right)}^{\mathrm{2}}+\mathrm{4}{\sum }_{i=\mathrm{1}}^n{d}_i{\sum }_{i=\mathrm{1}}^n{d}_i{y}_i^{\mathrm{2}}}.\end{array}（12）

下面给出ADMM-GCD算法的流程。

下面的命题给出ADMM-GCD算法的收敛性质。

命题1.1 当a>{\vartheta }^{-\mathrm{1}}+\mathrm{1}（SCAD）或a>{\vartheta }^{-\mathrm{1}}（MCP）时， ADMM-GCD算法生成的任意聚点({\mathit{\alpha }}^{(m+\mathrm{1})},{\mathit{\delta }}^{(m+\mathrm{1})},{\gamma }^{(m+\mathrm{1})},{\mathit{\eta }}^{(m+\mathrm{1})})都是\mathrm{ℒ}(\mathit{\alpha },\mathit{\delta },\gamma ,\mathit{\eta },{\mathit{\nu }}^{(m)})沿坐标方向的最小值，并且原始残差{\mathit{r}}^{(m)}=\mathit{\Delta }{\mathit{\alpha }}^{(m)}-{\mathit{\eta }}^{(m)}和对偶残差{\mathit{s}}^{(m)}=\vartheta {\mathit{\Delta }}^{\mathrm{T}}({\mathit{\eta }}^{(m+\mathrm{1})}-{\mathit{\eta }}^{(m)})满足\underset{m\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}{‖{\mathit{r}}^{(m)}‖}_{\mathrm{2}}^{\mathrm{2}}=\mathrm{0} 以及\underset{m\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}{‖{\mathit{s}}^{(m)}‖}_{\mathrm{2}}^{\mathrm{2}}=\mathrm{0}。

2 理论性质

定义{\mathrm{ℳ}}_{\mathrm{𝒢}}=\{\mathit{\alpha }\in {\mathbb{R}}^n|{\alpha }_i={\alpha }_j,\forall i,j\in {\mathrm{𝒢}}_k,\mathrm{1}\le k\le K\}。任意\mathit{\alpha }\in {\mathrm{ℳ}}_{\mathrm{𝒢}}，\mathit{\alpha }=\mathit{Z}\mathit{\tau }，其中，\mathit{Z}=\{z_{ik}\}是n\times K矩阵，当i\in {\mathrm{𝒢}}_k时，z_{ik}=\mathrm{1}，否则z_{ik}=\mathrm{0}。矩阵运算有{\mathit{Z}}^{\mathrm{T}}\mathit{Z}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(|{\mathrm{𝒢}}_{\mathrm{1}}|,\dots ,|{\mathrm{𝒢}}_K|)，其中，|{\mathrm{𝒢}}_k|表示{\mathrm{𝒢}}_k中元素的个数。

设{\mathit{\alpha }}^{\mathrm{*}},{\mathit{\delta }}^{\mathrm{*}},{\gamma }^{\mathrm{*}},和{\mathit{\tau }}^{\mathrm{*}}分别是参数\mathit{\alpha },\mathit{\delta },\gamma和\mathit{\tau }的真实值。定义\mathit{\theta }=({\theta }_{\mathrm{1}},\dots ,{\theta }_{K+p},{\theta }_{K+p+\mathrm{1}}{)}^{\mathrm{T}}:=({\mathit{\tau }}^{\mathrm{T}},{\mathit{\delta }}^{\mathrm{T}}{,\gamma )}^{\mathrm{T}}，对应真实值为{\mathit{\theta }}^{\mathrm{*}}=(({\mathit{\tau }}^{\mathrm{*}}{)}^{\mathrm{T}},({\mathit{\delta }}^{\mathrm{*}}{)}^{\mathrm{T}},{\gamma }^{\mathrm{*}}{)}^{\mathrm{T}}。

Tobit似然的形式展示出未删失观测值和删失观测值对应似然的差异。假设有n_{\mathrm{1}}个y_i>\mathrm{0}，对应响应变量为{\mathit{y}}_{\mathrm{1}}，有n_{\mathrm{0}}=n-n_{\mathrm{1}}个y_i\le \mathrm{0}，对应响应变量为{\mathit{y}}_{\mathrm{0}}。重新排列观察结果

\mathit{X}=\left[\begin{array}{c}{\mathit{X}}_{\mathrm{0}}\\ {\mathit{X}}_{\mathrm{1}}\end{array}\right],\mathit{y}=\left[\begin{array}{c}{\mathit{y}}_{\mathrm{0}}\\ {\mathit{y}}_{\mathrm{1}}\end{array}\right],\mathit{Z}=\left[\begin{array}{c}{\mathit{Z}}_{\mathrm{0}}\\ {\mathit{Z}}_{\mathrm{1}}\end{array}\right],

进而有

\mathit{ϵ}=\left[\begin{array}{c}{\mathit{ϵ}}_{\mathrm{0}}\\ {\mathit{ϵ}}_{\mathrm{1}}\end{array}\right],\mathit{\alpha }=\left[\begin{array}{c}{\mathit{\alpha }}_{\mathrm{0}}\\ {\mathit{\alpha }}_{\mathrm{1}}\end{array}\right],{\mathit{\alpha }}^{\mathrm{*}}=\left[\begin{array}{c}{\mathit{\alpha }}_{\mathrm{0}}^{\mathrm{*}}\\ {\mathit{\alpha }}_{\mathrm{1}}^{\mathrm{*}}\end{array}\right].

当真实组结构\mathrm{𝒢}已知时，\mathit{\alpha },\mathit{\delta }和\gamma的Oracle估计为

\begin{array}{l}({\widehat{\mathit{\alpha }}}^{\mathrm{o}\mathrm{r}},{\widehat{\mathit{\delta }}}^{\mathrm{o}\mathrm{r}},{\widehat{\gamma }}^{\mathrm{o}\mathrm{r}})=\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}{\mathrm{n}}_{\mathit{\alpha }\in {\mathrm{ℳ}}_{\mathrm{𝒢}},\mathit{\delta },\gamma }{\mathit{𝓁}}_n(\mathit{\alpha },\mathit{\delta },\gamma )\\ =\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}{\mathrm{x}}_{\mathit{\alpha }\in {\mathrm{ℳ}}_{\mathrm{𝒢}},\mathit{\delta },\gamma }\mathrm{l}\mathrm{n}{L}_n(\mathit{\alpha },\mathit{\delta },\gamma ).\end{array}（13）

相应地，

\begin{array}{l}{\widehat{\mathit{\theta }}}^{\mathrm{o}\mathrm{r}}=(({\widehat{\mathit{\tau }}}^{\mathrm{o}\mathrm{r}}{)}^{\mathrm{T}},({\widehat{\mathit{\delta }}}^{\mathrm{o}\mathrm{r}}{)}^{\mathrm{T}},{\widehat{\gamma }}^{\mathrm{o}\mathrm{r}}{)}^{\mathrm{T}}=\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}{\mathrm{n}}_{\mathit{\tau },\mathit{\delta },\gamma }{\mathit{𝓁}}_n(\mathit{Z}\mathit{\tau },\mathit{\delta },\gamma )\\ =\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}{\mathrm{x}}_{\mathit{\tau },\mathit{\delta },\gamma }\mathrm{l}\mathrm{n}{L}_n(\mathit{Z}\mathit{\tau },\mathit{\delta },\gamma ).\end{array}（14）

记\mathit{\zeta }=({\mathit{\alpha }}^{\mathrm{T}},{\mathit{\delta }}^{\mathrm{T}}{)}^{\mathrm{T}}，\tilde{\mathit{X}}=(\mathit{Z},\mathit{X})，仍然可以进行重新排列\tilde{\mathit{X}}=({\tilde{\mathit{X}}}_{\mathrm{0}}^{\mathrm{T}},{\tilde{\mathit{X}}}_{\mathrm{1}}^{\mathrm{T}}{)}^{\mathrm{T}}，其中，{\tilde{\mathit{X}}}_{\mathrm{0}}=({\tilde{\mathit{x}}}_{\mathrm{01}},\cdots ,{\tilde{\mathit{x}}}_{\mathrm{0}{n}_{\mathrm{0}}}{)}^{\mathrm{T}}。则

{\widehat{\mathit{\theta }}}^{\mathrm{o}\mathrm{r}}=(({\widehat{\mathit{\zeta }}}^{\mathrm{o}\mathrm{r}}{)}^{\mathrm{T}},{\widehat{\gamma }}^{\mathrm{o}\mathrm{r}}{)}^{\mathrm{T}}=\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}{\mathrm{x}}_{\mathit{\theta }}{S}_n(\mathit{\theta }),

其中，S_n(\mathit{\theta })={\sum }_{i=\mathrm{1}}^n{d}_i\left[\mathrm{l}\mathrm{n}\gamma -(\gamma y_i-{\tilde{\mathit{x}}}_i^{\mathrm{T}}{\mathit{\zeta })}^{\mathrm{2}}/\mathrm{2}\right]+(\mathrm{1}-d_i)\mathrm{l}\mathrm{n}\Phi (-{\tilde{\mathit{x}}}_i^{\mathrm{T}}\mathit{\zeta })。显然，上式对应重新参数化后的标准Tobit回归，其二阶偏导矩阵为

{\nabla }^{\mathrm{2}}{S}_n(\mathit{\theta })=\left[\begin{array}{cc}{\tilde{\mathit{X}}}_{\mathrm{0}}^{\mathrm{T}}D(\mathit{\zeta }){\tilde{\mathit{X}}}_{\mathrm{0}}& \mathrm{0}\\ \mathrm{0}& -n_{\mathrm{1}}{\gamma }^{-\mathrm{2}}\end{array}\right]+\left[\begin{array}{cc}-{\tilde{\mathit{X}}}_{\mathrm{1}}^{\mathrm{T}}{\tilde{\mathit{X}}}_{\mathrm{1}}& {\tilde{\mathit{X}}}_{\mathrm{1}}^{\mathrm{T}}{\mathit{y}}_{\mathrm{1}}\\ {\mathit{y}}_{\mathrm{1}}^{\mathrm{T}}{\tilde{\mathit{X}}}_{\mathrm{1}}& -{\mathit{y}}_{\mathrm{1}}^{\mathrm{T}}{\mathit{y}}_{\mathrm{1}}\end{array}\right],（15）

其中， D(\mathit{\zeta })=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(g^{\text{'}}(-{\tilde{\mathit{x}}}_{\mathrm{01}}^{\mathrm{T}}\mathit{\zeta }),\dots ,g^{\text{'}}(-{\tilde{\mathit{x}}}_{\mathrm{0}{n}_{\mathrm{0}}}^{\mathrm{T}}\mathit{\zeta }))，g(·)=\varphi (·)/\Phi (·)，\varphi (·)是标准正态分布的密度函数。

Amemiya^［11］说明了{\nabla }^{\mathrm{2}}{S}_n(\mathit{\theta })是负定阵，而{\widehat{\mathit{\theta }}}^{\mathrm{o}\mathrm{r}}满足\nabla S_n({\widehat{\mathit{\theta }}}^{\mathrm{o}\mathrm{r}})=\mathrm{0}，因此， {\widehat{\mathit{\theta }}}^{\mathrm{o}\mathrm{r}}就是标准Tobit回归的极大似然估计。

记函数\rho (t)=p(t,\lambda )/\lambda，\mathit{X}=({\mathit{X}}_{\mathrm{1}},\cdots ,{\mathit{X}}_p)，|{\mathrm{𝒢}}_{\mathrm{m}\mathrm{i}\mathrm{n}}|=\mathrm{m}\mathrm{i}{\mathrm{n}}_{\mathrm{1}\le k\le K}|{\mathrm{𝒢}}_k|。当K\ge \mathrm{2}时，用

b_n=\mathrm{m}\mathrm{i}{\mathrm{n}}_{i\in {\mathrm{𝒢}}_k,j\in {\mathrm{𝒢}}_{k^{\text{'}}},k\ne k^{\text{'}}}|{\alpha }_i^{\mathrm{*}}-{\alpha }_j^{\mathrm{*}}|=\mathrm{m}\mathrm{i}{\mathrm{n}}_{k\ne k^{\text{'}}}\left|{\tau }_k^{\mathrm{*}}-{\tau }_{k^{\text{'}}}^{\mathrm{*}}\right|

表示2个子群潜在公共截距的最小偏差。假设下列条件成立：

条件2.1　参数空间\{\mathit{\theta }=(\tilde{\mathit{\delta }},\gamma )|\tilde{\mathit{\delta }}\in {\mathbb{R}}^{K+p},\gamma >\mathrm{0}\}是紧集，包含真值{\mathit{\theta }}^{\mathrm{*}}的一个开邻域；

条件2.2　对任意\mathrm{1}\le j\le p，||{\mathit{X}}_j|{|}_{\mathrm{2}}=\sqrt[]{n}；对某些常数，||\mathit{X}|{|}_{\mathrm{\infty }}\le C_{\mathrm{1}}p，{\lambda }_{\mathrm{m}\mathrm{i}\mathrm{n}}({\tilde{\mathit{X}}}^{\mathrm{T}}\tilde{\mathit{X}})\ge C_{\mathrm{2}}\left|{\mathrm{𝒢}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\right|；

条件2.3　惩罚函数p(t,\lambda )关于t对称，\rho (\mathrm{0})=\mathrm{0}，且存在使得\rho (t)在[a\lambda ,+\mathrm{\infty })上是常数；{\rho }^{\text{'}}(t)存在，且只在有限个t上不连续，{\rho }^{\text{'}}(\mathrm{0}+)=\mathrm{1}；

条件2.4　调节参数\lambda满足

\lambda \gg {\left|{\mathrm{𝒢}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\right|}^{-\mathrm{1}}\mathrm{m}\mathrm{a}\mathrm{x}\{p/n_{\mathrm{0}}, p/n_{\mathrm{1}}, \sqrt[]{\mathrm{l}\mathrm{n}n}/n_{\mathrm{1}}\}.

条件2.1给出参数空间和参数真值的限制条件；条件2.2给出协变量\mathit{X}的限制条件（Amemiya^［9］）， 是参数估计大样本性质证明中常见的必要条件；在高维场景下， 通常假设条件2.3成立， 本文使用的SCAD和MCP惩罚函数满足条件2.3；条件2.4给出调节参数\lambda的限制条件， 在证明估计量的Oracle性质中， 用于对目标函数进行放缩。

引理2.1 （Amemiya^［9］） 在条件2.1和条件2.2下，若K=o(n)，p=o(n)，则式（13）的Oracle估计量{\widehat{\mathit{\theta }}}^{\mathrm{o}\mathrm{r}}是真值{\mathit{\theta }}^{\mathrm{*}}的强相合估计，并且

\sqrt[]{n}({\widehat{\mathit{\theta }}}^{\mathrm{o}\mathrm{r}}-{\mathit{\theta }}^{\mathrm{*}})\to N(\mathrm{0},\mathit{\Sigma }),（16）

其中， \mathit{\Sigma }=\underset{n\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}{\left[-{\nabla }^{\mathrm{2}}\mathrm{l}\mathrm{n}{L}_n({\mathit{\theta }}^{\mathrm{*}})/n\right]}^{-\mathrm{1}}。

定理2.1 在引理2.1下，有

（17）

其中，C=q\sqrt[]{\mathrm{2}{\lambda }_{\mathrm{m}\mathrm{a}\mathrm{x}}/\mathrm{\pi }}，q=K+p+\mathrm{1}，{\varphi }_n=\mathrm{1}/\sqrt[]{\mathrm{l}\mathrm{n}n}。

定理2.2 在定理2.1以及条件2.3和条件2.4下，假设K\ge \mathrm{2}。若b_n>a\lambda且\lambda \gg {\varphi }_n，其中，a是条件2.3给定的常数，{\varphi }_n在定理2.1中给出，那么对式（3）Q_n(\mathit{\alpha },\mathit{\delta },\gamma )，存在一个局部极小值点\widehat{\mathit{\theta }}=({\widehat{\mathit{\alpha }}}^{\mathrm{T}},{\widehat{\mathit{\delta }}}^{\mathrm{T}}{,\widehat{\gamma })}^{\mathrm{T}} 满足

P(\widehat{\mathit{\theta }}={\widehat{\mathit{\theta }}}^{\mathrm{o}\mathrm{r}})\to \mathrm{1}.（18）

3 模拟研究

从如下模型生成未删失响应变量

y_i^{\mathrm{*}}={\mu }_i+{\mathit{x}}_i^{\mathrm{T}}\mathit{\beta }+{ϵ}_i,i=\mathrm{1},\cdots ,n,（19）

其中，{\mathit{x}}_i=(x_{i\mathrm{1}},\cdots ,x_{i\mathrm{5}}{)}^{\mathrm{T}}从均值\mathrm{0}、方差1、协方差为0.3的多元正态分布生成，独立同分布的误差项{ϵ}_i~N(\mathrm{0,0}.{\mathrm{5}}^{\mathrm{2}})。假设响应变量观测值等概率地来自K=\mathrm{3}组。每组截距项设定为P({\tau }_{\mathrm{1}}=-\phi )=P({\tau }_{\mathrm{2}}=\mathrm{0})=P({\tau }_{\mathrm{3}}=\phi )=\mathrm{1}/\mathrm{3}，其中，\phi可取1，1.5和2；真实的系数向量\mathit{\beta }={(\mathrm{5,1},-\mathrm{2,0.5,0.1})}^{\mathrm{T}}。为得到删失比率为q的响应变量观测值y_i，对y_i^{\mathrm{*}}进行如下左删失操作：记\{y_i^{\mathrm{*}},i=\mathrm{1},\dots ,n\}的q分位数为c_q，删去在c_q左侧的y_i^{\mathrm{*}}，即得y_i=\mathrm{m}\mathrm{a}\mathrm{x}\{y_i^{\mathrm{*}},c_q\}。删失比率q可选20%和40%。

为得到目标函数调节参数\lambda的最优选择，本节使用针对高维数据提出的改进版BIC^［25］

\begin{array}{l}\mathrm{B}\mathrm{I}\mathrm{C}(\lambda )=-\mathrm{2}\mathrm{l}\mathrm{n}{L}_n(\widehat{\mathit{\alpha }}(\lambda ),\widehat{\mathit{\delta }}(\lambda ),\widehat{\gamma }(\lambda ))+\\ C_n\frac{\mathrm{l}\mathrm{n}n}{n}(\widehat{K}(\lambda )+p),\end{array}（20）

其中，C_n是依赖于n的正数。当解释变量个数d随着样本量增大而增大时，取C_n=\mathrm{l}\mathrm{n}(\mathrm{l}\mathrm{n}d)，其中， d=n+p ^［26］。本节取C_n=c\mathrm{l}\mathrm{n}(\mathrm{l}\mathrm{n}d)，并选择c=\mathrm{10}^［18］。

ADMM-GCD算法需要输入初始值才能启动，而初始值会影响算法是否收敛以及收敛速度，因此，一个重要的问题是选择恰当的初始值。Yan等^［27］受到文献［28］中Buckley-James迭代过程的启发，对成对截距项之差做L₂惩罚，用最终迭代得到的收敛值作为模拟初始值。应用到本节， 有({\mathit{\alpha }}^{(m+\mathrm{1})},{\mathit{\delta }}^{(m+\mathrm{1})},{\gamma }^{(m+\mathrm{1})})=

其中，\tilde{\mathit{\alpha }},\tilde{\mathit{\delta }},\tilde{\gamma }为当前值，{\lambda }^{\mathrm{*}}=\mathrm{0.001}。Wu等^［17］提出Logistic-Cox模型，把有限混合模型推广到删失数据。本节结合这两种方法选择初值。

本节从两方面评估模型的性能。1）参数估计的精确度：计算\widehat{\mathit{\mu }},\widehat{\mathit{\beta }}和\widehat{\sigma }的均方根误差（root mean squared error， RMSE）， 即{‖\widehat{\mathit{\mu }}-{\mathit{\mu }}^{\mathrm{*}}‖}_{\mathrm{2}}/\sqrt[]{n}，{‖\widehat{\mathit{\beta }}-{\mathit{\beta }}^{\mathrm{*}}‖}_{\mathrm{2}}/\sqrt[]{p}以及平均绝对误差（mean absolute error， MAE），即{\sum }_{i=\mathrm{1}}^n|{\widehat{\mu }}_i-{\mu }_i^{\mathrm{*}}|/n，{\sum }_{j=\mathrm{1}}^p|{\widehat{\beta }}_j-{\beta }_j^{\mathrm{*}}|/p，|\widehat{\sigma }-{\sigma }^{\mathrm{*}}|。2）子群识别的准确性：记录模拟结果中分组个数的估计值\widehat{K}，并计算\widehat{K}的正确估计率P(\widehat{K}=K)以及兰德指数（Rand index， RI）和调整兰德指数（adjusted Rand index， ARI）

\mathrm{R}\mathrm{I}=\frac{\mathrm{T}\mathrm{P}+\mathrm{T}\mathrm{N}}{\mathrm{T}\mathrm{P}+\mathrm{F}\mathrm{P}+\mathrm{F}\mathrm{N}+\mathrm{T}\mathrm{N}},

\mathrm{A}\mathrm{R}\mathrm{I}=\frac{\mathrm{2}(\mathrm{T}\mathrm{P}\times \mathrm{T}\mathrm{N}-\mathrm{F}\mathrm{P}\times \mathrm{F}\mathrm{N})}{(\mathrm{T}\mathrm{P}+\mathrm{F}\mathrm{N})(\mathrm{F}\mathrm{N}+\mathrm{T}\mathrm{N})+(\mathrm{T}\mathrm{P}+\mathrm{F}\mathrm{P})(\mathrm{F}\mathrm{P}+\mathrm{T}\mathrm{N})},

其中，真阳性（TP）决策意味着来自同一个真实子群观测值分到同一组，真阴性（TN）决策表示来自不同真实子群观测值分到不同的组，假阳性（FP）决策代表来自不同真实子群观测值分到同一组，而假阴性（FN）决策是来自同一个真实子群观测值被分到不同的组。RI和ARI的值在0～1之间，越接近1表示分组准确性越高。

在上述6种模型设定下，分别使用3种不同的惩罚函数SCAD、MCP和L₁。本节为每个模型设置重复模拟500次，记录每一种模拟的分组效果和参数估计指标的平均值和标准差，并与Oracle情况下的子群识别及参数估计结果做比较。具体结果见表1和表2。

根据表1和表2，本文所提方法在SCAD惩罚和MCP惩罚下的模拟结果有以下趋势：在相同的截距和惩罚函数下，删失比率越小，参数估计值的RMSE和MAE越小；在相同的删失比率和惩罚函数下，截距越大，识别出的子群数目越接近真实值，且衡量分组准确性的RI、ARI越大，说明子群识别效果越好。这些结果说明本文提出的方法是可行、有效的。使用SCAD和MCP惩罚都能较好地达到识别子群及估计参数的目标，然而，使用L₁惩罚无法识别子群。

接着在K=\mathrm{3},q=\mathrm{20}\mathrm{\%}的情况下，假设随机误差项服从混合正态分布，其密度函数为f(x)=\mathrm{0.1}\varphi (x|\mathrm{0,1})+\mathrm{0.9}\varphi (x|\mathrm{0,0}.{\mathrm{5}}^{\mathrm{2}})， 其中\varphi (·|\mu ,{\sigma }^{\mathrm{2}}) 为正态分布N(\mu ,{\sigma }^{\mathrm{2}})的密度函数， 考察SCAD惩罚和MCP惩罚的表现。具体结果见表3和表4。

表3和表4说明，本文所提方法在非正态误差项下也适用，能够正确识别子群效应，参数估计偏差较小。

最后进行K=\mathrm{5},q=\mathrm{20}\mathrm{\%}下的模拟。具体结果见表5和表6。表5和表6说明，本文所提方法在K=\mathrm{5}时也适用，能够正确识别子群效应，参数估计偏差较小。

4 实例分析

人类免疫缺陷病毒（human immunodeficiency virus， HIV）具有复制周期短、突变比率高的特点，使HIV感染患者在治疗过程中会迅速产生耐药性突变。本节数据来自斯坦福大学HIV耐药性数据库中艾滋病临床试验组进行的OPTIONS研究^［29］。数据由n=\mathrm{407}名参与者的12周随访评估组成。其中，q=\mathrm{35.6}\mathrm{\%}参与者的病毒载量低于或等于检测下限，无法被检测到；n_{\mathrm{1}}=\mathrm{360}名中度耐药性的参与者组成控制组，被随机分配使用或不使用核苷逆转录酶抑制剂；n_{\mathrm{2}}=\mathrm{47}名具有重度耐药性的参与者作为对照组，都使用了核苷逆转录酶抑制剂。

我们使用Jacobson和Zou^［30］筛选出的3个变量作为解释变量，分别是缬氨酸取代逆转录酶基因第 184位的蛋氨酸突变、基线病毒载量以及参与者是否正在服用一种整合酶链转移抑制剂（raltegravir， RAL），以 ln（HIV基线载量） 作为响应变量，使用SCAD和MCP惩罚识别OPTIONS研究中的控制组和对照组，并将结果与实际数据中的组别进行比较。具体结果见表7和表8。

根据表7， 使用SCAD和MCP惩罚均正确识别出2个子群，RI和ARI分别为0.924和0.771。

表8表明， 使用SCAD和MCP惩罚得到的参数估计量值与Oracle参数估计值接近。注意到，{\beta }_{\mathrm{2}}衡量的是RAL与病毒载量的关系， RAL这种抑制剂是目前抗病毒疗效较好的药物。理论上，使用了RAL的患者病毒载量会减少，对应本节的分析结果{\widehat{\beta }}_{\mathrm{2}}为负值。

最后比较SCAD和MCP惩罚的预测表现。把原始数据随机分成366个观测值的训练集和41个观测值的测试集，使用分层抽样确保训练集和测试集有相近的删失比率。重复此过程50次。在50个训练集上，分别使用SCAD、MCP惩罚和Oracle方法进行10折交叉验证，计算Tobit损失值和交叉验证的均方误差（mean square error，MSE）。具体结果见表9。文献［30］对OPTIONS数据集进行分析得到的Tobit损失是1.45， 由表9可知，本文所用方法的Tobit损失更小。

5 总结

本文提出一种惩罚似然推断方法来研究含有潜在子群的左删失数据。基于Tobit回归模型，用回归截距项解释子群效应，将成对融合惩罚的正则化方法拓展到Tobit似然框架。该方法可以同步进行子群识别和参数估计，无需事先确定子群的个数。理论分析和数值研究表明所提方法的有效性。本文所提出的方法可以进一步推广到其他类型删失数据，成对融合惩罚方法也可以应用到回归系数存在子群结构的情形^［19］，所考虑的回归模型也可以推广到半参数模型等。

参考文献

[1] Tobin J. Estimation of relationships for limited dependent variables［J］. Econometrica， 1958， 26（1）： 24. DOI： 10.2307/1907382 .

[2] Buckley J， James I. Linear regression with censored data［J］. Biometrika， 1979， 66（3）： 429-436. DOI： 10.1093/biomet/66.3.429 .

[3] Miller R G. Least squares regression with censored data［J］. Biometrika， 1976， 63（3）： 449-464. DOI： 10.1093/biomet/63.3.449 .

[4] Wei L J， Ying Z， Lin D Y. Linear regression analysis of censored survival data based on rank tests［J］. Biometrika， 1990， 77（4）： 845-851. DOI： 10.1093/biomet/77.4.845 .

[5] Prentice R L. Linear rank tests with right censored data［J］. Biometrika， 1978， 65（1）： 167-179. DOI： 10.1093/biomet/65.1.167 .

[6] 陈菲菲， 孙志华， 叶雪. 失效信息随机缺失时可加危险率模型的统计推断［J］. 中国科学院大学学报， 2016， 33（4）： 443-453. DOI： 10.7523/j．issn．2095-6134.2016.04.003 .

[7] Cuzick J. Asymptotic properties of censored linear rank tests［J］. The Annals of Statistics， 1985， 13（1）： 133-141. DOI： 10.1214/aos/1176346581 .

[8] Fay M P， Shaw P A. Exact and asymptotic weighted logrank tests for interval censored data： the interval r package［J］. Journal of Statistical Software， 2010， 36（2）： 1-34. DOI： 10.18637/jss.v036.i02 .

[9] Amemiya T. Regression analysis when the dependent variable is truncated normal［J］. Econometrica， 1973， 41（6）： 997. DOI： 10.2307/1914031 .

[10] Olsen R J. Note on the uniqueness of the maximum likelihood estimator for the Tobit model［J］. Econometrica， 1978， 46（5）： 1211. DOI： 10.2307/1911445 .

[11] Amemiya T. Tobit models： a survey［J］. Journal of Econometrics， 1984， 24（1/2）： 3-61. DOI： 10.1016/0304-4076（84）90074-5 .

[12] Turnbull B W. The empirical distribution function with arbitrarily grouped， censored and truncated data［J］. Journal of the Royal Statistical Society Series B： Statistical Methodology， 1976， 38（3）： 290-295. DOI： 10.1111/j.2517-6161.1976.tb01597.x .

[13] Zhou M. Empirical likelihood ratio with arbitrarily censored/truncated data by EM algorithm［J］. Journal of Computational and Graphical Statistics， 2005， 14（3）： 643-656. DOI： 10.1198/106186005x59270 .

[14] Everitt B S， Hand D J. Finite mixture distributions［M］. Dordrecht： Springer Netherlands， 1981. DOI： 10.1007/978-94-009-5897-5 .

[15] Wei S S， Kosorok M R. Latent supervised learning［J］. Journal of the American Statistical Association， 2013， 108（503）： 957-970. DOI： 10.1080/01621459.2013.789695 .

[16] Shen J， He X M. Inference for subgroup analysis with a structured logistic-normal mixture model［J］. Journal of the American Statistical Association， 2015， 110（509）： 303-312. DOI： 10.1080/01621459.2014.894763 .

[17] Wu R F， Zheng M， Yu W. Subgroup analysis with time-to-event data under a logistic-cox mixture model［J］. Scandinavian Journal of Statistics， 2016， 43（3）： 863-878. DOI： 10.1111/sjos.12213 .

[18] Ma S J， Huang J. A concave pairwise fusion approach to subgroup analysis［J］. Journal of the American Statistical Association， 2017， 112（517）： 410-423. DOI： 10.1080/01621459.2016.1148039 .

[19] Ma S J， Huang J， Zhang Z W， et al. Exploration of heterogeneous treatment effects via concave fusion［J］. The International Journal of Biostatistics， 2020， 16（1）： arXiv： 1607.03717. DOI： 10.1515/ijb-2018-0026 .

[20] Chen J X， Tran-Dinh Q， Kosorok M R， et al. Identifying heterogeneous effect using latent supervised clustering with adaptive fusion［J］. Journal of Computational and Graphical Statistics： a Joint Publication of American Statistical Association， Institute of Mathematical Statistics， Interface Foundation of North America， 2021， 30（1）： 43-54. DOI： 10.1080/10618600.2020.1763808 .

[21] 许赵辉， 郑泽敏， 吴捷. 基于加速失效模型的子群分析方法［J］. 应用概率统计， 2023， 39（5）： 765-780. DOI： 10.3969/j.issn.1001-4268.2023.05.010 .

[22] Fan J Q， Li R Z. Variable selection via nonconcave penalized likelihood and its oracle properties［J］. Journal of the American Statistical Association， 2001， 96（456）： 1348-1360. DOI： 10.1198/016214501753382273 .

[23] Zhang C H. Nearly unbiased variable selection under minimax concave penalty［J］. The Annals of Statistics， 2010， 38（2）： 894-942. DOI： 10.1214/09-aos729 .

[24] Yang Y， Zou H. An efficient algorithm for computing the HHSVM and its generalizations［J］. Journal of Computational and Graphical Statistics， 2013， 22（2）： 396-415. DOI： 10.1080/10618600.2012.680324 .

[25] Wang H S， Li R Z， Tsai C L. Tuning parameter selectors for the smoothly clipped absolute deviation method［J］. Biometrika， 2007， 94（3）： 553-568. DOI： 10.1093/biomet/asm053 .

[26] Wang H S， Li B， Leng C L. Shrinkage tuning parameter selection with a diverging number of parameters［J］. Journal of the Royal Statistical Society Series B： Statistical Methodology， 2009， 71（3）： 671-683. DOI： 10.1111/j.1467-9868.2008.00693.x .

[27] Yan X D， Yin G S， Zhao X Q. Subgroup analysis in censored linear regression［J］. Statistica Sinica， 2021： 1027-1054. DOI： 10.5705/ss.202018.0319 .

[28] Miller R， Halpern J. Regression with censored data［J］. Biometrika， 1982， 69（3）： 521-531. DOI： 10.1093/biomet/69.3.521 .

[29] Gandhi R T， Tashima K T， Smeaton L M， et al. Long-term outcomes in a large randomized trial of HIV-1 salvage therapy： 96-week results of AIDS clinical trials group a5241 （options）［J］. The Journal of Infectious Diseases， 2020， 221（9）： 1407-1415. DOI： 10.1093/infdis/jiz281 .

[30] Jacobson T， Zou H. High-dimensional censored regression via the penalized tobit likelihood［J］. Journal of Business & Economic Statistics， 2024， 42（1）： 286-297. DOI： 10.1080/07350015.2023.2182309 .

基金资助

国家自然科学基金(12171450)