A surrogate ℓ₀ sparse Cox's regression with applications to sparse high-dimensional massive sample size time-to-event data
一个代理 l ₀ 稀疏 Cox 回归，用于稀疏高维大规模样本量事件发生时间数据

1 INTRODUCTION 1 引言

Advancements in medical informatics tools and high-throughput biological experimentation are making large-scale data routinely accessible to researchers, administrators, and policymakers. This data deluge poses new challenges and critical barriers for quantitative researchers as existing statistical methods and software grind to a halt when analyzing these large-scale data sets, and calls for appropriate methods that can readily fit large-scale data. This paper primarily concerns survival analysis of sparse high-dimensional massive sample size (sHDMSS) data, a particular type of large-scale data with the following characteristics: (1) high-dimensional with a large number of covariates (p_n in thousands or tens of thousands), (2) massive in sample-size (n in thousands to hundreds of millions), (3) sparse in covariates with only a very small portion of covariates being nonzero for each subject, and (4) rare in event rate. An example of sHDMSS data is the pediatric trauma mortality data from the National Trauma Data Bank (NTDB) maintained by the American College of Surgeons.1 This data set includes 210 555 patient records of injured children under 15 collected over 5 years from 2006 to 2010. Each patient record includes 125 952 binary covariates that indicate the presence or absence of an attribute (ICD9 Codes, AIS codes, etc) as well as their two-way interactions. The data matrix is extremely sparse with less than 1% of the covariates being non zero. The event (mortality) rate is also very low at 2%. Another application domain where sHDMSS data are common is drug safety studies that use massive patient-level databases such as the US FDA's Sentinel Initiative (https://www.fda.gov/safety/fdassentinelinitiative/ucm2007250.htm) and the Observational Health Data Sciences and Informatics program (https://ohdsi.org/) to study rare adverse events with hundreds of millions of patient records and tens of thousands of patient attributes that are sparse in the covariates.
医学信息学工具和高通量生物实验的进步使研究人员、管理人员和政策制定者能够经常访问大规模数据。这种数据洪流给定量研究人员带来了新的挑战和关键障碍，因为现有的统计方法和软件在分析这些大规模数据集时停滞不前，并需要适当的方法来适应大规模数据。本文主要涉及稀疏高维大规模样本量 （sHDMSS） 数据的生存分析，这是一种具有以下特征的特殊类型的大规模数据：（1） 具有大量协变量的高维数据（p _n 以数千或数万计），（2） 样本量很大（n 以千万计到数亿），（3） 协变量稀疏，每个受试者只有很小一部分协变量为非零， （4）罕见事件发生率。sHDMSS 数据的一个例子是来自美国外科医师学会维护的国家创伤数据库 （NTDB） 的儿科创伤死亡率数据。1 本数据集包括2006年至2010年5年间收集的210 555份15岁以下受伤儿童的病历。每个患者记录包括 125 952 个二元协变量，这些协变量表示属性（ICD9 代码、AIS 代码等）的存在与否以及它们的双向交互。数据矩阵非常稀疏，只有不到 1% 的协变量不为零。事件（死亡率）也非常低，为2%。sHDMSS数据常见的另一个应用领域是使用大量患者级数据库的药物安全性研究，例如美国FDA的Sentinel Initiative（https://www.fda.gov/safety/fdassentinelinitiative/ucm2007250。htm）和观察性健康数据科学和信息学计划（https://ohdsi.org/），以研究具有数亿份患者记录和数以万计的患者属性的罕见不良事件，这些属性在协变量中是稀疏的。

The sHDMSS survival data present multiple challenges to quantitative researchers. First, not all of the thousands of covariates are expected to be relevant to an outcome of interest. It would also be practically undesirable to predict a patient outcome using thousands of covariates. Traditionally, researchers hand-pick subject characteristics to include in an analysis. However, hand picking can introduce not only bias, but also a source of variability between researchers and studies. Moreover, it would become impractical in large-scale evidence generation when hundreds or thousands of analyses are to be performed.2 Hence, automated sparse regression methods are desired. Secondly, the commonly used “divide and conquer” strategy for massive size data is deemed inappropriate for sHDMSS time-to-event data since each of the divided data would have too few events for a meaningful analysis. Third, sHDMSS data presents a critical barrier to the application of existing sparse survival regression methods, since most current methods and standard software become inoperable for large data sets due to high computational costs and large memory requirements. Although many sparse survival regression methods are available,3-10 to the best of our knowledge, only LASSO, Elastic Net11 and ridge regression have been adapted to fit sHDMSS time-to-event data. In particular, Mittal et al12 developed a tool, named CYCLOPS, for fitting LASSO and ridge Cox regression with sHDMSS time-to-event data by storing data in a sparse format, exploiting sparsity in the data and partial likelihood, and using multicore threading and vector processing, along with other high-performance computing techniques, which delivers > 10-fold speedup12 over its competitors. However, ridge Cox regression does not yield a sparse model and LASSO tends to select too many noise features and is biased for estimation.13, 14 Improved sparse Cox regression tools for sHDMSS time-to-event data are desired.
sHDMSS生存数据给定量研究人员带来了多重挑战。首先，并非所有数千个协变量都与感兴趣的结果相关。使用数千个协变量预测患者结果实际上也是不可取的。传统上，研究人员会亲自挑选要包含在分析中的受试者特征。然而，手工挑选不仅会带来偏见，还会成为研究人员和研究之间差异的来源。此外，当要进行数百或数千次分析时，在大规模证据生成中将变得不切实际。2 因此，需要自动稀疏回归方法。其次，对于大规模数据，常用的“分而治之”策略被认为不适用于 sHDMSS 事件时间数据，因为每个划分的数据都具有太少的事件，无法进行有意义的分析。第三，sHDMSS数据对现有稀疏生存回归方法的应用构成了一个关键的障碍，因为由于计算成本高和内存需求大，大多数当前方法和标准软件都无法用于大型数据集。尽管有许多稀疏生存回归方法可用，但据我们所知，只有 3-10 种，只有 LASSO、Elastic Net 11 和 ridge 回归经过调整以拟合 sHDMSS 事件发生时间数据。特别是，Mittal 等人 12 开发了一种名为 CYCLOPS 的工具，通过以稀疏格式存储数据、利用数据的稀疏性和部分似然性、使用多核线程和矢量处理以及其他高性能计算技术，将 LASSO 和 ridge Cox 回归与 sHDMSS 事件发生时间数据拟合，>竞争对手的加速速度提高了 10 倍 12。 然而，ridge Cox回归不会产生稀疏模型，并且LASSO倾向于选择过多的噪声特征，并且在估计时存在偏差。13， 14 需要改进的稀疏 Cox 回归工具，用于 sHDMSS 事件发生时间数据。

