生存数据应用在许多学科领域, 如生物医学、流行病学、公共卫生学、可靠性工程学、经济学和保险精算学等领域. 生存分析主要研究某特定事件的发生时间(例如: 死亡时间、疾病发生时间、系统失效时间、索赔时间等等) 与重要影响因素和协变量之间的关联. 对于生存数据的研究, 比例风险模型是应用最为广泛的统计半参数模型之一(Cox, 1972[1]; Andersen & Gill, 1982[2]; Lin, 1994[3]; Chen & Little, 1999[4]; Kalbfleisch & Prentice, 2002[5] 等等). 作为比例风险模型的一个重要替代, 加性风险模型假设基准风险函数与协变量效应之间具有一个加性结构. 在很多实际应用中, 加性风险模型往往能更好地拟合数据(Lin & Ying, 1994[3]; Mckeague & Sasieni, 1994[6]). 而当加性风险模型和比例风险模型均能较好地拟合数据时, 加性风险模型的回归参数更容易解释其实际意义(Lin & Ying, 1994[3]; Zeng & Cai, 2010[7]).
在许多大型队列研究中, 往往涉及对大量研究个体的长期追溯和随访. 当重要影响因素或协变量的采集非常昂贵时, 采用传统的简单随机抽样可能会导致实验过于昂贵而超过预算. 因此, 发展和研究能节约成本和提高效率的抽样机制具有非常重要的意义. 对于带有删失的生存数据, 病例队列设计(Case-Cohort design) 是应用最为广泛的有偏抽样机制之一. 其主要机制是: 首先, 从全队列中随机地抽取一个子队列(subcohort), 全队列中所有发生了感兴趣事件的个体称为病例(case). 然后, 子队列和子队列之外所有的病例组成病例队列样本. 最后, 仅对病例队列样本中的个体进行昂贵协变量的采集和观测. 对病例队列设计相关统计方法的研究已有大量和广泛的工作, 比例风险模型(Prentice, 1986[8]; Self & Prentice, 1988[9]; Chen & Lo, 1999[10]; Lin & Ying, 1993[11]; Kulich & Lin, 2004[12]), 加性风险模型(Kulich & Lin, 2000[13]; Sun et al, 2004[14]), 加速失效模型(Kong et al, 2004[15]; Lu & Tsiatis, 2006[16]). 近来, 基于多元失效时间的病例队列设计的研究越来越广泛(Kong & Cai, 2009[17]; Kang et al, 2013[18]; Kim et al, 2018[19]; Maitra et al, 2020[20]).
本文主要探讨病例队列设计中加性风险模型下参数的统计推断方法和应用. 首先, 我们探讨如何应用加性风险模型拟合由病例队列设计获取的生存数据, 考虑参数的一种加权估计方法并综述其渐近理论. 然后, 重点研究这种病例队列设计下的分析方法在实际中的应用问题. 一方面, 我们编写了可实现这种统计分析方法的R软件应用程序. 通过模拟研究展示了这种方法在有限样本量下的优良表现, 并评估了病例队列设计相较于简单随机抽样设计的有效性. 另一方面, 我们应用该方法分析了一个实际的乳腺癌数据, 展示了其成本效益与应用价值. 乳腺癌是女性最常见的恶性肿瘤之一, 全世界每年约有120万女性患上乳腺癌, 约有50万女性患者死亡. 中国癌症年发病数为160万, 现症病人200多万. 乳腺癌占女性恶性肿瘤发病率第一位, 每年约有4万女性死于乳腺癌[21]. 本文将基于来源于美国国家癌症研究所(National Cancer Institute)[22] 的乳腺癌数据, 探索影响乳腺癌患者生存时间的影响因素. 首先, 我们基于全队列数据应用加性风险回归方法分析数据, 探索出了一些对乳腺癌患者生存时间有显著性影响的因素. 进一步, 应用病例队列设计抽取样本, 并基于病例队列样本进行加性风险回归分析. 结果表明, 病例队列设计仅用了较小的样本量就达到了与全队列研究几乎一致的结果. 当协变量的测量非常昂贵时, 病例队列设计可提高研究效率和节约成本.
本文的主要结构为: 第2节介绍病例队列设计下加性风险模型参数的加权估计方法并综述其渐近理论, 第3节为模拟计算, 第4节为乳腺癌数据的生存分析.
假设一个大型队列研究包含$ N $个独立的研究个体. 我们用$ \widetilde{T}_i $表示第$ i $个个体的潜在失效时间, $ C_i $表示第$ i $个个体的删失时间或随访时间$ (i = 1, \cdots , N) $. 记观测时间为$ T_i = \min(\widetilde{T}_i, C_i) $, 右删失示性变量为$ \Delta_i = I(\widetilde{T}_i \leq C_i) $. 用$ Y_i(t) = I(T_i\geq t) $表示风险过程, $ N_i(t) = \Delta_iI(T_i\leq t) $表示计数过程, 其中$ I(\cdot) $表示示性函数. 记$ Z_i(t) $为第$ i $个个体在$ t $时刻处的$ p $维协变量. 记$ \tau $为事件终止时间.
我们考虑如下加性风险模型: 给定协变量$ Z_i(t) $时, 失效时间$ \widetilde{T}_i $的风险函数有如下形式:
其中$ \lambda_{0}(t) $为形式未知的基准风险函数, $ \beta $为感兴趣的$ p $维待估参数. 记基准累积风险函数为$ \Lambda_0(t) = \int_{0}^{t}\lambda_0(s)ds $. 在队列研究中, 当对每个个体的协变量进行观测时, 对$ \beta $和$ \Lambda_0(t) $的推断广泛地使用如下估计方程方法(Lin & Ying, 1994)[3]:
记$ \widehat{\beta} $和$ \widehat{\Lambda}_{0}(t) $为上述估计方程组的解, 有如下形式:
和
其中$ a^{\otimes2} = aa' $, 以及
在研究的病例队列设计下, 首先, 通过简单随机抽样方式从全队列中抽取一个样本容量为$ n_0 $的子队列. 子队列中的个体和子队列之外的所有病例个体组成病例队列样本, 记其样本量为$ n $. 然后, 仅对病例队列样本中的个体观测其协变量. 具体来说, 记$ \xi_i $为子队列示性变量, 即: $ \xi_i = 1 $表示第$ i $个个体被选入子队列, $ \xi_i = 0 $表示第$ i $个个体未被选入子队列. 假设每个个体入选子队列的概率$ P(\xi_i = 1) = \widetilde{\alpha} = n_0/N $. 因此, 病例队列设计下, 当$ \xi_i = 1 $或$ \Delta_i = 1 $时, 观测数据为$ (T_i, \Delta_i, Z_i) $; 当$ \xi_i = 0 $且$ \Delta_i = 0 $时, 观测数据为$ (T_i, \Delta_i) $.
由于病例队列设计中, 协变量不是对每一个个体都进行了观测, 因此(2) 中的估计方程方法不再适用, 需要提出新的推断方法. 受Horvitz & Thompson(1951)[23] 逆概率加权思想和Liang & Ziger(1986)[24] 广义估计方程思想的启发. 在病例队列的设计下, 对加性风险模型(1) 中的$ \beta $和$ \Lambda_0(t) $的推断可建立如下加权估计方程:
其中
上述加权估计方程有如下显式解:
这种逆概率加权的思想由Kalbfleisch & Lawless(1988)[25] 首次应用于生存分析数据. 基于多类型疾病的病例队列设计, Kang et al(2013)[18] 在加性风险模型下为多元失效时间发展了逆概率加权推断方法. 以下我们讨论和综述了上述$ \widehat{\beta}_w $的渐近性质. 设$ \beta_0 $为$ \beta $的真值. 假设如下正则条件成立:
(C1) $ \beta $的参数空间$ \mathbb{B} $是紧集, $ Z $的取值空间$ \mathbb{Z} $是紧集.
(C2) 给定$ Z_i $时, $ \widetilde{T_i} $与$ C_i $相互独立, $ \xi_i $与$ (T_i, \Delta_i, Z_i) $相互独立.
(C3) $ P(T_{i}{\geq}\tau)>0 $且$ \Lambda_{0}(t)<\infty $.
(C4) 存在某个$ \alpha \in (0, 1) $, 使得当$ N\rightarrow \infty $时, $ \widetilde{\alpha} = n_0/N \rightarrow \alpha $.
(C5) 矩阵
是非奇异矩阵, 其中
定理 1 ($ \widehat{\beta}_w $的渐近性质) \label{thm1} 在正则条件(C1)-(C5) 下, $ \widehat{\beta}_w $是$ \beta_0 $的相合估计, 即:
进一步, 有
这里
进一步, 为了实际应用中的计算问题, 我们为渐近方差$ A^{-1}\left\{\Sigma_1+\Sigma_2\right\}{A^{-1}} $提出了如下一种相合估计. 定义:
以及
易得
是$ A^{-1}\left\{\Sigma_1+\Sigma_2\right\}{A^{-1}} $的一个相合估计.
本节我们通过一系列模拟研究来展示上节中讨论的加权估计方法在有限样本量下的优良表现, 展示病例队列设计下加性风险回归方法的应用价值. 考虑如下加性风险模型:
设定参数真值$ \beta_1 = 0 $或$ 0.5 $, $ \beta_2 = 0 $或$ 0.5 $. 协变量$ Z_1 $分别从均匀分布$ U(0, 1) $和正态分布$ N(0, 1) $中生成, $ Z_2 $从成功率为0.5的Bernoulli分布中生成. 设定基准风险函数$ \lambda_{0}(t) = 1 $, 则失效时间$ \widetilde{T} $可以从风险率为$ \lambda_0(t)+\beta_1Z_1+\beta_2Z_2 $的指数分布中生成. 删失时间$ C $从均匀分布$ U[0, c] $中生成, 通过挑选$ c $的不同取值从而产生相应的删失率, 分别为$ \rho $ = 80%, 85% 和90%. 对于病例队列设计, 设定全队列样本总量$ N = 1000 $, 子队列样本量为$ n_0 = 200 $.
为了阐明问题, 我们比较以下几种方法:
Full: 基于全队列的估计方程方法($ \widehat{\beta}_F $);
SRS: 基于子队列的估计方程方法($ \widehat{\beta}_S $);
Naive: 基于与病例队列样本相同样本量的简单随机样本的估计方程方法($ \widehat{\beta}_N $);
CC: 病例队列设计下的逆概率加权估计法($ \widehat{\beta}_P $).
在每种参数设定下, 比较上述四种方法的参数估计值的均值(Mean), 估计值的样本标准差(SD), 标准差估计值的均值(SE), 95% 的正态区间覆盖率(CP) 以及估计的相对效率(RE), 结果均基于1000次的模拟结果计算获得. 模拟结果请见表 1和表 2.
在所有考虑的情况下, 关于$ \beta_1 $和$ \beta_2 $的四种估计都是无偏的, 标准误差估计的均值(SEs) 很好地估计了估计值的样本标准差(SDs), 置信区间覆盖率(CPs) 均约为$ 95\% $. 模拟结果表明, 较之传统简单随机抽样设计, 病例队列设计能有效地提高估计的效率. 在所有情况下, $ \widehat{\beta}_P $均比$ \widehat{\beta}_S $和$ \widehat{\beta}_N $更有效. 例如: 关于$ \beta_1 $的估计, 当$ \rho = 0.9, \ \beta_1 = \beta_2 = 0, \ Z_1\sim U(0, 1) $时, 相较于$ \widehat{\beta}_F $, $ \widehat{\beta}_S $, $ \widehat{\beta}_N $和$ \widehat{\beta}_P $的相对效率分别为$ 0.20 $, $ 0.28 $和$ 0.62 $. 这一方面说明, 在相同样本量下, 病例队列设计的效率约为简单随机抽样设计的$ 2.2 $倍. 另一方面说明, 病例队列设计仅用了全队列$ 28\% $的样本量可达到约$ 62\% $的效率. 再例如, 关于$ \beta_2 $的估计, 当$ \rho = 0.85, \ \beta_1 = 0, \ \beta_2 = 0.5, \ Z_1\sim N(0, 1) $时, 相较于$ \widehat{\beta}_F $, $ \widehat{\beta}_S $, $ \widehat{\beta}_N $和$ \widehat{\beta}_P $的相对效率分别为$ 0.20 $, $ 0.32 $和$ 0.54 $. 病例队列设计仅用了全队列$ 32\% $的样本量却达到了约$ 54\% $的效率, 且其效率为相同样本量下简单随机抽样设计的$ 1.7 $倍. 当删失率提高时, $ \widehat{\beta}_P $的效率更高.
总体来说, 在有限样本量下, 本文研究的病例队列设计下加性风险模型中的加权估计方法表现优异. 病例队列设计作为一种基于生存数据的有偏抽样机制, 通过将资源集中到认为包含有更多信息的群体上, 能够提高估计的效率并节约研究成本. 因此, 相比简单随机抽样设计有显著的成本效益. 尤其是感兴趣的事件是稀发事件时, 此时抽样设计非常高效.
本节我们研究乳腺癌相关数据. 数据来源于美国国家癌症研究所(National Cancer Institute)[22]. 我们选取了1975-2017年期间$ 40 $岁以上的女性乳腺癌患者的共$ 118763 $条数据. 我们基于加性风险模型对此数据集进行生存回归分析, 探索影响乳腺癌患者生存时间的主要影响因素. 进一步, 我们应用病例队列设计分析数据集, 展示此种有偏抽样机制的实际应用价值.
我们感兴趣的因变量是乳腺癌患者的生存时间, 而观测到的生存时间存在删失, 其删失率为$ 84.4\% $. 我们考虑如下6个潜在影响因素. 患者年龄(Age) 分为5组: 40-49岁(Age = 1), 50-59岁(Age = 2), 60-69岁(Age = 3), 70-79岁(Age = 4), 80岁以上(Age = 5). 种族(Race) 分为3种: 其他种族(Race = 1), 白种人(Race = 2), 黑种人(Race = 3). 癌细胞分化程度(Grade) 分为4种类型: 肿瘤高度分化(Grade = 1), 肿瘤中度分化(Grade = 2), 肿瘤低分化(Grade = 3), 肿瘤未分化(Grade = 4). 肿瘤直径大小的T分期(T) 分为5种类型: 肿瘤直径$ \leq 2cm $(T = 1), $ 2cm < $肿瘤直径$ \leq 5cm $(T = 2), 肿瘤直径$ >5cm $(T = 3), 肿瘤直接侵犯胸壁或皮肤(T = 4), 肿瘤无法评估(T = 5). 是否并发淋巴癌的N分期分为5种类型: 同侧腋窝无肿大淋巴结(N = 1), 同侧腋窝有尚可推动的肿大淋巴结(N = 2), 同侧腋窝肿大淋巴结融合或粘连(N = 3), 有同侧胸骨旁淋巴结转移(N = 4), 区域淋巴结无法评估(N = 5). 肿瘤是否转移的M分期(M) 分为2种类型: 肿瘤未转移(M = 1), 肿瘤转移(M = 2). 我们对数据中考虑的上述影响因素进行了描述性统计分析, 画出了其条形图(图 1) 及生存曲线图(图 2).
为了探究乳腺癌患者生存时间的影响因素. 我们建立如下加性风险模型:
首先, 基于全队列$ 118763 $条数据进行加性风险回归分析, 结果总结在表 3中(见Full-Cohort). 结果表明, 考虑的$ 6 $个协变量对乳腺癌患者的生存时间均有显著的影响. 具体来说, 年龄越大的患者死亡率越高. 黑种人死亡风险率最高, 白种人次之, 其他人种最低. 癌细胞分化程度越高, 肿瘤大小的T分期等级越高, 区域淋巴癌的N分期等级越高, 患者的生存率越低. 癌细胞有远处转移的患者比无远处转移的患者的死亡风险率高出$ 1.476\% $.
为了评估病例队列设计并展现其在实际应用中的可行性与有效性, 我们首先从全队列中随机抽取了一个容量为$ n_0 = 30000 $的子队列, 然后将子队列和子队列之外的病例($ n_c = 13903 $) 组成病例队列样本, 其样本容量为43903. 我们应用本文研究的病例队列设计下加性风险模型下的加权估计方法对数据进行分析, 其结果总结在表 3中(见Case-Cohort). 结果表明, 病例队列设计采用了较小的样本量, 但估计结果与基于全队列数据分析得到的估计结果十分接近. 病例队列设计仅用了全队列约37% 的样本量, 对于Grade, T, N的估计分别达到了约$ 52.6\% $, $ 41.2\% $和$ 45.8\% $的效率. 这说明在实际应用中, 当影响因素Grade, T和N的测量非常昂贵或耗时时, 采用病例队列设计会提高研究效率和节约成本.
基于前人的研究工作(Lin, 1994[3], Kang et al, 2013[18]), 我们在此给出定理1的证明.
定理1的证明 根据一致强大数定律(Pollard, 1990[26]), 对于$ t\in[0, \ \tau] $一致地有:
由此可得, 对于$ t\in[0, \ \tau] $一致地成立:
注意到: 在模型(2.1) 下, $ M^{(0)}_i(t) $是均值为0的随机过程. 因此, 我们有
由(A.1) 和(A.2) 式, 可得:
由此可得$ \widehat{\beta}_w $的相合性.
进一步, 我们有
对于(A.3) 式右侧的第二项, 由于零均值过程$ \omega_iM_{i}^{(0)}(t) $可以表示成两个单调过程的和, 因此, 根据单调随机过程的中心极限定理(van der Vaart & Wellner, 1996)[27], 我们有$ \frac{1}{\sqrt{N}}\sum_{i = 1}^{N}\omega_iM_{i}^{(0)}(t) $收敛到一个在$ [0, \tau] $上具有连续通道的紧高斯过程. 故而, 由Kulich & Lin (2000)[28] 中的引理A.1, 可得: (A.3) 式右边第二项收敛到零.
因此, 我们有
其中,
由中心极限定理, 易知,
由此可得,
基于(A.4) 式以及Slutsky定理, 我们可得如下渐近正态性: