数学杂志  2022, Vol. 42 Issue (5): 445-460   PDF    
扩展功能
加入收藏夹
复制引文信息
加入引用管理器
Email Alert
RSS
本文作者相关文章
张佳倩
邓立凤
丁洁丽
比例风险模型下病例队列设计中两种加权估计方法及其应用
张佳倩1, 邓立凤2, 丁洁丽1    
1. 武汉大学数学与统计学院, 湖北 武汉 430072;
2. 山东科技大学数学与系统科学学院, 山东 青岛 266590
摘要:对于大型队列研究或观察型研究, 基于生存数据的病例队列设计是一种能有效节约成本和提高效率的抽样机制. 这种抽样设计仅对一个随机抽取的子队列以及子队列之外所有经历了感兴趣事件的病例个体进行关键协变量的测量, 具有显著的成本效益. 本文研究如何应用比例风险模型拟合病例队列研究数据. 探讨逆概率加权和与时间相关加权这两种基于加权估计方程的统计推断方法和其渐近性质等理论结果. 通过一系列的统计模拟研究展示了病例队列设计的优良性以及相较于传统简单随机抽样设计的高效性. 进一步, 应用这两种推断方法分析了两个实际数据, 展示了其在实际中的应用价值和前景.
关键词病例队列设计    比例风险模型    逆概率加权法    与时间相关权法    
NFERENCE AND APPLICATION OF CASE-COHORT DESIGN UNDER THE PROPORTIONAL HAZARDS MODEL
ZHANG Jia-qian1, DENG Li-feng2, DING Jie-li1    
1. School of Mathematics and Statistics, Wuhan University, Wuhan, Hubei 430072, China;
2. College of Mathematics and Systems Science, Shandong University of Science and Technology, Qingdao, Shandong, 266590, China
Abstract: A case-cohort design is a cost-effective sampling scheme in large cohort studies. The key idea of such a design is to assemble the measurements of expensive covariates only on a subset of the entire cohort and all the subjects outside the subcohort that experience the event of interest. In this paper, we study the inference methods for case-cohort data under the Cox model. We consider two weighted estimating equation approaches, the inverse-probability and time-varying weighted methods. The asymptotic theories are established. A series of simulation studies are conducted to assess the finite-sample performance of the proposed methods and exhibit the superiority and efficiency of the case-cohort design. Some real data examples are analyzed to illustrate the application of the proposed methods.
Keywords: Case-Cohort design     proportional hazards model     inverse probability weight     time varying weight    
1 引言

流行病学、生物医学和遗传学等领域的观察型研究对了解人类疾病的影响因素至关重要. 在许多这样的研究中, 往往涉及对大量研究个体的长期追踪和随访. 因而, 这些研究的主要的预算和成本通常来自于昂贵协变量的采集与观测. 对于预算有限的大型观察型研究, 若采用传统的简单随机抽样可能会导致试验过于昂贵而超出预算. 因此, 发展和研究能节约成本和提高效率的抽样机制具有非常重要的理论意义和应用价值.

对于带有删失的失效时间数据, 病例队列设计(case-cohort design) 是应用最为广泛的有偏抽样机制之一. 其关键思想是: 首先, 从全队列中随机地抽取一个子队列(subcohort). 全队列中所有发生或者经历了感兴趣的事件(例如: 疾病, 死亡等) 的个体称为病例(case). 然后, 子队列和子队列之外所有的病例一起组成病例队列样本(case-cohort sample). 最后, 仅对病例队列样本中的个体进行昂贵协变量的采集和观测. 这种病例队列设计对于需要观测昂贵协变量的大型队列研究或观察型研究是具有成本效益的. 尤其是感兴趣的事件是稀发事件时, 此种抽样设计非常高效.

在Prentice (1986)[1]首次提出病例队列设计之后, 对于病例队列设计和相关统计方法的研究有了大量和广泛的工作, 这些工作主要沿着两条研究思路, 一条思路是基于似然函数的方法(例如: Prentice, 1986[1]; SelfPrentice, 1988[2]; ChenLo, 1999[3]; Kang et al, 2016[4]; Liu et al, 2018[5]等等), 一条思路是估计方程的方法(例如: LinYing, 1993[6]; KulichLin, 2000[7]; Sun et al, 2004[8]; CaiZeng, 2004[9], 2007[10]; KulichLin, 2004[11]; Kong et al, 2004[12]; LuTsiatis, 2006[13]; BreslowWellner, 2007[14]; Kang et al, 2013[15]; SteingrimssonStrawderman, 2017[16]等等). 在估计方程方法的研究中, 对基于不完全观测数据的逆概率加权思想的探讨尤其广泛, 例如: 加性风险模型(KulichLin, 2000[7]); 加速失效模型(KongCai, 2009[17]); 半参数转移模型(Kong et al, 2004[12]; LuTsiatis, 2006[13]). 在对带有多元失效时间的病例队列研究中, 一些加权估计方程应用了一种与时间相关的权函数并发展了相应的统计推断方法(KangCai, 2009[18]; Kang et al, 2013[15]; Yan et al, 2017[19]).

病例队列设计作为一种具有成本效益的有偏抽样机制, 其在各个领域的应用越来越广泛. 本文主要探讨病例队列研究中的逆概率加权和与时间相关加权这两种思想下的估计方程方法及其应用. 首先, 我们探讨如何用比例风险模型拟合病例队列数据, 考虑基于逆概率加权和与时间相关加权两种思想下的统计推断方法, 综述其渐近性质等理论结果. 然后, 我们重点研究上述两种方法在实际中的应用问题, 编写可实现这两种方法的操作性强的统计软件应用程序. 通过一系列的模拟研究来评估上述两种方法在有限样本下的表现, 展示了病例队列设计相较于简单随机抽样方法的优良性和有效性. 最后,我们将这两种方法应用于分析两个实际数据, 展示了该方法在实际中的应用价值与前景. 本文主要结构为: 第2节介绍比例风险模型和病例队列设计的抽样机制, 第3节探讨逆概率加权法和与时间相关的加权法这两种估计方法, 第4节为研究模拟计算, 第5节研究实际数据分析与应用.

2 模型与抽样设计

$ \widetilde{T} $表示潜在失效时间, $ C $表示删失时间. 记$ T=\min(\widetilde{T}, C) $为观测时间, $ \Delta=I(\widetilde{T}\leq C) $为右删失示性变量, 其中$ I(\cdot) $为示性函数. 记$ p $维协变量为$ {\mathit{\boldsymbol{Z}}}=(Z_1, Z_2, ..., Z_p)^T $. 假设研究队列包含$ N $个独立个体, $ \{T_i, \Delta_i, {\mathit{\boldsymbol{Z}}}_i: i=1, ..., N\} $为其观测数据. 记$ \tau $为实验终止时间. 我们考虑如下比例风险模型(Cox, 1972[20]), 即: 给定协变量$ {\mathit{\boldsymbol{Z}}} $时, 失效时间$ \widetilde{T} $的风险率函数为:

$ \begin{equation} \lambda(t|{\mathit{\boldsymbol{Z}}})=\lambda_0(t)\exp({\mathit{\boldsymbol{Z}}}^T{\mathit{\boldsymbol{\beta}}}), \end{equation} $ (2.1)

其中$ \lambda_0(t) $是未知的基准风险函数, $ {\mathit{\boldsymbol{\beta}}}=(\beta_1, ..., \beta_p)^T $为待估的$ p $维回归系数.

在队列研究中, 当协变量对每个个体均进行观测时, 对$ {\mathit{\boldsymbol{\beta}}} $的推断问题广泛地应用如下的偏似然函数法(Cox, 1972[20]; AndersenGill, 1982[21]):

$ \begin{equation} \ell_F({\mathit{\boldsymbol{\beta}}})=\sum\limits_{i=1}^n \Delta_i \left[ {{\mathit{\boldsymbol{Z}}}_i^T{\mathit{\boldsymbol{\beta}}}-\log\left\{{\sum\limits_{l=1}^N Y_l(T_i)\exp({\mathit{\boldsymbol{Z}}}_l^T{\mathit{\boldsymbol{\beta}}})}\right\}} \right]. \end{equation} $ (2.2)

$ {\mathit{\boldsymbol{\beta}}} $的估计可定义为: $ \mathit{\boldsymbol{{\widehat{\beta}}}}_F=\underset{{\mathit{\boldsymbol{\beta}}}}{\arg\max} \ell_F({\mathit{\boldsymbol{\beta}}}). $

在病例队列设计下, 首先, 通过简单随机抽样方式从全队列中抽取一个样本容量为$ {\widetilde{n}} $的子队列. 子队列中的个体和所有子队列之外的所有病例组成病例队列样本, 记其容量为$ n $. 然后, 仅对病例队列样本进行协变量的测量. 特别地, 记$ \xi_i $为子队列示性变量, 即: 如果第$ i $个个体被选入子队列, 则$ \xi_i=1 $; 否则, $ \xi_i=0 $. 假设每个个体入选子队列的概率为$ \text{P}\left( {\xi_i=1} \right) = \widetilde{\alpha} = {\widetilde{n}}/N. $因此, 病例队列设计下的观测数据结构可总结为:

$ \begin{align} \left\{ \begin{array}{ll} \text{当} \xi_i=1 \text{或} \Delta_i=1 \text{时}: \quad &(T_i, \, \Delta_i, \, Z_i);\\ \text{当} \xi_i=0 \text{且} \Delta_i=0 \text{时}: \quad &(T_i, \, \Delta_i).\ \end{array} \right. \end{align} $ (2.3)

由于病例队列设计中, 协变量不是对每个个体进行观测, 因此, 传统推断方法不再适用, 需要发展新的分析方法. 以下我们探讨两种基于加权估计方程的推断方法.

3 估计方法
3.1 逆概率加权法

受HorvitzThompson (1951)[22]思想的启发, 即: 对于不完全观测的数据进行逆概率加权(inverse-probability weight, 以下简称"IPW"), 在病例队列设计下, 对模型(2.1) 中$ {\mathit{\boldsymbol{\beta}}} $的推断可建立如下加权估计方程:

$ \begin{equation} U_w({\mathit{\boldsymbol{\beta}}})=\sum\limits_{i=1}^N \Delta_i \left[ {Z_i-\frac{S_w^{(1)}({\mathit{\boldsymbol{\beta}}}, T_i)}{S_w^{(0)}({\mathit{\boldsymbol{\beta}}}, T_i)}} \right], \end{equation} $ (3.1)

其中$ S_w^{(d)}({\mathit{\boldsymbol{\beta}}}, t)=\frac{1}{N}\sum\limits_{i=1}^N w_l Y_l(t)e^{{\mathit{\boldsymbol{Z}}}'_l{\mathit{\boldsymbol{\beta}}}}{\mathit{\boldsymbol{Z}}}_l^{\otimes d}, $这里$ a^{\otimes 0}=1, a^{\otimes 1}=a, a^{\otimes 2}=aa' $, 其中$ a $是一个向量. 权重$ w_i $定义为:

$ \begin{equation} w_i=\Delta_i+\frac{(1-\Delta_i)\xi_i}{\widetilde{\alpha}}, \qquad i=1, ..., N. \end{equation} $ (3.2)

注意到, 对于病例组中的个体, 无论是否选入子队列, 其权重均为$ 1 $; 而对于子队列中的对照组个体, 其权重为$ \widetilde{\alpha}^{-1} $. 这种逆概率加权的思想由KalbfleischLawless (1988)[23]首次应用于分析生存数据. Borgan et al (2000)[24]在分层病例队列研究中引入了相似的逆概率权. 对于两阶段抽样设计, KulichLin (2004)[11]也提出了类似的加权估计方法. 在比例风险模型下, KangCai (2009)[18]为多类型疾病的病例队列设计发展了基于多元失效时间的逆概率加权推断方法.

估计方程(3.1) 的解定义为逆概率加权估计, 记为$ \widehat{{\mathit{\boldsymbol{\beta}}}}_{IPW} $. 以下, 我们综述$ \widehat{{\mathit{\boldsymbol{\beta}}}}_{IPW} $的渐近性质. 记$ {\mathit{\boldsymbol{\beta}}}_0 $$ {\mathit{\boldsymbol{\beta}}} $的真值. 定义$ N_i(t) = \Delta_i I(T_i\leq t) $$ Y_i(t) = I(T_i\geq t) $分别为计数过程和风险过程. 定义$ M_i(t)= N_i(t)-\int_{0}^{t} Y_i(s)\lambda_0(s)e^{{\mathit{\boldsymbol{Z}}}'_i{\mathit{\boldsymbol{\beta}}}}ds $. 假设如下正则条件成立.

(A1) $ {\mathit{\boldsymbol{\beta}}} $的参数空间$ \mathcal{B} $是紧的.

(A2) 协变量$ {\mathit{\boldsymbol{Z}}} $的取值空间$ \mathcal{Z} $是紧的.

(A3) 给定$ Z_i $时, $ T_i $$ C_i $相互独立. $ \xi_i $$ (T_i, \delta_i, Z_i) $相互独立.

(A4) 存在某个$ \alpha\in (0, 1) $, 使得$ \widetilde{\alpha}=\widetilde{n}/N \to \alpha $.

(A5) $ \int_0^{\tau} \lambda_0 (t)dt < \infty $ (A6) $ \text{P}(I(T_1\geq t)=1 $, 对任意$ t\in[0, \tau]) > 0 $.

(A7) 对$ d=0, 1, 2 $, $ \sup\limits_{{\mathit{\boldsymbol{\beta}}}\in\mathcal{B}, t\in [0, \tau]}\left\Vert{S_w^{(d)}({\mathit{\boldsymbol{\beta}}}, t)-s^{(d)}({\mathit{\boldsymbol{\beta}}}, t)}\right\Vert\stackrel{P}{\to}0. $其中$ s^{(d)}({\mathit{\boldsymbol{\beta}}}, t)=E\left[ {Y_1(t)e^{{\mathit{\boldsymbol{Z}}}'_1{\mathit{\boldsymbol{\beta}}}}{\mathit{\boldsymbol{Z}}}_1^{\otimes d}} \right] $. $ s^{(d)}({\mathit{\boldsymbol{\beta}}}, t) $$ t\in[0, \tau] $上一致的关于$ {\mathit{\boldsymbol{\beta}}} $连续, 并且$ s^{(0)}({\mathit{\boldsymbol{\beta}}}, t) $下方有界且大于零.

(A8) 矩阵

$ \Sigma({\mathit{\boldsymbol{\beta}}}_0)=\int_0^{\tau}\left[ {\frac{s^{(2)}({\mathit{\boldsymbol{\beta}}}, t)}{s^{(0)}({\mathit{\boldsymbol{\beta}}}_0, t)} -\left\{{\frac{s^{(1)}({\mathit{\boldsymbol{\beta}}}, t)}{s^{(0)}({\mathit{\boldsymbol{\beta}}}_0, t)}}\right\}^{\otimes 2}} \right]s^{(0)}({\mathit{\boldsymbol{\beta}}}_0, t)\lambda_0(t)dt, $

是正定的.

下面定理给出了$ \widehat{{\mathit{\boldsymbol{\beta}}}}_{IPW} $的渐近性质.

定理1 ($ \widehat{{\mathit{\boldsymbol{\beta}}}}_{IPW} $的渐近性质)   在正则条件(A1)-(A8)下, $ \widehat{{\mathit{\boldsymbol{\beta}}}}_{IPW} $依概率收敛于$ {\mathit{\boldsymbol{\beta}}}_0 $, 即:

$ \begin{equation*} \widehat{{\mathit{\boldsymbol{\beta}}}} {\stackrel{\text{P}}{\longrightarrow}} {\mathit{\boldsymbol{\beta}}}_0. \end{equation*} $

进一步,

$ \begin{equation} \sqrt{N}(\widehat{{\mathit{\boldsymbol{\beta}}}}_{IPW}-{\mathit{\boldsymbol{\beta}}}_0){\stackrel{d}{\longrightarrow}} N\left( {0, \Omega({\mathit{\boldsymbol{\beta}}}_0)} \right), \end{equation} $ (3.3)

其中: 渐近方差矩阵$ \Omega({\mathit{\boldsymbol{\beta}}}_0)= \Sigma^{-1}({\mathit{\boldsymbol{\beta}}}_0)\left\{{\Sigma_1({\mathit{\boldsymbol{\beta}}}_0)+\Sigma_2({\mathit{\boldsymbol{\beta}}}_0)}\right\}\Sigma^{-1}({\mathit{\boldsymbol{\beta}}}_0) $, $ \Sigma_1({\mathit{\boldsymbol{\beta}}}_0) =E\left( {G_1({\mathit{\boldsymbol{\beta}}}_0)^{\otimes 2}} \right), $ $ \Sigma_2({\mathit{\boldsymbol{\beta}}}_0) =\frac{1-\alpha}{\alpha}E\left[ {(1-\Delta_1)G_1({\mathit{\boldsymbol{\beta}}}_0)^{\otimes 2}} \right], $这里

$ \begin{equation*} G_1({\mathit{\boldsymbol{\beta}}}) =\int_0^{\tau}\left[ {Z_1-\frac{s^{(1)}({\mathit{\boldsymbol{\beta}}}, t)}{s^{(0)}({\mathit{\boldsymbol{\beta}}}, t)}} \right]dM_1(t). \end{equation*} $

此定理基于与KangCai (2009)[18]相似的讨论和证明思路即可证明.

进一步, 为了实际应用中的计算问题, 我们为渐近方差$ \Omega({\mathit{\boldsymbol{\beta}}}_0) $提出如下一种相合估计. 定义:

$ \begin{align*} \widehat{\Sigma}({\mathit{\boldsymbol{\beta}}}) &= \frac{1}{N}\sum\limits_{i=1}^N \Delta_i\left[ {\frac{S_w^{(2)}({\mathit{\boldsymbol{\beta}}}, T_i)}{S_w^{(0)}({\mathit{\boldsymbol{\beta}}}, T_i)} -\left( {\frac{S_w^{(1)}({\mathit{\boldsymbol{\beta}}}, T_i)}{S_w^{(0)}({\mathit{\boldsymbol{\beta}}}, T_i)}} \right)^{\otimes 2}} \right], \\ \widehat{G}_i({\mathit{\boldsymbol{\beta}}}) &= \Delta_i\left[ {Z_i-\frac{S_w^{(1)}({\mathit{\boldsymbol{\beta}}}, T_i)}{S_w^{(0)}({\mathit{\boldsymbol{\beta}}}, T_i)}} \right] -\frac{1}{N}\sum\limits_{l=1}^N\frac{w_l\Delta_lI(T_l\leq T_i)e^{{\mathit{\boldsymbol{Z}}}'_i{\mathit{\boldsymbol{\beta}}}}}{S_w^{(0)}({\mathit{\boldsymbol{\beta}}}, T_l)} \left[ {{\mathit{\boldsymbol{Z}}}_i-\frac{S_w^{(1)}({\mathit{\boldsymbol{\beta}}}, T_l)}{S_w^{(0)}({\mathit{\boldsymbol{\beta}}}, T_l)}} \right], \\ \widehat{\Sigma}_1({\mathit{\boldsymbol{\beta}}}) &= \frac{1}{N}\sum\limits_{i=1}^N w_i\widehat{G}_i({\mathit{\boldsymbol{\beta}}})^{\otimes 2}, \\ \widehat{\Sigma}_2({\mathit{\boldsymbol{\beta}}}) &= \frac{1-\widetilde{\alpha}}{\widetilde{\alpha}}\frac{1}{N}\sum\limits_{i=1}^N w_i(1-\Delta_i)\widehat{G}_i({\mathit{\boldsymbol{\beta}}})^{\otimes 2}. \end{align*} $

$ \Omega({\mathit{\boldsymbol{\beta}}}_0) $可由$ \widehat{\Omega}({\widehat{\mathit{\boldsymbol{\beta}}}_{IPW}}) = \widehat{\Sigma}^{-1}({\widehat{\mathit{\boldsymbol{\beta}}}_{IPW}})\left( {\widehat{\Sigma}_1({\widehat{\mathit{\boldsymbol{\beta}}}_{IPW}})+\widehat{\Sigma}_2({\widehat{\mathit{\boldsymbol{\beta}}}_{IPW}}} \right)\widehat{\Sigma}^{-1}({\widehat{\mathit{\boldsymbol{\beta}}}_{IPW}}) $估计.

3.2 时间相关权法

在逆概率加权思想基础上, Barlow (1994)[25]考虑了一种随时间变化的加权思想(time-varying weight, 以下简称"TVW“). Borgan et al (2000)[24]和KulichLin (2004)[11]将这种与时间相关的权分别应用于分层病例队列设计和两阶段抽样设计. 还有一些研究工作发展和建立了分析多类型疾病的病例队列数据的与时间相关加权(TVW) 法(KangCai, 2009[18]; KangCai, 2013[15]; Yan et al, 2017[19]). 受到上述研究的启发, 我们探讨如下与时间相关的权函数:

$ \begin{equation} w_i(t)=\Delta_i+\frac{(1-\Delta_i)\xi_i}{\widehat{\alpha}(t)}, \end{equation} $ (3.4)

其中

$ \begin{equation*} \widehat{\alpha}(t)=\frac{\sum_{i=1}^N(1-\Delta_i)\xi_iY_i(t)}{\sum_{i=1}^N(1-\Delta_i)Y_i(t)}. \end{equation*} $

在TVW的加权思想下, 所有病例组的个体其权重均为$ 1 $; 而被选入子队列中的对照组个体权重为$ \widehat{\alpha}(t) $. 注意到$ \widehat{\alpha}(t) $的定义为: 在时刻$ t $, 入选子队列并仍在风险集中的对照组个体数与所有仍在风险集中的对照组个体个数之比. 故而, $ \widehat{\alpha}(t) $为真实的子队列入选概率$ \widetilde{\alpha} $的一个估计. 基于上述TVW权, 在病例队列下, 对模型(2.1) 中$ {\mathit{\boldsymbol{\beta}}} $的估计可建立如下加权得分方程,

$ \begin{equation} \widetilde{U}_W({\mathit{\boldsymbol{\beta}}})=\sum\limits_{i=1}^N \Delta_i \left[ {Z_i-\frac{{\widetilde{S}_w}^{(1)}({\mathit{\boldsymbol{\beta}}}, T_i)}{{\widetilde{S}_w}^{(0)}({\mathit{\boldsymbol{\beta}}}, T_i)}} \right], \end{equation} $ (3.5)

其中$ {\widetilde{S}_w}^{(d)}({\mathit{\boldsymbol{\beta}}}, t)=\frac{1}{N}\sum\limits_{i=1}^N w_l(t) Y_l(t)e^{{\mathit{\boldsymbol{Z}}}'_l{\mathit{\boldsymbol{\beta}}}}{\mathit{\boldsymbol{Z}}}_l^{\otimes d}. $由方程(3.5) 求解得到的$ {\mathit{\boldsymbol{\beta}}} $的估计, 我们记为$ \widehat{{\mathit{\boldsymbol{\beta}}}}_{TVW} $. 以下, 我们研究和综述$ \widehat{{\mathit{\boldsymbol{\beta}}}}_{TVW} $的渐近理论. 我们修改条件(A7) 为如下的条件(A7'):

(A7') 对$ d=0, 1, 2 $,

$ \sup\limits_{{\mathit{\boldsymbol{\beta}}}\in\mathcal{B}, t\in [0, \tau]}\left\Vert{\widetilde{S}_w^{(d)}({\mathit{\boldsymbol{\beta}}}, t)-s^{(d)}({\mathit{\boldsymbol{\beta}}}, t)}\right\Vert\stackrel{P}{\to}0. $

其中$ s^{(d)}({\mathit{\boldsymbol{\beta}}}, t)=E\left[ {Y_l(t)e^{{\mathit{\boldsymbol{Z}}}'_l{\mathit{\boldsymbol{\beta}}}}{\mathit{\boldsymbol{Z}}}_l^{\otimes d}} \right] $. $ s^{(d)}({\mathit{\boldsymbol{\beta}}}, t) $$ t\in[0, \tau] $上一致的关于$ {\mathit{\boldsymbol{\beta}}} \in \mathcal{B} $连续, 并且$ s^{(0)}({\mathit{\boldsymbol{\beta}}}, t) $下方有界且大于零.

下面定理给出了$ \widehat{{\mathit{\boldsymbol{\beta}}}}_{TVW} $的渐近性质.

定理2 ($ \widehat{{\mathit{\boldsymbol{\beta}}}}_{TVW} $渐近性质)   在正则条件(A1)-(A6), (A7') 和(A8) 下, $ \widehat{{\mathit{\boldsymbol{\beta}}}}_{TVW} $依概率收敛于$ {\mathit{\boldsymbol{\beta}}}_0 $. 进一步,

$ \begin{equation} \sqrt{N}(\widehat{{\mathit{\boldsymbol{\beta}}}}_{TVW}-{\mathit{\boldsymbol{\beta}}}_0){\stackrel{d}{\longrightarrow}} N\left( {0, \widetilde{\Omega}({\mathit{\boldsymbol{\beta}}}_0)} \right), \end{equation} $ (3.6)

其中渐近方差矩阵$ \widetilde{\Omega}({\mathit{\boldsymbol{\beta}}}_0) =\Sigma^ {-1}({\mathit{\boldsymbol{\beta}}}_0)\left\{{\Sigma_1({\mathit{\boldsymbol{\beta}}}_0)+\widetilde{\Sigma}_2({\mathit{\boldsymbol{\beta}}}_0)}\right\}\Sigma^{-1}({\mathit{\boldsymbol{\beta}}}_0) $, $ \widetilde{\Sigma}_2({\mathit{\boldsymbol{\beta}}})=\frac{1-\alpha}{\alpha}V_1({\mathit{\boldsymbol{\beta}}}_0), $此处

$ \begin{equation*} V_1({\mathit{\boldsymbol{\beta}}}) =Var\left( {(1-\Delta_1)\int_0^{\tau}\left[ {\widetilde{R}_1({\mathit{\boldsymbol{\beta}}}, t) -\frac{Y_1(t)E\left\{{(1-\Delta_1)\widetilde{R}_1({\mathit{\boldsymbol{\beta}}}, t)}\right\}}{E\left\{{(1-\Delta_1)Y_1(t)}\right\}}} \right]d\Lambda_0(t)} \right), \end{equation*} $

其中$ \widetilde{R}_i({\mathit{\boldsymbol{\beta}}}, t)=Y_i(t)\left[ {Z_i(t)-\frac{s^{(1)}({\mathit{\boldsymbol{\beta}}}, t)}{s^{(0)}({\mathit{\boldsymbol{\beta}}}, t)}} \right]e^{{\mathit{\boldsymbol{\beta}}}'Z_i(t)}. $

此定理证明可参考KangCai (2013)[15].

进一步, 我们提出渐近方差矩阵$ \widetilde{\Omega}({\mathit{\boldsymbol{\beta}}}_0) $的一种估计方法. 我们采用非参数自助法(Hjort, 1985[26]; Efron & Tibshirani, 1993[27]; Burr, 1994[28]) 来估计$ \widehat{{\mathit{\boldsymbol{\beta}}}}_{TVW} $的渐近方差. 其基本思想是通过对观测数据的重复抽样建立经验分布函数. 当相关分布未知时, 非参数自助法是一种应用广泛的计算估计值方差或标准差的数值计算方法. 具体地, 记$ Y_{obs}=\{{\mathit{\boldsymbol{X}}}_1, ..., {\mathit{\boldsymbol{X}}}_N\} $为一组观测数据, 其中$ {\mathit{\boldsymbol{X}}}_i=(T_i, \Delta_i, {\mathit{\boldsymbol{Z}}}_i) $. 从观测数据$ Y_{obs}=\{{\mathit{\boldsymbol{X}}}_1, ..., {\mathit{\boldsymbol{X}}}_N\} $中有放回地随机抽取一组新样本$ Y_{obs}^\ast=\{{\mathit{\boldsymbol{X}}}_1^\ast, ..., {\mathit{\boldsymbol{X}}}_N^\ast\} $, 其中每个$ {\mathit{\boldsymbol{X}}}_i^\ast $的入样概率均为$ 1/N $. 基于抽取的Bootstrap样本$ Y_{obs}^\ast $, 我们计算加权估计$ \widehat{{\mathit{\boldsymbol{\beta}}}}_{TVW}^\ast $. 独立地重复上述过程$ B $次, 可得到$ B $个Bootstrap估计$ \{\widehat{{\mathit{\boldsymbol{\beta}}}}_{TVW}^\ast(b)\}^B_{b=1} $, 其中$ \widehat{{\mathit{\boldsymbol{\beta}}}}_{TVW}^\ast(b) = (\widehat{\beta}_1^\ast(b), ..., \widehat{\beta}_p^\ast(b))^T $. 因此, $ \widehat{{\mathit{\boldsymbol{\beta}}}}_{TVW} $的方差的第$ k $个分量可由如下的样本方差估计:

$ \begin{equation*} \widehat{ \text{Var}}(\widehat{\beta}_k)=\frac{1}{B-1}\sum\limits_{b=1}^B\left[ {\widehat{\beta}_k^\ast(b) -\frac{1}{B}\sum\limits_{b=1}^B\widehat{\beta}_k^\ast(b)} \right]^2, \quad k=1, ..., p. \end{equation*} $

进一步地, 可以得到$ \widehat{\beta}_k $的Wald型Bootstrap置信区间. 当$ \{\widehat{\beta}^\ast_k(b)\}_{b=1}^B $近似正态分布时, $ \widehat{\beta}_k $$ (1-\alpha)100\% $置信区间为$ [\widehat{\beta}_k - z_{\alpha/2}\cdot \widehat{ \text{se}}(\widehat{\beta}_k), \widehat{\beta}_k + z_{\alpha/2}\cdot \widehat{ \text{se}}(\widehat{\beta}_k)], $其中$ z_{\alpha} $表示标准正态分布的上$ \alpha $分位数, $ \widehat{ \text{se}}(\widehat{\beta}_k)=\sqrt{\widehat{ \text{Var}}(\widehat{\beta}_k)} $.

4 模拟研究

本节通过一系列模拟研究来评估这两种方法在有限样本量下的表现, 展示其在实际中的可行性和应用性. 我们考虑如下比例风险模型:

$ \begin{equation*} \lambda(t|Z_{1}, Z_{2})=\lambda_{0}(t)\exp\{\beta_{1}Z_{1}+\beta_{2}Z_{2}\}. \end{equation*} $

我们设定参数真值为$ \beta_1=0, 0.693 $$ \beta_2=0, -0.5 $. 协变量$ Z_1 $从成功概率为$ 0.5 $的Bernoulli分布中生成, $ Z_2 $从标准正态分布中生成. 设定基准风险函数$ \lambda_{0}(t) $$ 1 $$ 2t $. 由此, 失效时间$ \widetilde{T} $可以分别由风险率为$ \exp\{\beta_{1}Z_{1}+\beta_{2}Z_{2}\} $的指数分布和形状参数为$ 2 $、尺度参数为$ \left[ {\exp\{\beta_{1}Z_{1}+\beta_{2}Z_{2}} \right]^{-1/2} $的 Weibull分布中随机生成. 删失时间$ C $从均匀分布$ U[0, \, c] $中生成, 通过挑选$ c $的不同取值从而产生相应的删失率, 分别为$ \rho=80\% $$ 90\% $. 对于病例队列设计, 设定子队列的样本量分别为$ \widetilde{n}=200 $$ 300 $, 从样本总量为$ N=1000 $的全队列中随机地抽取获得.

我们关注以下几个问题: 第一, 使用病例队列设计替代简单随机抽样设计能提高多少效率? 第二, 逆概率加权法和与时间相关权法这两种推断方法在有限样本量下表现如何? 第三, 两种推断方法的估计效率是否不同? 因此, 我们比较了以下几种方法:

Full: 全队列下基于比例风险模型的偏似然估计法($ {\widehat{\mathit{\boldsymbol{\beta}}}_{Full}} $);

Naive: 基于与病例队列样本相同样本量的简单随机抽样下的偏似然估计法($ {\widehat{\mathit{\boldsymbol{\beta}}}_{Naive}} $);

IPW: 病例队列设计下的逆概率加权法($ {\widehat{\mathit{\boldsymbol{\beta}}}_{IPW}} $);

TVW: 病例队列设计下与时间相关权法($ {\widehat{\mathit{\boldsymbol{\beta}}}_{TVW}} $).

在每种参数设定下, 比较上述四种方法参数估计值相较于真值的偏差(Bias), 估计值的样本标准差(SD), 标准差估计值的均值(SE), $ 95\% $的正态区间覆盖率(CP), 以及估计的相对效率(RE). 所有结果均基于$ 1000 $次的模拟获得. 模拟结果请见表 1-4.

表 1 参数$ \beta_1 $$ \beta_2 $的模拟结果, 其中子队列样本量为$ \widetilde{n}=200 $, 基准风险函数$ \lambda_0(t)=1 $.

表 2 参数$ \beta_1 $$ \beta_2 $的模拟结果, 其中子队列样本量为$ \widetilde{n}=300 $, 基准风险函数$ \lambda_0(t)=1 $.

表 3 参数$ \beta_1 $$ \beta_2 $的模拟结果, 其中子队列样本量为$ \widetilde{n}=200 $, 基准风险函数$ \lambda_0(t)=2t $.

表 4 参数$ \beta_1 $$ \beta_2 $的模拟结果, 其中子队列样本量为$ \widetilde{n}=300 $, 基准风险函数$ \lambda_0(t)=2t $.

在所有考虑的情况下, 关于$ \beta_1 $$ \beta_2 $的四种估计都是无偏的, 标准误差估计的均值(SEs) 很好地估计了估计值的样本标准差(SDs), 置信区间覆盖率(CPs) 均约为$ 95\% $. 从结果来看, $ {\widehat{\mathit{\boldsymbol{\beta}}}_{IPW}} $$ {\widehat{\mathit{\boldsymbol{\beta}}}_{TVW}} $的表现相似, $ {\widehat{\mathit{\boldsymbol{\beta}}}_{IPW}} $$ {\widehat{\mathit{\boldsymbol{\beta}}}_{TVW}} $的效率稍高一点. 这说明IPW方法由于依据的是真实抽样概率$ \widetilde{\alpha} $, 其估计要比依据真实抽样概率估计$ \widehat{\alpha}(t) $的TVW方法更有效一些, 因为$ \widehat{\alpha}(t) $仅通过删失个体信息得到. 特别地, 对于高删失率的情况, 两者的差距要更小一些. 从相对效率(REs) 来看, $ {\widehat{\mathit{\boldsymbol{\beta}}}_{IPW}} $$ {\widehat{\mathit{\boldsymbol{\beta}}}_{TVW}} $均比$ {\widehat{\mathit{\boldsymbol{\beta}}}_{Naive}} $更有效. 这说明相较于简单随机抽样设计, 病例队列设计能有效地提高估计的效率. 例如: 关于$ \beta_1 $的估计, 当$ \widetilde{n}=200, \lambda_0(t)=1, \beta_1=0, \beta_2=-0.5, \rho=90\% $时, 相较于$ {\widehat{\mathit{\boldsymbol{\beta}}}_{Full}} $, $ {\widehat{\mathit{\boldsymbol{\beta}}}_{Naive}} $, $ {\widehat{\mathit{\boldsymbol{\beta}}}_{IPW}} $$ {\widehat{\mathit{\boldsymbol{\beta}}}_{TVW}} $的相对效率分别为$ 0.30 $, $ 0.61 $$ 0.70 $, 即: 一方面, 相同样本量下, 病例队列设计的效率约为简单随机抽样设计下的$ 2 $倍; 另一方面, 病例队列设计仅用了约全队列$ 30\% $的样本量, 却达到了$ 60\% $多的效率. 再比如, 关于$ \beta_2 $的估计, 当$ \widetilde{n}=300, \lambda_0(t)=2t, \beta_1=0.693, \beta_2=-0.5, \rho=90\% $时, 相较于$ {\widehat{\mathit{\boldsymbol{\beta}}}_{Full}} $, $ {\widehat{\mathit{\boldsymbol{\beta}}}_{Naive}} $, $ {\widehat{\mathit{\boldsymbol{\beta}}}_{IPW}} $$ {\widehat{\mathit{\boldsymbol{\beta}}}_{TVW}} $的相对效率分别为$ 0.35 $, $ 0.68 $$ 0.67 $, 即病例队列设计仅用约全队列$ 35\% $的样本量实现了相较于简单随机抽样设计的两倍效率.

总体来说, 删失率较高$ (\rho=90\% v.s. 80\%) $或子队列样本量更大$ (\widetilde{n}=200\% v.s. 300) $时, $ {\widehat{\mathit{\boldsymbol{\beta}}}_{IPW}} $$ {\widehat{\mathit{\boldsymbol{\beta}}}_{TVW}} $的效率更高. 另外, 当参数真值更接近于$ 0 $$ (\beta_1=0.693 v.s. \beta_1=0 $, $ \beta_2=-0.5 v.s. \beta_2=0 $), 两种估计方法的效率也有小幅度的提高, 其中$ {\widehat{\mathit{\boldsymbol{\beta}}}_{TVW}} $相较于$ {\widehat{\mathit{\boldsymbol{\beta}}}_{IPW}} $提高得尤为明显. 例如: 当$ \widetilde{n}=200, \rho=90\%, \lambda_0(t)=1 $时, 相较于全样本估计, $ {\widehat{\mathit{\boldsymbol{\beta}}}_{TVW}} $$ \beta_1=0 $$ \beta_2=0 $时的效率比在$ \beta_1=0.693 $$ \beta_2=-0.5 $时的效率提高了约$ 20\% $.

5 实际数据分析

本节通过分析一个肾母细胞瘤的病例队列数据和一个员工离职管理的队列数据, 来展示病例队列设计以及其所探讨的两种加权估计推断方法在实际中的应用.

5.1 肾母细胞瘤

肾母细胞瘤是一种多发于幼儿的罕见肾癌. 美国国家肾母细胞瘤研究组织(NWTSG) 进行了一系列研究来探索患儿的肿瘤组织学类型与发病时间之间的联系. 我们研究的数据来自于NWTSG的第三次、第四次临床试验的数据(BreslowChatterjee, 1999[29]; D'Angio et al, 1989[30]; Green et al, 1998[31]), 其中包含$ 4028 $个幼儿的记录.

患儿的肿瘤组织学类型是最重要的暴露因素. 地方医疗机构的病理学家在患儿最初接受治疗时, 对其肿瘤组织学类型进行初步评估. 然后, 来自NWTSG病理中心的有经验的病理学家对其组织学类型进行确定性的评估. 后者的评估结果更加的准确, 但其费用也更加的昂贵且其过程更加的耗时. 因此, BreslowChatterjee (1999)[29]为此研究提出了一个病例队列设计方案. 具体地, 从全队列$ 4028 $个患儿中随机地抽取了$ 668 $个患儿作为子队列. 病例队列样本由子队列以及子队列之外所有经历了发病或者死亡的患儿组成. 在NWTSG病理中心中仅对病例队列样本中的患儿肿瘤确定其组织类型. 我们应用本文研究的两种加权估计方程方法来分析这样一个病例队列数据.

我们感兴趣的因变量是患儿的发病时间, 而观测到的发病时间带有右删失, 其删失率约为$ 85.8\% $. 我们考虑三个潜在影响因素. 肿瘤组织学类型(Histype) 分为两种类型: 一种是肿瘤由被称为“组织结构不良型” (Histype $ =1 $) 的罕见细胞类型组成, 另一种是肿瘤由“组织结构良好型” (Histype $ =0 $) 的细胞组成. 疾病分期(Stage)分为以下四种类型: 肿瘤广泛分布于肾脏并被完全切除(Stage $ =1 $), 肿瘤超出肾脏但完全切除(Stage $ =2 $), 腹腔中有残余肿瘤或淋巴结中有肿瘤(Stage $ =3 $), 癌细胞转移到肺或者肝脏(Stage $ =4 $). 确诊年龄(Age) 以月为单位, 我们对其进行了中心标准化处理. 我们对数据中所有考虑的影响因素进行了描述性统计分析, 画出了组织学类型条形图, 疾病分期饼图以及确诊年龄直方图, 见图 1, 23.

图 1 组织学类型条形图

图 2 疾病分期条形图

图 3 确诊年龄直方图

我们的目的是评估肿瘤组织学类型, 疾病分期和确诊年龄对患儿发病时间(T) 的影响. 我们考虑如下比例风险模型:

$ \begin{equation*} \lambda(t|Z) = \lambda_{0}(t) \exp\{\beta_{1}\mathrm{Histype}+\beta_{2}\mathrm{Stage}+\beta_{3}\mathrm{Age}\}. \end{equation*} $

在此模型框架下, 我们应用逆概率加权法(IPW), 与时间相关权法(TVW), 以及基于子队列样本的偏似然估计法(SRS)分析了肾母细胞瘤研究数据. 分析结果总结在表 5中.

表 5 肾母细胞瘤研究数据分析结果

总体来说, 三种方法分析得到的估计结果是一致的. 结果表明, 肿瘤组织学类型对发病时间的影响是显著的. 组织结构不良型的患儿较之组织结构良好型的患儿发病或者死亡的风险更高, $ e^{1.3489}=3.9 $倍(SRS), $ e^{1.3336}=3.8 $倍(IPW) 或$ e^{1.3428}=3.8 $倍(TVW). 结果也显示, 疾病分期程度越高的患儿发病和死亡的风险越高, 确诊年龄越小的患儿发病和死亡的风险越高. 在三种估计中, 病例队列设计下的IPW法和TVW法比简单随机抽样下的SRS法更有效. 考虑到肿瘤组织学类型的评估十分昂贵且耗时, 采用病例队列抽样设计能提高研究效率和节约成本.

5.2 员工离职管理

在现代商业社会, 研究员工自身的特点与企业文化是否相适应对于稳定企业员工团队至关重要. 我们研究的数据来源于国内某大型商业银行人力资源部门[32], 共有$ 1300 $个样本, 每个样本对应于该银行当年的一位销售人员. 我们关注的因变量是这些员工的工作年限$ T $ (单位: 月), 在观测时间内有的员工已离职而有的员工仍在继续工作, 因此观测数据带有右删失, 删失率约为$ 22.69\% $. 为了展示病例队列抽样设计, 我们从全队列$ 1300 $个样本中随机地抽取了$ 300 $个作为子队列, 子队列以及子队列之外所有已离职的员工数据组成病例队列样本. 我们考虑比例风险模型, 研究可能与员工工作年限相关的三个潜在影响因素: 一个是员工户籍$ X_1 $, $ X_1=0 $表示本地员工, $ X_1=1 $表示异地员工; 一个是员工性别$ X_2 $, $ X_2=0 $表示性别为男, $ X_2=1 $表示性别为女; 还有一个是员工年龄$ X_3 $, 我们对其作了中心标准化处理. 对数据中上述所有协变量进行了描述性统计分析, 结果见表 6.

表 6 员工离职数据描述性统计分析表.

我们应用逆概率加权法(IPW), 与时间相关权法(TVW) 以及基于相同样本量的简单随机抽样下的偏似然估计法(Naive) 分析了员工离职数据. 分析结果总结在表 7中.

表 7 员工离职研究数据分析结果.

总体来说, 三种分析方法得到的估计结果是一致的. 结果表明, 员工的户籍、性别和年龄特征都与其工作年限显著相关. 具体地, 异地员工相比于本地员工离职的风险更高, 分别约为$ e^{0.2662}=1.3 $倍(Naive), $ e^{0.2730}=1.3 $倍(IPW) 或$ e^{0.2979}=1.3 $倍(TVW); 男员工相比于女员工离职的风险更高, 分别约为$ e^{0.1995}=1.2 $倍(Naive), $ e^{0.2868}=1.3 $倍(IPW) 或$ e^{0.3718}=1.5 $倍(TVW); 年龄越小的员工离职的风险越高. 在三种估计中, 病例队列设计下的IPW法和TVW法比简单随机抽样下的Naive法更有效.

6 总结

病例队列设计作为应用最为广泛的有偏抽样机制之一, 能有效节约研究成本和提高研究效率. 本文研究了病例队列设计在比例风险模型下的应用, 探讨了逆概率加权和与时间相关加权这两种思想下的估计方程方法, 并综述了其渐近理论. 然后, 我们重点研究上述两种方法在实际中的应用问题, 为上述两种方法编写了一套操作性强的计算方法与可实现这两种方法的应用程序. 模拟研究结果表明在病例队列设计下, 上述两种方法在有限样本下均表现优异, 其估计效率均明显高于传统简单随机抽样设计. 最后, 实际数据分析结果也展示了估计方法在实际中的应用价值.

为了进一步提高估计效率, 未来的工作包括研究基于最优权的推断方法. 相关的判别准则有赤池信息准则(AIC) (Akaike, 1973[33]); 贝叶斯信息准则(BIC) (Schwarz, 1978[34]) 或广义交叉验证(CGV) (Craven & Wanba, 1978[35]). 当感兴趣事件发生率较低时, 病例队列设计的效果更好. 对于删失率较低的生存数据, 相关的研究包括广义病例队列设计(Cai & Zeng, 2007[10]) 和基于因变量抽样设计(Ding et al, 2014[36]; Yu et al, 2016[37]) 等等. 这些也将是我们未来的主要工作之一.

参考文献
[1] Prentice R L. A case-cohort design for epidemiologic cohort studies and disease prevention trials[J]. Biometrika, 1986, 73: 1–11. DOI:10.1093/biomet/73.1.1
[2] Self S G, Prentice R L. Asymptotic distribution theory and efficiency results for case-cohort studies[J]. The Annals of Statistics, 1988, 16: 64–81.
[3] Chen Kani, Lo Shaw-Wha. Case-cohort and case-control analysis with Cox's model[J]. Biometrika, 1999, 86: 755–764. DOI:10.1093/biomet/86.4.755
[4] Kang Suhyun, Lu Wenbin, Liu Mengling. Efficient estimation for accelerated failure time model under case-cohort and nested case-control sampling[J]. Biometrics, 2016, 73: 114–123.
[5] Liu Dandan, Cai Tianxi, Lok A, Yingye. Nonparametric maximum likelihood estimators of time-dependent accuracy measures for survival outcome under two-stage sampling designs[J]. Journal of the American Statistical Association, 2018, 113: 522.
[6] Lin D Y, Ying Z. Cox regression with incomplete covariate measurements[J]. Journal of the American Statistical Association, 1993, 88: 1341–1349. DOI:10.1080/01621459.1993.10476416
[7] Kulich M, Lin D Y. Additive hazards regression for case-cohort studies[J]. Biometrika, 2000, 87: 73–87. DOI:10.1093/biomet/87.1.73
[8] Sun J, Sun L, Flournoy N. Addictive hazards model for competing risks analysis of the case-cohort design[J]. Communication in Statistics - Theory and Method, 2004, 33: 351–366. DOI:10.1081/STA-120028378
[9] Cai Jianwen, Zeng Donglin. Sample size/power calculation for case-cohort studies[J]. Biometric, 2004, 60: 1015–1024. DOI:10.1111/j.0006-341X.2004.00257.x
[10] Cai Jianwen, Zeng Donglin. Power calculation for case-cohort studies with nonrare events[J]. Biometrics, 2007, 63: 1288–1295. DOI:10.1111/j.1541-0420.2007.00838.x
[11] Kulich M, Lin D Y. Improving the efficiency of relative-risk estimation in case-cohort studies[J]. Journal of the American Statistical Association, 2004, 99: 832–844. DOI:10.1198/016214504000000584
[12] Kong Lan, Cai Jianwen, Sen Pranab K. Weighted estimating equations for semiparametric transformation models with censored data from a case-cohort design[J]. Biometrika, 2004, 91: 305–319. DOI:10.1093/biomet/91.2.305
[13] Lu Wenbin, Tsiatis Anastasios A. Semiparametric transformation models for the case-cohort study[J]. Biometrika, 2006, 93: 207–214. DOI:10.1093/biomet/93.1.207
[14] Breslow N E, Wellner J A. Weighted likelihood for semiparametric models and two-phase stratified samples, with application to cox regression[J]. Scandinavian Journal of Statistics, 2007, 34: 86–102. DOI:10.1111/j.1467-9469.2006.00523.x
[15] Kang Sangwook, Cai Jianwen, Chambless L. Marginal additive hazards model for case-cohort studies with multiple disease outcomes: an application to the Atherosclerosis Risk in Communities (ARIC) study[J]. Biostatistics, 2013, 14: 28–41. DOI:10.1093/biostatistics/kxs025
[16] Steingrimsson J A, Strawderman R L. Estimation in the semiparametric accelerated failure time model with missing covariates: improving efficiency through augmentation[J]. Journal of the American Statistical Association, 2017, 112: 519.
[17] Kong Lan, Cai Jianwen. Case-cohort analysis with accelerated failure time model[J]. Biometrics, 2009, 65: 135–142. DOI:10.1111/j.1541-0420.2008.01055.x
[18] Kang Sangwook, Cai Jianwen. Marginal hazard model for case-cohort studies with multiple disease outcomes[J]. Biometrika, 2009, 96: 887–901. DOI:10.1093/biomet/asp059
[19] Yan Ying, Zhou Haibo, Cai Jianwen. Improving efficiency of parameter estimation in case-cohort studies with multivariate failure time data[J]. Biometrics, 2017, 73: 1042–1052. DOI:10.1111/biom.12657
[20] Cox D R. Regression models and life-tables[J]. Journal of the Royal Statistical Society. Series B (Methodological), 1972, 34: 187–220. DOI:10.1111/j.2517-6161.1972.tb00899.x
[21] Andersen P K, Gill R D. Cox's regression model for counting processes: a Large sample study[J]. The Annals of Statistics, 1982, 10: 1100–1120.
[22] Horvitz D G, Thompson D J. A generalization of sampling without replacement from a finite universe[J]. Journal of the American Statistical Association, 1951, 47: 663–685.
[23] Kalbfleisch J D, Lawless J F. Likelihood analysis of multi-state models for disease incidence and mortality[J]. Statistics in Medicine, 1988, 7: 147–160.
[24] Borgan O, Langholz B, Samuelsen S, Goldstein L, Pogoda J. Exposure stratified case-cohort designs[J]. Lifetime Data Analysis, 2000, 6: 86–102.
[25] Barlow W. Robust variance estimation for the case-cohort design[J]. Biometrics, 1994, 50: 1064–72. DOI:10.2307/2533444
[26] Hjort N. Bootstrapping Cox's Regression Model[R]. California: Stanford University, Dept. of Statistics, 1985.
[27] Efron B, Tibshirani R. An Introduction to the Bootstrap[J]. Journal of Educational and Behavioral Statistics, 1993, 22: 245–245.
[28] Burr D. A Comparison of Certain Bootstrap Confidence Intervals in the Cox Model[J]. Journal of the American Statistical Association, 1994, 89: 1290–1302. DOI:10.1080/01621459.1994.10476869
[29] Breslow N E, Chatterjee N. Design and analysis of two-phase studies with binary outcome applied to Wilms tumor prognosis[J]. Journal of the Royal Statistical Society, Series C, 1999, 48: 457–468. DOI:10.1111/1467-9876.00165
[30] D'Angio G J, Breslow N, Bechwith J B, et al. Treatment of Wilms' tumor, results of the third national Wilms' tumor study[J]. Cancer, 1989, 64: 349–360. DOI:10.1002/1097-0142(19890715)64:2<349::AID-CNCR2820640202>3.0.CO;2-Q
[31] Green D M, Breslow N E, Bechwith J B, et al. Comparison between single-dose and divided-dose administration of dactinomycin and doxorubicin for patients with Wilms' tumor: a report from the national Wilms' tumor study group[J]. Journal of Clinical Oncology, 1998, 16: 237–245. DOI:10.1200/JCO.1998.16.1.237
[32] 王汉生. 商务数据分析与应用[M]. 北京: 中国人民大学出版社, 2011.
[33] Akaike H. Information theory and an extension of the maximum likelihood principle[J]. In B. N. Petrov and F. Caski (Eds.), Proceedings of the Second International Symposium on Information Theory,, 1973: 267–281.
[34] Schwarz G. Estimating the dimension of a model[J]. The Annals of Statistics, 1978, 6: 461–464.
[35] Craven P, Wahba G. Smoothing noisy data with spline functions[J]. Numerische mathematik, 1978, 31: 377–403.
[36] Ding J, Zhou H, Liu Y, Cai J, Longnecker M P. Estimating effect of environmental contaminants on women's subfecundity for the MoBa study data with an outcome-dependent sampling scheme[J]. Biostatistics, 2014, 15: 636–650.
[37] Yu J, Liu Y, Cai J, Sandler D P, Zhou H. Outcome-dependent sampling design and inference for Cox's proportional hazards Model[J]. Journal of Statistical Planning and Inference, 2016, 178: 24–36. DOI:10.1016/j.jspi.2016.05.001