The purpose of this paper is to develop a surrogate ℓ₀-based sparse Cox regression method and adapt it to sHDMSS time-to-event data. It is well known that ℓ₀-penalized regression is natural for variable selection and parameter estimation with some optimal properties.15-18 On the other hand, it is also known to have some pitfalls such as instability19 and being unscalable to even moderate dimensional covariates. The broken adaptive ridge (BAR) estimator, defined as the limit of an iteratively reweighted ℓ₂-penalization algorithm, was introduced to approximate the ℓ₀-penalization problem and has been recently shown to possess some desirable selection, estimation, and clustering properties under the linear model and several other model settings.10, 20-22 It is also computationally scalable to high-dimensional covariates and stable for variable selection as discussed later in Remark 2 of Section 2. However, the BAR method has yet to be rigorously studied for the Cox model. Moreover, current BAR algorithms have only been implemented for densely-represented covariates and are unsuitable for sHDMSS data due to high computational costs, high memory requirements, and numerical instability. Computation of the Cox partial likelihood and its derivatives is particularly demanding for massive sample size data since the required number of operations grows at the rate of O(n²). The key contributions of this paper are twofold. First, we rigorously extend the BAR methodology to the Cox model. Specifically, we establish the selection consistency, an oracle property for parameter estimation, and a grouping property of highly correlated covariates for the Cox model. It is worth noting that the theoretical extension of the BAR methodology to Cox model is nontrivial and notably different from other models because the log-partial likelihood for the Cox model is not the sum of independent terms and the standard martingale central limit theorem used to derive the asymptotic theory for Cox's model with a fixed number of covariates is no longer applicable when the number of parameters diverges. Furthermore, because BAR involves performing an infinite number of penalized regressions, the derivations of its selection consistency and oracle property for estimation are substantially different from those for a single-step oracle estimator in the literature. The second key contribution of this paper is to develop an efficient implementation of BAR for Cox regression with sHDMSS time-to-event data by leveraging existing efficient massive ℓ₂-penalized Cox regression techniques,12 which include employing a column relaxation with logistic loss (CLG) algorithm using one-dimensional updates and a one-step Newton-Raphson approximation as well as exploiting the sparsity in the covariate structure and the Cox partial likelihood to reduce the number of operations from O(n²) to O(n).
本文旨在开发一种基于代理 l ₀ 的稀疏 Cox 回归方法，并使其适应 sHDMSS 事件发生时间数据。众所周知，l ₀ 惩罚回归对于具有某些最优属性的变量选择和参数估计是很自然的。15- 18 另一方面，众所周知，它也存在一些缺陷，例如不稳定性 19 并且无法扩展到甚至中等维协变量。断开自适应岭 （BAR） 估计器定义为迭代重加权 l ₂ 惩罚算法的极限，用于近似 l ₀ 惩罚化问题，最近已被证明具有一些理想的选择、估计和聚类属性在线性模型和其他几个模型设置下。10， 20- 22 它还可以在计算上扩展到高维协变量，并且对于变量选择是稳定的，如第 2 节后面的注释 2 所述。然而，BAR方法尚未针对Cox模型进行严格的研究。此外，目前的BAR算法只针对密集表示的协变量实现，由于计算成本高、内存要求高、数值不稳定，不适合sHDMSS数据。Cox 偏似然及其导数的计算对于大量样本量数据的要求特别高，因为所需的运算数量以 O（n ² ） 的速率增长。本文的主要贡献是双重的。首先，我们将BAR方法严格扩展到Cox模型。具体来说，我们为 Cox 模型建立了选择一致性、用于参数估计的预言机属性以及高度相关协变量的分组属性。 值得注意的是，BAR方法对Cox模型的理论扩展是非平凡的，并且与其他模型明显不同，因为Cox模型的对数偏似然不是独立项的总和，并且当参数数量发散时，用于推导具有固定数量协变量的Cox模型的渐近理论的标准马丁格尔中心极限定理不再适用。此外，由于 BAR 涉及执行无限数量的惩罚回归，因此其用于估计的选择一致性和预言机属性的推导与文献中单步预言机估计器的推导有很大不同。本文的第二个主要贡献是通过利用现有的高效大规模 l ₂ 惩罚 Cox 回归技术 12，开发一种使用 sHDMSS 事件发生时间数据的 BAR 的高效实现，其中包括使用一维更新和一步牛顿-拉夫森近似的逻辑损失列松弛 （CLG） 算法，以及利用协变量结构中的稀疏性和 Cox 偏似然来降低从 O（n ² ） 到 O（n） 的操作次数。

In Section 2, we formally define the BAR estimator, state its theoretical properties for variable selection, parameter estimation, and grouping highly correlated covariates for the Cox model, and describe an efficient implementation for sHDMSS survival data. We also discuss how to adapt BAR as a postscreening sparse regression method for ultrahigh dimensional Cox regression with relatively small sample size. In Section 3, we present simulation studies to demonstrate the performance of the CoxBAR estimator with both moderate and massive sample size in various low and high-dimensional settings. We provide a real data example using the pediatric trauma mortality data12 in Section 4. Lastly, we give closing remarks in Section 5. The appendix collects proofs of the theoretical results and regularity conditions needed for the derivations. An R package has been developed for BAR and made available at https:github.com/OHDSI/BrokenAdaptiveRidge.
在第 2 节中，我们正式定义了 BAR 估计器，陈述了其对变量选择、参数估计和 Cox 模型高度相关协变量分组的理论属性，并描述了 sHDMSS 生存数据的有效实现。我们还讨论了如何将 BAR 作为筛选后稀疏回归方法，用于样本量相对较小的超高维 Cox 回归。在第 3 节中，我们提出了仿真研究，以证明 CoxBAR 估计器在各种低维和高维设置下具有中等和大样本量的性能。我们在第 4 节中使用儿科创伤死亡率数据 12 提供了一个真实的数据示例。最后，我们在第 5 节中作结束语。附录收集了推导所需的理论结果和正则性条件的证明。已为 BAR 开发了一个 R 包，可在 https:github.com/OHDSI/BrokenAdaptiveRidge 上获得。

2 METHODOLOGY 2 方法论

2.1 Cox's BAR regression and its large sample properties
2.1 Cox's BAR回归及其大样本性质

2.1.1 The data structure, model, and estimator
2.1.1 数据结构、模型和估计器

Suppose that one observes a random sample of right-censored time-to-event data consisting of n independent and identically distributed triplets ${\{(X_i,{\delta }_i,{\mathbf{z}}_i(·))\}}_{i=1}^n$, where for subject i, $X_i=\min (T_i,C_i)$ is the observed event time, δ_i=I(T_i ≤ C_i) is the censoring indicator, T_i is the event time of interest, and C_i is a censoring time that is conditionally independent of T_i given a p_n-dimensional, possibly time-dependent, covariate vector ${\mathbf{z}}_i(·)={(z_{i1}(·),\text{…},z_{i{p}_n}(·))}^{\prime }$.
假设观察到一个由 n 个独立且分布相同的三元 ${\{(X_i,{\delta }_i,{\mathbf{z}}_i(·))\}}_{i=1}^n$ 组组成的右删失事件时间数据的随机样本，其中对于受试者 i， $X_i=\min (T_i,C_i)$ 是观察到的事件时间，δ=I（T ≤ C） 是删失指标，T 是感兴趣的事件时间，C 是条件独立于 T 的删失时间，给定 p _n 维， 可能与时间相关的协变量向量 ${\mathbf{z}}_i(·)={(z_{i1}(·),\text{…},z_{i{p}_n}(·))}^{\prime }$ 。

Assume the Cox23 proportional hazard model
假设 Cox 23 比例风险模型

$$h\left\{t|\mathbf{z}(t)\right\}=h_0(t)\exp \left\{\mathbf{z}{(t)}^{\prime }\mathit{\beta }\right\},$$ (1)

where h{t | z(t)} is the conditional hazard function of T_i given {z(u), 0 ≤ u ≤ t,}, h₀(t) is an unspecified baseline hazard function, and $\mathit{\beta }=({\beta }_1,\text{…},{\beta }_{p_n})$ is a vector of time-independent regression coefficients. Denote by β₁ and β₂ the first q_n and remaining p_n−q_n components of β, respectively, and define ${\mathit{\beta }}_0={\left({\mathit{\beta }}_{01}^{\prime },{\mathit{\beta }}_{02}^{\prime }\right)}^{\prime }$ as the true values of β, where, without loss of generality, ${\mathit{\beta }}_{01}=({\beta }_{01}\text{…},{\beta }_{0q_n})$ is a vector of q_n nonzero values and β₀₂=0 is a p_n−q_n dimensional vector of zeros. Further technical assumptions for β₀ and p_n are given later in condition (C6) of Section S4 of the Supplementary Material. For simplicity, we work on the time interval s∈[0, 1] as in the work of Andersen and Gill,24 which can be extended to any time interval [0, τ] for 0<τ<∞. Using the standard counting process notation, the log-partial likelihood for the Cox model is defined as
其中 h{t | z（t）} 是给定 {z（u）， 0 ≤ u ≤ t，} 的条件风险函数，h ₀ （t） 是未指定的基线风险函数， $\mathit{\beta }=({\beta }_1,\text{…},{\beta }_{p_n})$ 是与时间无关的回归系数的向量。分别用 β ₁ 和 β ₂ 表示 β 的前 q _n 和剩余的 p _n −q _n 分量，并定义为 ${\mathit{\beta }}_0={\left({\mathit{\beta }}_{01}^{\prime },{\mathit{\beta }}_{02}^{\prime }\right)}^{\prime }$ β 的真值，其中在不损失普遍性的情况下， ${\mathit{\beta }}_{01}=({\beta }_{01}\text{…},{\beta }_{0q_n})$ 是 q _n 非零值的向量，β ₀₂ =0 是零的 p _n −q _n 维向量。β ₀ 和 p _n 的进一步技术假设将在补充材料第 S4 节的条件 （C6） 中给出。为简单起见，我们按照 Andersen 和 Gill 的工作 24 处理时间间隔 s∈[0， 1]，它可以扩展到任何时间间隔 [0， τ] 0<τ<∞。使用标准计数过程表示法，Cox 模型的对数部分似然定义为

$$l_n(\mathit{\beta })={\displaystyle \sum _{i=1}^n}\underset{0}{\overset{1}{\int }}{\mathit{\beta }}^{\prime }{\mathbf{z}}_i(s)d{N}_i(s)-\underset{0}{\overset{1}{\int }}\log \left[{\displaystyle \sum _{j=1}^n}{Y}_j(s)\exp \left\{{\mathit{\beta }}^{\prime }{\mathbf{z}}_j(s)\right\}\right]d\overline{N}(s),$$ (2)

where, for subject i, Y_i(s)=I(X_i ≥ s) is the at-risk process and N_i(s)=I(X_i ≤ s,δ_i=1) is the counting process of the uncensored event with intensity process $h_i(t|\mathit{\beta })=h_0(t)Y_i(t)\exp \{{\mathbf{z}}_i{(t)}^{\prime }\mathit{\beta }\}$ and $\overline{N}={\sum }_{i=1}^n{N}_i$.
其中，对于受试者 i，Y（s）=I（X ≥ s） 是风险过程，N（s）=I（X ≤ s，δ=1） 是具有强度过程 $h_i(t|\mathit{\beta })=h_0(t)Y_i(t)\exp \{{\mathbf{z}}_i{(t)}^{\prime }\mathit{\beta }\}$ 的未经审查事件的计数过程，并且 $\overline{N}={\sum }_{i=1}^n{N}_i$ 。

Our Cox's BAR estimation of β starts with an initial Cox ridge regression estimator25
我们对 β 的 Cox BAR 估计从初始 Cox 岭回归估计器 25 开始

$${\widehat{\mathit{\beta }}}^{(0)}=\arg \underset{\mathit{\beta }}{\min }\left\{-2l_n(\mathit{\beta })+{\xi }_n{\displaystyle \sum _{j=1}^{p_n}}{\beta }_j^2\right\},$$ (3)

which is updated iteratively by a reweighed ℓ₂-penalized Cox regression estimator
它由重新加权的 l ₂ 惩罚 Cox 回归估计器迭代更新

$${\widehat{\mathit{\beta }}}^{(k)}=\arg \underset{\mathit{\beta }}{\min }\left\{-2l_n(\mathit{\beta })+{\lambda }_n{\displaystyle \sum _{j=1}^{p_n}}\frac{{\beta }_j^2}{{\left({\widehat{\beta }}_j^{(k-1)}\right)}^2}\right\},k\ge 1,$$ (4)

where ξ_n and λ_n are nonnegative penalization tuning parameters. The BAR estimator is defined as
其中 ξ _n 和 λ _n 是非负惩罚调整参数。BAR估计器定义为

$$\widehat{\mathit{\beta }}=\underset{k\to \infty }{\lim }{\widehat{\mathit{\beta }}}^{(k)}.$$ (5)

Since ℓ₂-penalization yields a nonsparse solution, defining the BAR estimator as the limit is necessary to produce sparsity. Although λ_n is fixed at each iteration, it is weighted inversely by the square of the ridge regression estimates from the previous iteration. Consequently, coefficients whose true values are zero will have larger penalties in the next iteration, whereas penalties for truly nonzero coefficients will converge to a constant. We will show later in Theorem 1 that, under certain regularity conditions, the estimates of the truly zero coefficients shrink toward zero while the estimates of the truly nonzero coefficients converge to their oracle estimates with probability tending to 1. As illustrated by a small simulation in Section S2 (Figure S1) of the Supplementary Material, the signal (nonzero coefficients) and noise (zero coefficients) can be quickly separated within a few BAR iterations, although more iterations may be necessary in some scenarios to improve estimation of the nonzero coefficients.
由于 l ₂ -惩罚化产生非稀疏解，因此将 BAR 估计器定义为极限对于生成稀疏性是必要的。尽管 λ 在每次迭代中 _n 都是固定的，但它的权重与上一次迭代的岭回归估计值的平方成反比。因此，真实值为零的系数在下一次迭代中将受到更大的惩罚，而真正非零系数的惩罚将收敛到一个常数。我们将在定理 1 的后面说明，在某些正则性条件下，真正零系数的估计值向零收缩，而真正非零系数的估计值收敛到它们的预言机估计值，概率趋于 1。如补充材料第S2节（图S1）中的小型仿真所示，信号（非零系数）和噪声（零系数）可以在几次BAR迭代中快速分离，尽管在某些情况下可能需要更多迭代来改进对非零系数的估计。

Remark 1.(Computational aspects of BAR) For moderate size data, one may calculate ${\widehat{\mathit{\beta }}}^{(k)}$ in 4 using the Newton-Raphson method as in the work of Frommlet and Nuel,26 who outlined an iterative reweighted ridge regression for generalized linear models. It appears at first sight that 4 will encounter numerical overflow as some of the coefficients ${\widehat{\beta }}_j^{(k-1)}$ will go to zero as k increases. However, it can be shown that after some simple algebraic manipulation the Newton-Raphson updating formula will only involve multiplications, instead of divisions, by ${\widehat{\beta }}_j^{(k-1)}$s, and thus numerical overflow can be avoided. Further details are provided in Section S1 (Equation (3)) of the Supplementary Material. We also note that because the limit of the BAR algorithm cannot be numerically achieved at any finite iteration step, an extra thresholding rule for small coefficients will be required to numerically obtain a sparse solution. However, this thresholding level can be set arbitrarily small (by default, we set the threshold value to 10⁻⁶ in our implementation) since it is simply used for numerical convergence to zero and has minimal impact on the resulting BAR estimator. Furthermore, Equation (3) of Section S1 of the Supplementary Material implies that, once a ${\widehat{\beta }}_j^{(k-1)}$ becomes zero, it will remain as zero in subsequent iterations. Thus, one only needs to update ${\widehat{\mathit{\beta }}}^{(k)}$ within the reduced nonzero parameter space, an appealing computational advantage for high-dimensional settings.
备注 1.（BAR的计算方面）对于中等大小的数据，可以使用 Newton-Raphson 方法在 4 中进行计算 ${\widehat{\mathit{\beta }}}^{(k)}$ ，就像 Frommlet 和 Nuel 的工作一样，26 他们概述了广义线性模型的迭代重加权岭回归。乍一看，4 似乎会遇到数值溢出，因为随着 k 的增加，一些系数 ${\widehat{\beta }}_j^{(k-1)}$ 将变为零。然而，可以证明，经过一些简单的代数运算后，牛顿-拉夫森更新公式将只涉及乘法，而不是除 ${\widehat{\beta }}_j^{(k-1)}$ 法，乘以 s，因此可以避免数值溢出。补充材料的S1部分（等式（3））提供了更多详细信息。我们还注意到，由于BAR算法的极限在任何有限迭代步骤中都无法数值实现，因此需要对小系数的额外阈值规则进行数值化以获得稀疏解。但是，这个阈值水平可以任意设置（默认情况下，我们在实现 ⁻⁶ 中将阈值设置为 10），因为它只是用于数值收敛到零，并且对生成的 BAR 估计器的影响最小。此外，补充材料 S1 部分的等式 （3） 意味着，一旦 a ${\widehat{\beta }}_j^{(k-1)}$ 变为零，它将在随后的迭代中保持为零。因此，只需要在减少的非零参数空间内进行更新 ${\widehat{\mathit{\beta }}}^{(k)}$ ，这对于高维设置来说是一个有吸引力的计算优势。

For massive-size data with large n and p_n, the Newton-Raphson procedure, which at each iteration, calls for the calculation of both the gradient and Hessian can become practically infeasible due to high computational costs, memory requirements, and numerical instability. In Section 2.2, we will discuss how to adapt an efficient algorithm for massive ℓ₂-penalized Cox regression via cyclic coordinate descent and exploit the sparsity of the covariate structure to make BAR scalable to sHDMSS data.
对于具有大 n 和 p _n 的海量数据，Newton-Raphson 过程在每次迭代中都要求计算梯度和 Hessian 由于计算成本高、内存要求和数值不稳定性而变得几乎不可行。在第 2.2 节中，我们将讨论如何通过循环坐标下降调整一种有效的算法，以实现大规模 l ₂ 惩罚的 Cox 回归，并利用协变量结构的稀疏性使 BAR 可扩展到 sHDMSS 数据。

Remark 2.(Broken adaptive ridge versus best subset selection) The BAR method can be viewed as a performing a sequence of surrogate ℓ₀-penalizations, where each reweighted ℓ₂ penalty serves as a surrogate ℓ₀-penalty and the approximation of ℓ₀-penalization improves with each iteration. Consequently, BAR enjoys the best of ℓ₀- and ℓ₂-penalized regressions. For example, we establish in the next two sections that BAR possesses the oracle properties for estimation and selection consistency (an ℓ₀ property) as well as a grouping property (an ℓ₂ property). Numerically, for a fixed tuning parameter value, BAR is a surrogate to ℓ₀-penalization is not expected to be identical, but can be similar to the exact global ℓ₀-penalization solution where the latter must be solved using the best subset search (BSS). We illustrate this in Section S3 of the Supplementary Material (Figures S2 and S3) using a small simulation study. It is worth emphasizing that BAR overcomes some shortcomings of BSS; for example, BSS is computationally NP-hard and can be unstable for variable selection,19 whereas BAR is scalable to high-dimensional covariates and is more stable for variable selection as demonstrated in Figures S2 and S3 in Section S3 of the Supplementary Material.
备注 2.（破碎的自适应脊与最佳子集选择）BAR 方法可以看作是执行一系列代理 l ₀ 惩罚，其中每个重新加权的 l ₂ 惩罚都充当代理 l ₀ 惩罚，并且 l ₀ 惩罚的近似值随着每次迭代而提高。因此，BAR 享有 l ₀ 和 l ₂ 惩罚回归的最佳效果。例如，在接下来的两节中，我们将确定 BAR 具有用于估计和选择一致性的预言机属性（l ₀ 属性）以及分组属性（l ₂ 属性）。从数值上讲，对于固定的调优参数值，BAR 是 l ₀ 的替代项，预计不会完全相同，但可以类似于精确的全局 l ₀ 惩罚化解决方案，后者必须使用最佳子集搜索 （BSS） 求解。我们在补充材料的S3部分（图S2和S3）中使用一个小型模拟研究来说明这一点。值得强调的是，BAR克服了BSS的一些缺点;例如，BSS 在计算上是 NP 困难的，对于变量选择可能不稳定 19，而 BAR 可扩展到高维协变量，并且对于变量选择更稳定，如补充材料第 S2 节中的图 S3 和 S3 所示。

2.1.2 Oracle properties 2.1.2 Oracle属性

We establish the oracle properties for the BAR estimator for simultaneous variable selection and parameter estimation where we allow both q_n and p_n to diverge to infinity.
我们为 BAR 估计器建立了预言机属性，用于同时进行变量选择和参数估计，其中我们允许 q _n 和 p _n 发散到无穷大。

Theorem 1. (Oracle properties)Assume the regularity conditions (C1) to (C6) in Section S4.1 of the Supplementary Material hold. Let ${\widehat{\mathit{\beta }}}_1$ and ${\widehat{\mathit{\beta }}}_2$ be the first q_n and the remaining p_n−q_n components of the BAR estimator $\widehat{\mathit{\beta }}$, respectively. Then, as n→∞, with probability tending to one,

• (a) the BAR estimator $\widehat{\mathit{\beta }}=({\widehat{\mathit{\beta }}}_1,{\widehat{\mathit{\beta }}}_2)$ exists and is unique, where ${\widehat{\mathit{\beta }}}_2=\mathbf{0}$;
• (b)  （二） $\sqrt{n}{\mathbf{b}}_n^{\prime }\mathrm{\sum }{({\mathit{\beta }}_0)}_{11}^{-1/2}({\widehat{\mathit{\beta }}}_1-{\mathit{\beta }}_{01})\stackrel{D}{\to }N(0,1)$, for any q_n-dimensional vector b_n such that ‖b_n‖₂ ≤ 1 and where ∑(β₀)₁₁ is the first q_n×q_n submatrix of ∑(β₀), where ∑(β₀) is defined in Condition (C4).
  ，对于任何 q _n 维向量 b _n ，使得 ‖b _{n 2} ‖ ≤ 1，并且 ∑（β ₀ ） ₁₁ 是 ∑（β ₀ 的第一个 q _n ×q _n 子矩阵），其中 ∑（β ₀ ） 在条件 （C4） 中定义。

Theorem 1(a) establishes selection consistency of the BAR estimator. Part (b) of the theorem essentially states that the nonzero component of the BAR estimator is asymptotically normal and equivalent to the weighted ridge estimator of the oracle model, as shown in the proof provided in Section S4.2 of the Supplementary Material.
定理1（a）建立了BAR估计器的选择一致性。该定理的 （b） 部分基本上指出，BAR 估计器的非零分量是渐近正态的，等效于预言机模型的加权脊估计器，如补充材料第 S4.2 节中提供的证明所示。

Remark 3. (Ultrahigh-dimensional covariates setting)Although we allow p_n to diverge, the asymptotic properties of the BAR estimator in the Section 28 are derived for p_n<n. In an ultrahigh-dimensional setting where the number of covariates far exceeds the number of observations (p_n≫n), one may couple a sure screening27 method with the BAR estimator to obtain a two-step estimator with desirable selection and estimation properties. The orders of q_n, p_n, and n and their relationships depend on the employed screening procedure. For example, coupling the BAR estimator with the sure joint screening procedure28 has been explored in the work of Kawaguchi.29
备注 3.（超高维协变量设置）尽管我们允许 p _n 发散，但第 28 节中 BAR 估计器的渐近性质是针对 p _n ≫n） 的超高维环境中，可以将确定筛选 27 方法与 BAR 估计器耦合，以获得具有理想选择和估计特性的两步估计器。q _n 、p _n 和 n 的顺序及其关系取决于所采用的筛选程序。例如，在Kawaguchi的工作中已经探索了BAR估计器与确定联合筛选程序28的耦合。29

2.1.3 A grouping property
2.1.3 分组属性

When the true model has a group structure, it is desirable for a variable selection method to either retain or drop all variables that are clustered within the same group. It is well known that ridge regression possesses the grouping property for highly correlated covariates.11 Because the BAR estimator is based on an iterative ridge regression, we show that BAR also possesses a grouping property for highly correlated covariates as stated in following theorem.
当真实模型具有组结构时，变量选择方法最好保留或删除聚类在同一组中的所有变量。众所周知，岭回归具有高度相关协变量的分组特性。11 由于 BAR 估计器基于迭代岭回归，因此我们表明 BAR 还具有高度相关的协变量的分组属性，如以下定理所述。

Theorem 2.Let λ_n, ${\{(X_i,{\delta }_i,{\mathbf{z}}_i)\}}_{i=1}^n$ be given and assume that $Z=({\mathbf{z}}_i^{\prime },\text{…}{\mathbf{z}}_n^{\prime })$ is standardized. That is, for all j=1,…,p_n, ${\displaystyle {\text{∑}}_{i=1}^n}{z}_{ij}=0,{\mathbf{z}}_{[,j]}^{\prime }{\mathbf{z}}_{[,j]}=n-1,$ where z_[, j] is the j^th column of Z. Suppose the regularity conditions (C1) to (C6) in Section S4.1 of the Supplementary Material hold and let $\widehat{\mathit{\beta }}$ be the BAR estimator. Then, for any ${\widehat{\beta }}_i\ne 0$ and ${\widehat{\beta }}_j\ne 0$,
定理 2.设 λ _n ， ${\{(X_i,{\delta }_i,{\mathbf{z}}_i)\}}_{i=1}^n$ 并假设 $Z=({\mathbf{z}}_i^{\prime },\text{…}{\mathbf{z}}_n^{\prime })$ 它是标准化的。也就是说，对于所有 j=1,...,p _n ， ${\displaystyle {\text{∑}}_{i=1}^n}{z}_{ij}=0,{\mathbf{z}}_{[,j]}^{\prime }{\mathbf{z}}_{[,j]}=n-1,$ 其中 z _[, j] 是 Z 的 j ^th 列。 假设补充材料第 S4.1 节中的正则性条件 （C1） 至 （C6） 成立，并设 $\widehat{\mathit{\beta }}$ 为 BAR 估计器。然后，对于任何 ${\widehat{\beta }}_i\ne 0$ 和 ${\widehat{\beta }}_j\ne 0$ ，

$$|{\widehat{\beta }}_i^{-1}\pm {\widehat{\beta }}_j^{-1}|\le \frac{1}{{\lambda }_n}\sqrt{2\left\{(n-1)(1\pm r_{ij})\right\}}\sqrt{n{(1+d_n)}^2},$$ (6)

with probability tending to one, where $d_n={\sum }_{i=1}^n{\delta }_i$, and $r_{ij}=\frac{1}{n-1}{\mathbf{z}}_{[,i]}^{\prime }{\mathbf{z}}_{[,j]}$ is the sample correlation of z_[,i] and z_[, j].
概率趋于 1，其中 $d_n={\sum }_{i=1}^n{\delta }_i$ 和 $r_{ij}=\frac{1}{n-1}{\mathbf{z}}_{[,i]}^{\prime }{\mathbf{z}}_{[,j]}$ 是 z _[,i] 和 z _[, j] 的样本相关性。

The proof is provided in Section S4.3 of the Supplementary Material. It is seen from 6 that, as r_i j→1, the absolute difference between ${\widehat{\beta }}_i$ and ${\widehat{\beta }}_j$ approaches 0, implying that the estimated coefficients of two highly positively correlated variables will be similar in magnitude. Similarly, the estimated coefficients of two highly negatively correlated variables are also similar in magnitude with a sign change.
该证明在补充材料的第 S4.3 节中提供。从图 6 可以看出，当 r _i j →1 时，和 之间的 ${\widehat{\beta }}_i$ ${\widehat{\beta }}_j$ 绝对差值接近 0，这意味着两个高度正相关变量的估计系数在大小上相似。同样，两个高度负相关变量的估计系数在量级上也相似，并具有符号变化。

2.1.4 Selection of tuning parameters
2.1.4 调优参数的选择

Model complexity depends critically on the choice of the tuning parameters. The BAR estimator depends on two tuning parameters, ie, ξ_n for the initial ridge estimator in 3 and λ_n for the iterative ridge step in 4. Our simulations in Section 3.1 illustrate that, while fixing λ_n, the BAR estimator is insensitive to the choice of ξ_n over a wide interval (Figure 1) and thus practically only optimization with respect to λ_n is needed.
模型的复杂性很大程度上取决于调优参数的选择。BAR 估计器依赖于两个调谐参数，即 3 中的初始岭估计器为 ξ，4 _n 中的迭代岭步长为 λ _n 。我们在 3.1 节中的模拟表明，在固定 λ 时 _n ，BAR 估计器对宽区间内 ξ _n 的选择不敏感（图 1），因此实际上只需要对 λ _n 进行优化。

We optimize with respect to λ_n in a similar manner to currently used penalization methods. A popular strategy for tuning parameter selection is to perform optimization with respect to a data-driven selection criterion such as cross-validation (CV),30, 31 Akaike information criterion,15 and Bayesian information criterion (BIC).16, 17, 32 Although CV has been used extensively in the literature, it has been known to asymptotically overfit models with a positive probability.33, 34 Recent theoretical work has shown that, for penalized Cox models that possess the oracle property, BIC-based tuning parameter selection identifies the true model with probability tending to one.32 Further discussion on selecting λ_n for BAR is provided in the last paragraph of Section 3.2.
我们 _n 以与当前使用的惩罚方法类似的方式对 λ 进行优化。调整参数选择的常用策略是针对数据驱动的选择标准（例如交叉验证 （CV）、30、31 赤池信息标准、15 和贝叶斯信息标准 （BIC））执行优化。16， 17， 32 尽管 CV 在文献中被广泛使用，但已知它会以正概率渐近过拟合模型。33， 34 最近的理论研究表明，对于具有预言机属性的惩罚 Cox 模型，基于 BIC 的调参数选择可以识别概率趋向于 1 的真实模型。32 关于 _n 为 BAR 选择 λ 的进一步讨论在第 3.2 节的最后一段中提供。

3 SIMULATIONS 3 模拟

This section presents three simulation studies. First, we demonstrate in Section 3.1 that, for fixed λ_n, the BAR estimator is insensitive to the tuning parameter ξ_n of its initial ridge estimator and does well in terms of performing variable selection and correcting possible bias of the initial ridge estimator. Then, in Section 3.2, we evaluate and compare the operating characteristics of BAR with some popular penalized Cox regression methods, where we only consider settings with moderate sample sizes due to the fact that most of the competing methods are inoperable for massive sample size data. Finally, in Section 3.3, we use a sHDMSS setting to illustrate the performance of BAR over its closest competitor.
本节介绍三个仿真研究。首先，我们在第 3.1 节中证明，对于固定 λ _n ，BAR 估计器对其初始岭估计器的调谐参数 ξ _n 不敏感，并且在执行变量选择和校正初始岭估计器的可能偏差方面表现良好。然后，在第 3.2 节中，我们评估并比较了 BAR 的操作特性与一些流行的惩罚 Cox 回归方法，其中我们只考虑具有中等样本量的设置，因为大多数竞争方法都不适用于大量样本量数据。最后，在第 3.3 节中，我们使用 sHDMSS 设置来说明 BAR 相对于其最接近的竞争对手的性能。

Sections 3.1 and 3.2 employ the same simulation structure. Event times are drawn from an exponential proportional hazards model with baseline hazard h₀(t)=1, ${\mathit{\beta }}_0=(0.40,0,0.45,0,0.50,0.55,0,0,0.70,0.80,{\mathbf{0}}_{p_n-10})$, representing q_n=6 small to moderate effect sizes; the design matrix $Z=({\mathbf{z}}_1^{\prime },\text{…},{\mathbf{z}}_n^{\prime })$ is generated from a p_n-dimensional normal distribution with mean zero and covariance matrix ∑=(σ_i j) with an autoregressive structure such that σ_i j=0.5^|i−j| and independent censoring times are generated from uniform distribution U(0,u_max), where u_max is chosen to achieve different percentages of censoring. We describe how we simulate sHDMSS time-to-event data in Section 3.3.
第 3.1 节和第 3.2 节采用相同的模拟结构。事件时间是从基线风险 h ₀ （t）=1 的指数比例风险模型中得出的， ${\mathit{\beta }}_0=(0.40,0,0.45,0,0.50,0.55,0,0,0.70,0.80,{\mathbf{0}}_{p_n-10})$ 表示 q _n =6 个小到中等效应量;设计矩阵 $Z=({\mathbf{z}}_1^{\prime },\text{…},{\mathbf{z}}_n^{\prime })$ 由均值为零的 p _n 维正态分布和协方差矩阵 ∑=（σ _i j ） 生成，具有自回归结构，使得 σ _i j =0.5 ^|i−j| 和独立删失时间由均匀分布 U（0，u _max ） 生成，其中选择 U _max 以实现不同的删失百分比。我们将在第 3.3 节中描述如何模拟 sHDMSS 事件发生时间数据。

3.3 Sparse high-dimensional massive sample size data
3.3 稀疏高维海量样本量数据

In this simulation, we simulate a sHDMSS time-to-event data set with n=200,000, p_n=20,000, and q_n=80. Event times are generated from an exponential hazards model with baseline hazard h₀(t)=1, regression coefficients ${\mathit{\beta }}_0=({\mathbf{0}\mathbf{.}\mathbf{7}}_{10},{\mathbf{0}\mathbf{.}\mathbf{5}}_{10},{\mathbf{0}\mathbf{.}\mathbf{8}}_{10},{\mathbf{1}}_{10},{\mathbf{-}\mathbf{0}\mathbf{.}\mathbf{7}}_{10},{\mathbf{-}\mathbf{0}\mathbf{.}\mathbf{5}}_{10},{\mathbf{-}\mathbf{0}\mathbf{.}\mathbf{8}}_{10},{\mathbf{-}\mathbf{1}}_{10},{\mathbf{0}}_{p_n-80})$, and a censoring rate of 95%. The covariates for each subject are simulated such, on average, 2% are assigned a nonzero value. The amount of memory used to store this dense design matrix would require over 16 GB, which exceeds the functional capacity of most statistical software packages on standard hardware. To overcome this difficulty, we efficiently store the information in a coordinate list fashion and compare our massive Cox's regression for BAR (mBAR) with the massive sparse Cox's regression for LASSO (mCox-LASSO) using the Cyclops package,12, 38 which, to the best of our knowledge, is the fastest software available today that exploits the sparsity of sHDMSS time-to-event data for efficient computing and offers > 10-fold speedup12 over its competitors such as CoxNet7 and FastCox.47 For LASSO, CV (mCox-LASSO (CV)), combined with a nonconvex optimization technique which is more efficient than the classical grid search approach, and BIC score minimization (mCox-LASSO (BIC)), implemented with the classical grid search approach, were used to find the optimal value for the tuning parameter. For the mBAR method, we implement BIC score minimization using a grid search and two prefixed tuning parameters ${\lambda }_n=0.5\log (p_n)$ and $\log (p_n)$ for comparative purposes. We report the bias ( ${\Vert \widehat{\mathit{\beta }}-{\mathit{\beta }}_0\Vert }_2$), number of FP, FN, and BIC score ( $-2l_n(\widehat{\mathit{\beta }})+\log (n){\sum }_{j}I({\widehat{\mathit{\beta }}}_j\ne 0)$) in Table 2.
在此仿真中，我们模拟了 n=200,000、p _n =20,000 和 q _n =80 的 sHDMSS 事件发生时间数据集。事件时间由指数风险模型生成，基线风险 h ₀ （t）=1，回归系数 ${\mathit{\beta }}_0=({\mathbf{0}\mathbf{.}\mathbf{7}}_{10},{\mathbf{0}\mathbf{.}\mathbf{5}}_{10},{\mathbf{0}\mathbf{.}\mathbf{8}}_{10},{\mathbf{1}}_{10},{\mathbf{-}\mathbf{0}\mathbf{.}\mathbf{7}}_{10},{\mathbf{-}\mathbf{0}\mathbf{.}\mathbf{5}}_{10},{\mathbf{-}\mathbf{0}\mathbf{.}\mathbf{8}}_{10},{\mathbf{-}\mathbf{1}}_{10},{\mathbf{0}}_{p_n-80})$ ，删失率为 95%。模拟每个主题的协变量，平均分配 2% 的非零值。用于存储此密集设计矩阵的内存量将需要超过 16 GB，这超出了标准硬件上大多数统计软件包的功能容量。为了克服这一困难，我们以坐标列表的方式有效地存储信息，并使用 Cyclops 软件包 12， 38 将 BAR 的大量 Cox 回归 （mBAR） 与 LASSO 的大量稀疏 Cox 回归 （mCox-LASSO） 进行比较，据我们所知，它是当今最快的软件，它利用 sHDMSS 事件时间数据的稀疏性进行高效计算，并提供> 10 倍的加速 12 比其竞争对手（如 CoxNet 7） 快 7 倍和 FastCox。47 对于LASSO来说，使用CV（mCox-LASSO（CV））结合比经典网格搜索方法更有效的非凸优化技术，以及使用经典网格搜索方法实现的BIC分数最小化（mCox-LASSO（BIC））来找到调谐参数的最优值。对于mBAR方法，我们使用网格搜索和两个前缀调优参数 ${\lambda }_n=0.5\log (p_n)$ 来实现BIC分数最小化，并 $\log (p_n)$ 用于比较目的。我们在表 2 中报告了偏差 （ ${\Vert \widehat{\mathit{\beta }}-{\mathit{\beta }}_0\Vert }_2$ ）、FP、FN 和 BIC 评分 （ $-2l_n(\widehat{\mathit{\beta }})+\log (n){\sum }_{j}I({\widehat{\mathit{\beta }}}_j\ne 0)$ ） 的数量。

Table 2. (Sparse high-dimensional and massive sample size) Estimation and variable selection results for massive Cox regression with broken adaptive ridge (mBAR) and LASSO penalty (mCox-LASSO12

) for a simulated sHDMSS data set with n=200 000, p_n=20 000, and q_n=80. (Bias = ${\Vert \widehat{\mathit{\beta }}-{\mathit{\beta }}_0\Vert }_2$; FP= number of false positives; FN = number of false negatives)
表 2.（稀疏高维和海量样本量）对于n=200 000、p _n =20 000和q _n =80的模拟sHDMSS数据集，具有破碎自适应脊（mBAR）和LASSO惩罚（mCox-LASSO 12）的大规模Cox回归的估计和变量选择结果。（偏差 = ${\Vert \widehat{\mathit{\beta }}-{\mathit{\beta }}_0\Vert }_2$ ;FP=误报数;FN = 假阴性数）

Method 方法	Bias 偏见	FP	FN	BIC score BIC 评分
mBAR ( ${\lambda }_n=0.5\log (p_n)$) mBAR （ ${\lambda }_n=0.5\log (p_n)$ ）	1.19	0	3	83 313.02
mBAR ( ${\lambda }_n=\log (p_n)$) mBAR （ ${\lambda }_n=\log (p_n)$ ）	2.02	0	10	83 573.96
mBAR (BIC) mBAR （BIC）	0.97	5	0	83 266.47
mCox-LASSO (BIC) mCox-LASSO （BIC）	2.93	12	0	84 479.47
mCox-LASSO (CV) mCox-LASSO （CV）	2.12	963	0	93 770.58

• Abbreviations: BIC, Bayesian information criterion; CV, cross-validation.
  缩写：BIC，贝叶斯信息准则;简历，交叉验证。

We observe that both mCox-LASSO methods have retained all 80 true nonzero coefficients together with a moderate to large number of noise variables (12 for BIC and 967 for CV). In contrast, mBAR (BIC) chooses a sparser model selecting all 80 nonzero coefficients and 5 noise variables. As expected, mBAR (BIC) is less biased (0.82) than mCox-LASSO (2.49 for BIC and 2.02 for CV) and has a much lower BIC score when compared to both mCox-LASSO methods. We also notice that mBAR with the two prefixed λ_n tends to underestimate the true model, ie, fixing ${\lambda }_n=\log (p_n)$ results in estimating a model that is too sparse, whereas ${\lambda }_n=0.5\log (p_n)$ produces a model that is closer to the oracle model.
我们观察到，两种 mCox-LASSO方法都保留了所有 80 个真正的非零系数以及中等到大量的噪声变量（BIC 为 12 个，CV 为 967 个）。相比之下，mBAR （BIC） 选择稀疏模型，选择所有 80 个非零系数和 5 个噪声变量。正如预期的那样，mBAR （BIC） 的偏差 （0.82） 低于 mCox-LASSO（BIC 为 2.49，CV 为 2.02），并且与两种 mCox-LASSO 方法相比，BIC 分数要低得多。我们还注意到，带有两个前缀 λ 的 mBAR _n 往往会低估真实模型，即修复 ${\lambda }_n=\log (p_n)$ 导致估计模型过于稀疏，而 ${\lambda }_n=0.5\log (p_n)$ 生成更接近预言机模型的模型。

We further examined the solution paths of mCox-LASSO and mBAR in Figure 2. The vertical solid and dashed lines in the mCox-LASSO solution path plot (Figure 2A) represent the estimates at the optimal tuning parameter obtained via CV and BIC minimization, respectively. We can see that the mCox-LASSO solution path changes rapidly as its tuning parameter varies and shows severe bias. In contrast, the mBAR solution path plot (Figure 2B) with respect to λ_n changes very slowly where the vertical line represents the estimates at the optimal tuning parameter selected by BIC minimization and selects a model with estimates that are less biased than mCox-LASSO (see Table 2). Furthermore, the optimal value of λ_n that minimizes the BIC score for mBAR roughly corresponds to $0.3\log (p_n)$. Since our empirical results suggest that the optimal value for λ_n generally lies within some constant of $\log (p_n)$, we recommend that a coarse grid search within $c\log (p_n)$ where c∈(0,1] can be used. This is further corroborated by additional simulations in the Supplementary Material (Tables S1 to S3).
我们进一步检查了图2中mCox-LASSO和mBAR的求解路径。mCox-LASSO解决方案路径图（图2A）中的垂直实线和虚线分别表示通过CV和BIC最小化获得的最佳调谐参数的估计值。我们可以看到，mCox-LASSO求解路径随着其调谐参数的变化而迅速变化，并显示出严重的偏差。相比之下，mBAR解路径图（图2B）相对于λ _n 的变化非常缓慢，其中垂直线表示通过BIC最小化选择的最佳调谐参数的估计值，并选择估计值比mCox-LASSO偏差小的模型（见表2）。此外，最小化 mBAR 的 BIC 分数的 λ _n 的最优值大致对应于 $0.3\log (p_n)$ 。由于我们的实证结果表明，λ 的最优值 _n 通常位于 的 $\log (p_n)$ 某个常数内，我们建议在 的范围内进行粗略网格搜索 $c\log (p_n)$ ，其中 c∈（0,1）。补充材料中的其他模拟进一步证实了这一点（表S1至S3）。

For the mBAR method, we also made a solution path plot with respect to ξ_n, while fixing ${\lambda }_n=\log (p_n)$ in Figure 2C. It shows that the mBAR estimates are very stable over a large range of ξ_n, affirming our observation in Section 3.1 with small scale data that mBAR is generally insensitive ξ_n.
对于mBAR方法，我们还绘制了关于ξ的求解路径图 _n ，同时 ${\lambda }_n=\log (p_n)$ 固定在图2C中。它表明mBAR估计值在ξ的大范围内非常稳定 _n ，证实了我们在第3.1节中用小尺度数据观察到的mBAR通常是不敏感的ξ _n 。

4 PEDIATRIC TRAUMA MORTALITY
4 小儿创伤死亡率

For an application of mBAR regression in the sHDMSS setting, we consider a subset of the NTDB, a trauma database maintained by the American College of Surgeons.1 This data set was previously analyzed by Mittal et al12 as an example for efficient massive Cox regression with mCox-LASSO and ridge regression to sHDMSS data. The data set includes 210 555 patient records of injured children under 15 that were collected over 5 years (2006 to 2010). Each patient record includes 125 952 binary covariates, which indicate the presence or absence of an attribute (ICD9 Codes, AIS codes, etc) as well as the two-way interactions. The outcome of interest is mortality after time of injury. The data is extremely sparse, with less than 1% of the covariates being nonzero and has a censoring rate of 98%. We randomly split the data into training and test sets of 168 000 and 42 555, respectively. The mortality rate of both sets were approximately equal to the combined rate. Similar to Section 3.3, we were unable to load the training set (n=168 000, p_n=125 000) into other popular oracle procedures due to the memory requirements needed to support a dense design matrix of that size and compare mBAR to mCox-LASSO. The BIC-score minimization over a penalization path of 10 tuning parameters was used to select the final model for both mBAR (fixing ${\xi }_n=\log (p_n)$) and mCox-LASSO. In addition, we perform mCox-LASSO using CV and mBAR with fixed tuning parameters ${\lambda }_n=0.5\log (p_n)$ and $\log (p_n)$. The BIC score based on the training data is used to compare selection performance between models and discriminatory performance is measured using Harrell's c-statistic48, 49 based on the test data.
对于mBAR回归在sHDMSS环境中的应用，我们考虑NTDB的一个子集，NTDB是由美国外科医师学会维护的创伤数据库。1 Mittal 等人 12 之前分析了该数据集，作为使用 mCox-LASSO 进行有效大规模 Cox 回归和 sHDMSS 数据的岭回归的示例。该数据集包括210 555份15岁以下受伤儿童的病历，这些病历是在5年内（2006年至2010年）收集的。每个患者记录包括 125 952 个二元协变量，这些协变量表示属性（ICD9 代码、AIS 代码等）的存在与否以及双向交互。感兴趣的结果是受伤后的死亡。数据非常稀疏，只有不到 1% 的协变量为非零，删失率为 98%。我们将数据随机分为训练集和测试集，分别为168 000和42 555。两组死亡率大致等于合并死亡率。与第 3.3 节类似，我们无法将训练集 （n=168 000， p _n =125 000） 加载到其他流行的预言机程序中，因为支持该大小的密集设计矩阵并将 mBAR 与 mCox-LASSO 进行比较所需的内存要求。在10个调谐参数的惩罚路径上，BIC分数最小化用于选择mBAR（固定 ${\xi }_n=\log (p_n)$ ）和mCox-LASSO的最终模型。此外，我们使用 CV 和 mBAR 执行 mCox-LASSO，并具有固定的调谐参数 ${\lambda }_n=0.5\log (p_n)$ 和 $\log (p_n)$ .基于训练数据的BIC分数用于比较模型之间的选择性能，并使用基于测试数据的Harrell's c统计量48,49来测量判别性能。

Table 3 summarizes the findings for our example, which reflect what we observe in Section 3.3. Massive Cox's regression for BAR, using BIC minimization, selects fewer covariates than both mCox-LASSO methods. Both model selection and discriminatory performance are similar to slightly superior for mBAR (BIC) over both mCox-LASSO methods. Again, mBAR with prefixed λ_n selects far fewer covariates than mBAR (BIC); however, the overall high c-index for both methods suggest that the strong predictors for pediatric trauma are still retained in the model. In terms of runtime, mBAR (BIC) is more time consuming than LASSO (BIC) as expected, but BAR with a prefixed tuning parameter value can help to reduce the runtime with a comparable prediction performance.
表 3 总结了我们示例的结果，这些结果反映了我们在第 3.3 节中观察到的内容。使用 BIC 最小化的 BAR 的 Massive Cox 回归选择的协变量比两种 mCox-LASSO方法都少。mBAR （BIC） 的模型选择和判别性能都略优于 mCox-LASSO 方法。同样，前缀为 λ 的 mBAR _n 选择的协变量比 mBAR （BIC） 少得多;然而，两种方法的总体高 C 指数表明，模型中仍保留了儿科创伤的强预测因子。在运行时方面，mBAR （BIC） 比 LASSO （BIC） 更耗时，但具有前缀调优参数值的 BAR 可以帮助减少运行时间，并具有相当的预测性能。

5 DISCUSSION 5 讨论

We have extended the BAR methodology to Cox's model as a new sparse Cox regression method and rigorously established that it is selection consistent, oracle for parameter estimation, stable, and has a grouping property for highly correlated covariates. We illustrate through empirical studies that the BAR estimator has satisfactory performance for variable selection and parameter estimation. We have also extended the application of BAR to the sHDMSS domain by taking advantage of the fact that the BAR algorithm allows us to easily adapt existing high performance algorithms and software for massive ℓ₂-penalized Cox regression.12
我们已将 BAR 方法扩展到 Cox 模型，作为一种新的稀疏 Cox 回归方法，并严格确定它是选择一致的、用于参数估计的预言机、稳定的，并且具有高度相关的协变量的分组属性。通过实证研究证明，BAR估计器在变量选择和参数估计方面具有令人满意的性能。我们还利用BAR算法允许我们轻松调整现有的高性能算法和软件，以实现大规模l ₂ 惩罚Cox回归这一事实，将BAR的应用扩展到sHDMSS域。12

Our surrogate ℓ₀-based BAR method and theory can be easily extended to a surrogate ℓ_d-based BAR method for any d∈[0,1], by replacing ${({\widehat{\beta }}_j^{(k-1)})}^2$ with ${|{\widehat{\beta }}_j^{(k-1)}|}^{2-d}$ in 4. We have observed empirically that, as d increases toward 1, the resulting estimator becomes less sparse, and the average number of FP as well as estimation bias tend to increase, especially for larger p_n, while the average number of FN tends to decrease. In practice, d can be used as a resolution tuning parameter.
我们基于代理 l ₀ 的 BAR 方法和理论可以很容易地扩展到任何 d∈[0,1] 的代理 l 的 BAR 方法，只需将 ${|{\widehat{\beta }}_j^{(k-1)}|}^{2-d}$ in 4 替换 ${({\widehat{\beta }}_j^{(k-1)})}^2$ 为代理 l _d 。我们根据经验观察到，当 d 向 1 增加时，得到的估计器变得不那么稀疏，FP 的平均数和估计偏差趋于增加，尤其是对于较大的 p _n ，而 FN 的平均数趋于减少。在实践中，d 可以用作分辨率调谐参数。

Our theoretical and empirical results have established the BAR method as a valid and viable tool for variable selection and parameter estimation under the p_n<n setting although p_n is allowed to diverge with n. Theoretical properties of the BAR estimator for the high-dimensional setting (p_n≫n) remain to be investigated. Furthermore, as pointed out by a referee, although BAR is selection consistent and oracle, it is subject to the same postselection inference issues as other variable selection methods.50, 51 Lastly, although iteratively performing reweighted ℓ₂-penalizations allows BAR to enjoy the best of ℓ₀- and ℓ₂-penalized regressions and to readily adopt an existing efficient implementation of ℓ₂-penalization for sHDMSS data, its iterative nature does present another layer of computational complexity. While this added layer of computational complexity is not a practical concern for moderate size data, it can considerably increase the runtime in a large data setting when both n and p are large. As illustrated in our real data example, trying a prefixed tuning parameter value based on the extended BIC ${\lambda }_n=c\log (p_n)$ can reduce the runtime of BAR with reasonably good performance. To further improve its computational efficiency, we are currently developing some modified BAR algorithms including a cyclic coordinatewise BAR algorithm, which will have comparable computational complexity and runtime to other popular variable methods such as LASSO. This line of further developments is beyond the scope of this paper and will be fully studied in a sequel paper.
我们的理论和实证结果表明，尽管p _n 允许与n发散，但BAR方法在p _n ≫n）的BAR估计器的理论性质仍有待研究。此外，正如裁判所指出的，尽管BAR是选择一致的，并且是预言机，但它与其他变量选择方法一样，会受到相同的选择后推断问题的影响。50， 51 最后，尽管迭代执行重加权 l ₂ 惩罚化允许 BAR 享受 l ₀ 和 l ₂ 惩罚回归的最佳效果，并随时采用现有的 l ₂ 惩罚化 sHDMSS 数据的有效实现，但其迭代性质确实带来了另一层计算复杂性。虽然对于中等大小的数据来说，这种增加的计算复杂性层并不是一个实际问题，但当 n 和 p 都很大时，它可以大大增加大型数据设置中的运行时间。如我们的真实数据示例所示，尝试基于扩展 BIC ${\lambda }_n=c\log (p_n)$ 的前缀调优参数值可以减少 BAR 的运行时间，并具有相当好的性能。为了进一步提高其计算效率，我们目前正在开发一些改进的BAR算法，包括循环坐标BAR算法，该算法的计算复杂度和运行时间将与其他流行的变量方法（如LASSO）相当。这一系列的进一步发展超出了本文的范围，将在后续论文中进行全面研究。

CONFLICT OF INTEREST 利益冲突

The authors declare no potential conflict of interest.
作者声明没有潜在的利益冲突。

DATA AVAILABILITY STATEMENT
数据可用性声明

The data that support the findings of this study are available from the American College of Surgeons.1 Restrictions apply to the availability of these data.
支持这项研究结果的数据可从美国外科医师学会获得。1 这些数据的可用性受到限制。