竞争风险下的置信区间：CIF 与 Kaplan–Meier 的统计分析与实现规范

1问题界定：SAP 与 Shell 到底要求了什么

SAP 第 5.10.1 节把转移发生时间的分析拆成主分析与支持性分析两条线，并明确规定了两条线各自的估计量。主分析写道，Cumulative Incidence Function 用于估计转移发生时间及各时间点的累积发生率，且「death prior to the development of metastatic disease treated as a competing risk」；未发生转移且未死亡者，在末次转移评估日删失。该段随后要求「25th percentile, median and 75th percentile of the time to occurrence of metastatic disease (years) along with 2-sided 95% confidence intervals (CIs) will be provided」，并要求 1、2、3、4、5 年的累积发生率「along with 2-sided 95% CIs」。支持性分析则写得明确得多：分位数使用 Brookmeyer–Crowley method with log-log transformation，年度无病率使用 Greenwood method with log-log transformation。

把这两段并置起来，最关键的缺口就出现了：支持性分析（Kaplan–Meier 分支）写明了两处 CI 的构造方法，而主分析（CIF 分支）只写了「要输出 95% CI」，没有写这些 CI 用什么方差、施加什么变换。Shell 层面同样没有补上：表 14.2.1.2 的脚注只定义了时间变量的口径，没有 CI 方法的脚注；表 14.2.1.3 的三个分位数单元格则带脚注 [2]，明确写着「95% CI is computed using the Brookmeyer–Crowley method」。也就是说，同一个 Shell 文件里的两张表，一张有 CI 方法、一张没有。这不是排版疏漏，而是一个必须在编程前用方法学决定填补的空白，因为它直接决定 8 处 CI 的数值。

2为什么 CIF 与 KM 的 CI 不能互换

两种方法的差异不在估计公式的写法，而在对「转移前死亡」这一观测的处置。理解这一点之后，后续所有 CI 的差异都可以从中推导出来，而不必逐条记忆。

2.1 风险集的口径：谁在分母里，谁不在

在讨论「死亡是否计入分母」之前，必须先固定风险集的定义。定义一旦固定，两条线的分母是否相同就不再是需要争论的问题——它们用的是同一个风险集。这是本文件后续所有结论的地基，因此单独说明。

记 Y(t) 为时点 t 处的风险集人数，它由在 t 时刻仍在随访、且此前未发生任何事件的受试者构成。这里「事件」不区分类型：转移、转移前死亡、以及任何已定义的终点，只要发生过，受试者在事件时点之后即退出风险集；删失同理，受试者在自己的删失时点仍留在集合内，此后退出。Wolbers 等在其图注中把这个集合写得非常直白：「the number of subjects who are still at risk, i.e. under follow-up and free of any event, at each time point」¹²。

为什么不能用「最终结局」筛分母

第 ② 条容易被反向理解：既然这些人最终会死于转移之外的原因，为什么不干脆把他们从分母里整批剔除？三条理由，任何一条都足以否掉这个做法。

• 可预测性。计数过程理论要求风险集在时点 t 处仅由 t 之前已观测到的信息决定，即必须是可预测的（predictable）。「此人最终会死于其他原因」在 t 时刻属于未来信息；一旦用它筛分母，估计量的概率含义就不再是「真实世界的累积发生率」。
• 目标量会被换掉。整批剔除会使分母被系统性削小，算出来的不再是 P(T ≤ t, 原因 = 1) 的 Aalen–Johansen 估计，而是一个条件于未来的量。它既没有对应的临床问题，也无法在 SAP 中定义清楚。
• 恒等式会破。只有两个累积发生率共用同一个全因风险集时，才有 1 − Ŝ(t) = F̂₁(t) + F̂₂(t)。换了分母，第 2.2 节那条最省事的 QC 判据立刻失效。

KM 与 CIF 的分母完全相同

这是本节最需要记住的一句。SAS 官方文档在 CIF 的估计式里只定义了一个风险集：设 t₁ < … < tL 为不同的未删失时间，对每个 l 令 Yl 为「the number of subjects at risk at tl」，令 djl 为 tl 处原因 j 的失败数；而估计式中出现的生存函数被明确定义为「the Kaplan-Meier estimator that would have been obtained by assuming that all failure causes are of the same type」⁴——正是 KM 分支所用的、把所有原因合并后的那一条 KM。同一个 Yl、同一条全因 KM，唯一的差别在分子怎么记。

分母相同这一点，在软件输出里还有一处直接证据：OUTCIF= 数据集同时给出 AtRisk（「the number of subjects at risk just before the specified time」）、Event（该时点关注原因的失败数）与 AllEventTypes（该时点任何原因的失败数）三列²。一张表里只有一个 AtRisk，却有两个失败数，正是「分母共享、分子分列」的软件层写照。（列名请以 proc contents 实测为准。）

2.2 恒等式：KM 的高估量恰好等于竞争事件的累积发生率

把这个恒等式写清楚，两种估计量的关系就不再是「一个高估另一个低估」的模糊印象了。设转移为第 1 类事件、转移前死亡为第 2 类事件，则 Aalen–Johansen 估计量满足 F̂₁(t) + F̂₂(t) = 1 − Ŝ(t)，其中 Ŝ 是把所有原因合并后的 Kaplan–Meier 估计。于是 1 − Ŝ(t) − F̂₁(t) 恒等于 F̂₂(t)，也就是竞争事件的累积发生率。Kaplan–Meier 分支所报告的无病率损失之所以系统性偏大，偏大的部分正是把本应显式建模的死亡算成了「转移的可能性」。这一结论不依赖任何模拟数据，可由代码直接验证，因此也应当成为交付前的一项例行 QC。^9,13,26

2.3 方差结构：为什么 CIF 不能沿用 Greenwood 方差

方差结构上的差异同样值得单独说明。Ŝ(t) 的方差由 Greenwood 公式给出，其推导只用到「单一事件类型」的乘积结构；而 F̂₁(t) 是 Ŝ(u−) 与各时点转移增量 ΔN₁(u)/Y(u) 的加权和，因而其方差必须同时承载两部分的随机性，并考虑两者之间的协方差。这就是为什么 CIF 的置信区间不能沿用 Greenwood 方差，也是 PROC LIFETEST 为什么把方差选项 ERROR= 明确限定为「calculation of the variance of the CIF estimator」¹，而生存函数一侧另有一套固定机制。

3置信区间的两种构造机制

本项目中出现的全部 CI 都可归入两种机制之一：对估计值施加变换后再反变换的变换法，以及把逐点置信限倒过来用的反演法。把这两种机制分开看，Shell 里 8 处 CI 就不再是一堆彼此无关的公式。

3.1 变换法：方差 + 变换，作用于时间点

变换法的通式是先把估计量 θ̂ 送进一个单调变换 g，在变换后的尺度上用 delta 法近似正态，取 ± z 倍标准误，再反变换回原尺度：

CI = [ g⁻¹( g(θ̂) − z · σ̂_g ) ,  g⁻¹( g(θ̂) + z · σ̂_g ) ] ,
其中  σ̂_g = | g′(θ̂) | · σ̂(θ̂)      （delta 法）
θ̂ = Ŝ(t)（生存函数一侧）或 F̂₁(t)（CIF 一侧）
z  = 1 − α/2 分位点（α = 0.05 时为 1.959964）

SAS 用 CONFTYPE= 这一个选项来切换 g。五种取值在文档中被逐条列出，默认值为 LOGLOG¹；对 CIF 一侧，OUTCIF= 数据集中的 CONFTYPE 列被明确定义为「the name of the transformation that is applied to the CIF to compute the confidence intervals for the CIF」²。表 3 给出五种变换的完整对照。

3.2 反演法：把逐点置信限倒过来用，作用于时间点

分位数 CI 不引入新的方差公式，而是把 3.1 节的置信限反过来当作检验统计量使用。设 F(t) 为所研究的分布函数（KM 分支取 1 − S(t)，CIF 分支取 F₁(t)），l(t)、u(t) 为其逐点置信下限与上限，则 p 分位数的置信集定义为

{ t : l(t) ≤ p ≤ u(t) }
在 l(t)、u(t) 同为单调阶梯函数时，该集合为区间 [ qL , qU )，其中
qL = inf{ t : u(t) ≥ p }     ← 由置信上限首次触及 p 确定
qU = inf{ t : l(t) ≥ p }     ← 由置信下限首次触及 p 确定

这一构造即 Brookmeyer 与 Crowley 提出的方法¹¹，也正是 SAS 在 Quartiles 表中给出分位数 CI 的依据：CONFTYPE= 的文档写得很清楚，该选项除了作用于生存函数的逐点 CI 与置信带之外，「in addition to the confidence intervals for the quartiles of the survival times」¹。换言之，SAP 对 KM 分支写明的「Brookmeyer–Crowley + log-log」，在 SAS 中就是 CONFTYPE=LOGLOG 下的 Quartiles 表；而同一个反演逻辑可以逐字搬到 CIF 上，只是把 S 换成 F₁。

实现时有两点需要写进程序注释。第一，l(t)、u(t) 与 F̂(t) 必须来自同一次 PROC LIFETEST 调用、同一组选项，否则三个阶梯不自洽，qL ≤ q̂ ≤ qU 可能被破坏。第二，点估计与上下界在判据上统一采用「首次达到或超过」的写法，即可保证三者的单调次序；严格的 > 与 ≥ 之别只可能在阶梯跳变恰好落在 p 上的极罕见情形下产生一个时间点的差异，SAS 内建 Quartiles 表以自身实现为准，因此KM 分支应当以 Quartiles 表作为主口径，把自行反演的结果作为交叉核对；CIF 分支则没有这个选择，只能由自行反演给出，并在 SAP 中把规则写死。

这张表的用途不是推导数值，而是说明一件事：「置信区间」这三个字在本项目里不是唯一确定的，它取决于一个必须被写进 SAP 的选项。因此本文件在第 7 节给出的增补文字里，把 CONFTYPE=LOGLOG、ERROR=AALEN、ALPHA=与 ALPHAQT= 都写到了参数层面。

3.3 两类机制之外的第三种：同时置信带

为完整起见需说明，OUTSURV= 还可以输出 同时置信带（CONFBAND=EP 的等精度带、CONFBAND=HW 的 Hall–Wellner 带、CONFBAND=ALL 两者）¹。它与前述逐点 CI 的区别在于推断对象是整条曲线而不是单个时间点，覆盖概率的控制对象不同，因此数值必然更宽，不能与逐点 CI 混用，也不能通过 OUTCIF= 获取（该选项只提供逐点限）。本项目 Shell 未要求置信带，故不作为报表口径；仅在需要讨论曲线整体比较时才有意义。

4METHOD= 参数族：每一个选项究竟改变了什么

SAS 中「方法」并不是一个单一参数，而是分布在不同层次上的几个选项。把它们混在一起理解是常见的误解来源，最典型的一例是把 SAP 里写的「Greenwood method」当成某个 SAS 选项去查找——它其实不是选项，而是默认行为。

4.1 关于 NELSON：为什么它不在 METHOD= 的取值里

这是一个很容易记混的地方，且不只一个人记混。SAS 官方论坛上就有人写着「SAS references say to specify METHOD=NELSON within proc lifetest」而发帖求助；而文档里真正可用的写法是直接把 NELSON 当作 PROC LIFETEST 语句上的一个独立选项。两者结果上差别不小：写错了会报错，而不是静默给出一个不同的估计量。

也就是说，两个选项作用在不同的对象上：METHOD= 决定「用哪个估计量取估计生存函数」，NELSON 则决定「累积风险列用哪个估计量」。后者不修改 S(t)、不修改 S(t) 的置信区间、也不修改分位数，它只把输出中的——如果有的话——累积风险列从 −log S(t) 换成 Nelson–Aalen 估计值。

4.2 AALEN 与 DELTA：两个方差估计量在小样本下并不等价

由于 SAP 在 CIF 分支上没有规定方差，ERROR= 的默认值 AALEN 会在无人察觉的情况下成为项目口径，这一隐性默认值得评估。Aalen 型方差来自计数过程的鞅理论；delta 法则把 CIF 估计量视为多项分布的矩，用一阶展开近似。Braun 与 Yuan 在六种方差估计量与 bootstrap 之间做了系统的模拟比较，结论是：基于多项分布矩的估计量（Dinse 型 / delta 型）与 bootstrap 的表现接近，即使在 20 例的样本中仍相当准确；而基于鞅理论的估计量除一个例外，在少于 100 例的样本中倾向于系统性高估或低估经验方差¹⁸。Choudhury 的独立模拟同样支持 Dinse–Larson 型方差配合 log(−log) 变换¹⁷。

这一结论对本研究的直接含义是：若各队列的转移事件数较少（这在 registry 研究中很常见），ERROR=AALEN（即 SAS 默认）给出的 CIF 置信限可能偏离名义覆盖水平。稳健的做法是把 ERROR=DELTA 作为预先规定的敏感性分析同时跑一遍，比较两组 CIF_STDERR 的相对差异；若差异可忽略，则在 SAP 中以默认值定稿并记录该核对结论；若差异明显，则应把方差选择本身作为一个待申办方与统计审阅方确认的决策项，而不是留作默认。两个方差的具体公式见 SAS 文档「Estimation of the CIF」一节⁴，其中 delta 法那一项的写法被作者描述为「in the spirit of Greenwood's formula」。

5CIF 分位数为何没有现成输出

这是本项目最容易卡住一处：SAP 要求输出 CIF 的 25、50、75 分位数及 95% CI，而 SAS 无论如何设置选项，都不会给出这三个数。原因分三层，只有第一层是软件层面的，后两层是统计层面的。

5.1 软件层：Quartiles 表根本不属于 CIF

SAS 官方文档列出了 PROC LIFETEST 产生的全部 ODS 表及其适用语句⁵。与竞争风险相关的表只有三张：CIF（描述为「Cumulative incidence function estimates」）、FailureSummary（「Summary of failure outcomes for competing-risks data」）与 GrayTest（「Results of k-sample test of Gray (1988) comparing CIFs」）。而 Quartiles 表的描述是「Quartiles of the survival times」，其适用条件被限定为 METHOD=PL | B | FH——也就是说，它挂的是生存函数这条线，与 CIF 无关。OUTCIF= 数据集的列清单里只有风险集与失败数三列（AtRisk、Event、AllEventTypes，见 2.1 节）加上 CIF、CIF_STDERR、ALPHA、CONFTYPE、CIF_LCL、CIF_UCL，没有分位数列²。图形侧的证据一致：PLOTS=CIF 与 PLOTS=CIF(CL) 生成 cifPlot，但它不传递 Median、LowerMedian、UpperMedian 这几个动态变量，而 SurvivalPlot 会传⁶。软件没有欠缺什么，它只是从未把 CIF 分位数当作一个输出对象。

5.2 统计层之一：CIF 不是真分布函数，分位数可能不存在

累积发生率函数是非真（improper）分布函数：F₁(∞) = P(发生转移) < 1，其渐近水平（平台值）等于该组最终发生转移的比例。因此它的逆函数只在平台值以下有定义。文献对此的表述很直接：CIF 的渐近上界小于 1，意味着它是非真函数，定义均值没有意义（恒为无穷），而「分位数有限、且可能从观测数据中识别」，因此被推荐作为刻画 CIF 曲线的摘要量¹⁹；p 分位数定义为 F₁⁻¹(p) = inf{ t : F₁(t) ≥ p }¹⁹。反过来看，当平台值低于 p 时，该分位数在数学上根本不存在，任何软件都不可能给出——包括 SAS。这解释了为什么 Shell 会给分位数留出三个单元格却没有任何提示：在本研究的规模下，「中位转移时间为 NE」是完全可以预期的正常结果，而不是数据质量问题。

5.3 统计层之二：CIF 分位数的推断方法本身就成熟得较晚

还有一层历史原因值得知道：CIF 分位数的非参数推断直到 Peng 与 Fine 的工作才被系统建立，他们提出基于累积发生率函数定义分位数，并证明了该分位数函数非参数估计量的一致性与弱收敛，由此得到置信区间、置信带与两样本检验¹⁹。Lee 与 Fine 随后给出简化版非参数推断与参数化版本，并在经验研究中比较两者，发现在模型误设不严重时参数法均方误差更小²⁰。也就是说，这是一个方法学上较晚成熟的领域，通用统计软件没有把它做成默认输出，是有其背景的，而不是实现疏漏。这也意味着：如果项目要使用超出「逐点置信限反演」的精巧方法，就必须在 SAP 中明确引用方法学文献，而不能期望软件默认承担。

5.4 四条可行路径及其取舍

6参考代码：统计分析部分

以下代码只包含统计分析逻辑，不含报表宏、页眉页脚、输出层格式化与骨架程序。输入为项目的 ADTTE（paramcd = 'TTMET'），其中 aval 为距父研究首次给药的时长、单位已是「年」（avalu = 'YEARS'），cnsrscd 取 1（转移事件）、2（转移前死亡，即竞争事件）、0（删失），anlgrp 为分组变量。所有选项均显式写出，不使用默认值的隐式依赖。基准环境为 SAS 9.4M7（TS1M7）/ SAS/STAT 15.2（第 8 节）；OUTCIF= 的列名请先用 proc contents 核一遍再改列，不要凭记忆编码。

/*------------------------------------------------------------------
  输入：项目 ADTTE（adtte.sas Step 6 的输出）
        paramcd = 'TTMET'
        aval    = 距父研究首次给药的时长，单位已为「年」(AVALU='YEARS')
        cnsrscd = 1 转移事件 / 2 转移前死亡(竞争事件) / 0 删失
  输出：tte
------------------------------------------------------------------*/
data tte;
    set adam.adtte;
    where paramcd = 'TTMET' and aval ne .;
    evt = (cnsrscd = 1);          /* KM 分支的事件指示：0 与 2 同为删失 */
    keep usubjid anlgrp aval cnsrscd evt;
run;

/* 方差：METHOD=KM 的默认 Greenwood（SAS 无 ERROR=GREENWOOD，无需也无法指定）
   变换：SAP 指定 log-log，即 CONFTYPE=LOGLOG（恰为默认，仍显式写出）
   分位数 CI：由 CONFTYPE 驱动的 Brookmeyer-Crowley 反演               */
proc lifetest data = tte
        method   = km            /* KM | BRESLOW | FH | ACT/LIFE/LT        */
        conftype = loglog        /* 施加于 S(t) 的变换；默认值              */
        alpha    = 0.05          /* 逐点 CI 显著性水平                      */
        alphaqt  = 0.05          /* 分位数 CI 显著性水平（独立于 ALPHA=）    */
        outsurv  = km_surv       /* SURVIVAL / SDF_STDERR / SDF_LCL / SDF_UCL */
        noprint;
    time aval  * cnsrscd(0 2);   /* 死亡(2)在 KM 分支按删失处理              */
    strata anlgrp;
run;

/* 2a. 1–5 年无病率(%)及其 95% CI
       阶梯右连续：取 time <= t 的最后一条观测                        */
proc sql;
    create table km_rate as
    select g.anlgrp,
           y.yr                 as yr,
           s.survival           as rate      label = 'KM 无病率',
           s.sdf_lcl            as rate_lcl  label = '95% CI 下限',
           s.sdf_ucl            as rate_ucl  label = '95% CI 上限'
    from (select distinct anlgrp from km_surv)                         as g,
         (select 1 as yr union all select 2 as yr union all select 3 as yr
          union all select 4 as yr union all select 5 as yr)           as y,
         km_surv                                                       as s
    where s.anlgrp = g.anlgrp
      and s.aval  <= y.yr
    group by g.anlgrp, y.yr
    having s.aval = max(s.aval);
quit;

/* 基准环境：SAS 9.4M7（TS1M7）/ SAS/STAT 15.2
   方差：ERROR= 仅作用于 CIF；默认 AALEN，DELTA 为可选（见块 5）
   变换：CONFTYPE= 施加于 CIF 以计算 CIF 的置信区间（SAS 文档原话）
   TIME 语句：eventcode=1 指定关注事件；cnsrscd 的其余非零值视为竞争事件
              15.1 及以后版本在文档中列作 FAILCODE / EVENTCODE，同一选项
   输出列：CIF / CIF_STDERR / ALPHA / CONFTYPE / CIF_LCL / CIF_UCL
            以及风险集计数 AtRisk / Event / AllEventTypes
   注意：t=0 的初始观测按 CIF(0)=0 解释（文档该句疑为复制笔误，见 2.1 节） */
proc lifetest data = tte
        conftype = loglog        /* 默认；SAP 未规定，此处显式写死          */
        error    = aalen         /* AALEN(默认) | DELTA                     */
        alpha    = 0.05
        outcif   = cif_est       /* CIF / CIF_STDERR / ALPHA / CONFTYPE /   */
        noprint;                 /* CIF_LCL / CIF_UCL                       */
    time aval  * cnsrscd(0) / eventcode = 1;
    strata anlgrp;
run;

/* 3a. 审计留痕：回读本次分析实际使用的变换与水平，写入 LOG */
data _null_;
    set cif_est(obs = 1 keep = conftype alpha);
    put 'NOTE: [QC] CIF CI transformation actually applied = ' conftype
        ', alpha = ' alpha;
run;

/* 3b. 1–5 年累积发生率(%)及其 95% CI */
proc sql;
    create table cif_rate as
    select g.anlgrp,
           y.yr                as yr,
           s.cif               as cif       label = 'CIF 累积发生率',
           s.cif_lcl           as cif_lcl   label = '95% CI 下限',
           s.cif_ucl           as cif_ucl   label = '95% CI 上限'
    from (select distinct anlgrp from cif_est)                         as g,
         (select 1 as yr union all select 2 as yr union all select 3 as yr
          union all select 4 as yr union all select 5 as yr)           as y,
         cif_est                                                       as s
    where s.anlgrp = g.anlgrp
      and s.aval  <= y.yr
    group by g.anlgrp, y.yr
    having s.aval = max(s.aval);
quit;

/* 3c. Event / Death / Censored 计数（表头三行） */
proc freq data = tte noprint;
    tables anlgrp * cnsrscd / out = cif_counts;
run;

/* 4.1 把两条线统一成同构的阶梯三元组 (f, l, u)
       KM  ：f = 1 - S(t) 非增，故 f 的下限 = 1 - S(t) 的上限
       CIF ：f = F1(t)    非减，直接取 OUTCIF 的两个限                */
data step_km;
    set km_surv(rename = (survival = s sdf_lcl = sl sdf_ucl = su));
    f = 1 - s;                  /* 补生存函数(非增)      */
    l = 1 - su;                 /* f 的下限 = 1 - S 上限 */
    u = 1 - sl;                 /* f 的上限 = 1 - S 下限 */
    keep anlgrp aval f l u;
run;

data step_cif;
    set cif_est(rename = (cif = f cif_lcl = l cif_ucl = u));
    keep anlgrp aval f l u;
run;

data step_all;
    set step_km(in = in_km) step_cif(in = in_cif);
    length src $3;
    if in_km  then src = 'KM';
    if in_cif then src = 'CIF';
run;

proc sort data = step_all; by src anlgrp aval; run;

/* 4.2 Brookmeyer-Crowley 反演
       q  = inf{ t : f(t) >= p }
       qL = inf{ t : u(t) >= p }        u = f 的置信上限
       qU = inf{ t : l(t) >= p }        l = f 的置信下限
       依据 { t : l(t) <= p <= u(t) } = [ qL , qU )
       自检条件：qL <= q <= qU 恒成立                                   */
data quant;
    set step_all;
    by src anlgrp aval;
    array pct{3} _temporary_ (0.25 0.50 0.75);
    array qq {3} q25  q50  q75;
    array ql {3} l25  l50  l75;
    array qu {3} u25  u50  u75;
    array fq {3} $1 fq25 fq50 fq75;   /* 点估计是否可估 Y/N */
    array fl {3} $1 fl25 fl50 fl75;
    array fu {3} $1 fu25 fu50 fu75;
    retain q25-q75 l25-l75 u25-u75 fq25-fq75 fl25-fl75 fu25-fu75;
    if first.anlgrp then do j = 1 to 3;
        qq{j} = .; ql{j} = .; qu{j} = .;
        fq{j} = 'N'; fl{j} = 'N'; fu{j} = 'N';
    end;
    do j = 1 to 3;
        if missing(qq{j}) and f >= pct{j} then do; qq{j} = aval; fq{j} = 'Y'; end;
        if missing(ql{j}) and u >= pct{j} then do; ql{j} = aval; fl{j} = 'Y'; end;
        if missing(qu{j}) and l >= pct{j} then do; qu{j} = aval; fu{j} = 'Y'; end;
    end;
    if last.anlgrp then output;
    keep src anlgrp q25-q75 l25-l75 u25-u75 fq25-fq75 fl25-fl75 fu25-fu75;
run;

/* 4.3 自检：三者的单调次序与「不可估」标记 */
data quant_chk;
    set quant;
    array qq{3} q25  q50  q75;
    array ql{3} l25  l50  l75;
    array qu{3} u25  u50  u75;
    array fq{3} $1 fq25 fq50 fq75;
    array fl{3} $1 fl25 fl50 fl75;
    array fu{3} $1 fu25 fu50 fu75;
    do j = 1 to 3;
        if fq{j} = 'Y' and fl{j} = 'Y' and ql{j} > qq{j} then
            put 'WARNING: [QC] qL > q, src=' src ' anlgrp=' anlgrp ' j=' j;
        if fq{j} = 'Y' and fu{j} = 'Y' and qu{j} < qq{j} then
            put 'WARNING: [QC] qU < q, src=' src ' anlgrp=' anlgrp ' j=' j;
    end;
run;

/* 5.1 KM 分支交叉核对：LIFETEST 内建的 Quartiles 表是 Brookmeyer-Crowley 的
       官方实现，应以它为主口径，与 4.2 的自行反演结果比对。
       该表列名随版本略有差异，先取 CONTENTS 再写比对逻辑，切勿凭记忆编码。 */
ods output quartiles = km_quart;
proc lifetest data = tte
        method   = km
        conftype = loglog
        alphaqt  = 0.05
        noprint;
    time aval  * cnsrscd(0 2);
    strata anlgrp;
run;
ods output close;

proc contents data = km_quart varnum; run;   /* 确认列名后再比对 */

/* 5.2 敏感性分析：改用 delta 法方差，比较 CIF_STDERR 与置信限的差异。
       参照 Braun & Yuan (2007)：Aalen 型方差在 n < 100 时可能系统性偏倚，
       Dinse / delta 型更接近 bootstrap。                                     */
proc lifetest data = tte
        conftype = loglog
        error    = delta
        alpha    = 0.05
        outcif   = cif_delta
        noprint;
    time aval  * cnsrscd(0) / eventcode = 1;
    strata anlgrp;
run;

data cif_var_cmp;
    merge cif_est  (in = a keep = anlgrp aval cif cif_stderr cif_lcl cif_ucl
                            rename = (cif_stderr = se_a cif_lcl = l_a
                                      cif_ucl = u_a cif = cif_a))
          cif_delta(in = b keep = anlgrp aval cif cif_stderr cif_lcl cif_ucl
                            rename = (cif_stderr = se_d cif_lcl = l_d
                                      cif_ucl = u_d cif = cif_d));
    by anlgrp aval;
    if a and b;
    se_ratio = se_d / se_a;                     /* 方差口径的相对差异 */
    w_a      = u_a - l_a;                       /* Aalen 区间宽度     */
    w_d      = u_d - l_d;                       /* Delta 区间宽度     */
    w_ratio  = w_d / w_a;
run;

proc means data = cif_var_cmp n min p10 median p90 max maxdec = 3;
    var se_ratio w_ratio;
run;
/* 5.3 风险集（分母）一致性核对
       两条分支的风险集是同一个 Y(t)（见 2.1 节），因此 AtRisk 应当等于
       「观察时间 >= 该时点」的人数。两者不等，问题一定在分析集、
       时点取数或分层定义，而不是「两种方法的风险集算法不同」。
       说明：OUTSURV= 不含 at-risk 列，只能用推算值去比，不能拿另一列去凑。
       t=0 那一行是最好的烟点：AtRisk 应等于分析集总人数。   */
proc sql;
    create table atrisk_chk as
    select c.anlgrp,
           c.aval                        as aval,
           c.atrisk                      as atrisk_sas,
           count(t.usubjid)              as atrisk_chk,
           c.atrisk - count(t.usubjid)   as diff
    from cif_est as c
         left join tte as t
              on  t.anlgrp = c.anlgrp
              and t.aval  >= c.aval
    group by c.anlgrp, c.aval, c.atrisk
    having c.atrisk ne count(t.usubjid);   /* 应输出 0 行；有行即说明分母口径有误 */
quit;

data _null_;
    if 0 then set atrisk_chk nobs = n;
    put 'NOTE: [QC] atrisk mismatch rows = ' n;   /* 需为 0 */
run;

/* 5.4 初始观测断言：CIF(0) 必须为 0（文档 t=0 那句写作 1，疑为笔误） */
data _null_;
    set cif_est(where = (aval = 0));
    if cif ne 0 then put 'WARNING: [QC] CIF(0) ne 0, cif=' cif;
run;

6.1 可运行的 worked example：以 CIF 数据集列表为输入做反演

上面的块 4 是生产代码，输入是整个分析集、需先跑 PROC LIFETEST 生成阶梯函数。这里再给一个自包含、可手算核对的独立示例：直接从一个已经算好的 CIF 数据集（即 OUTCIF= 归一化后的形态——四列 time / F / F_lcl / F_ucl）出发，走完整条反演链路。它能单独提交、单独验证，用来向审阅者逐位证明「反演」这一步做了什么，也可作为新统计程序员的上手样板。完整文件见 cif_inversion_example.sas。

输入：一份固定的 CIF 阶梯数据集

下面 8 行就是 OUTCIF= 等价形态（列名已归一化：F 即 CIF、F_lcl 即 CIF_LCL、F_ucl 即 CIF_UCL）。时间单位取「月」便于读数。注意三条序列都非降，这是反演成立的前提；首行 t=0 的缺失 CIF 行已在上游剔除。

反演三式

对非降的 F，三条曲线的「第一次达到 p 的时点」就是反演结果；方向与生存函数分支相反（F 递增，故下界由置信上限驱动）：

因为三条曲线都是非降阶梯函数，每个 inf{ t : 曲线 ≥ p } 都退化为「从上往下扫，取第一次满足条件的行」——无需求根、无需插值。这正是 SAS 自带的 KM Quartiles 表（Brookmeyer–Crowley 反演）同一套机制，只是这里把「单组、固定 3 分位」平铺成能逐行读的 data step。

完整 SAS 程序（自包含、可直接提交）

/*=======================================================================
  CIF PERCENTILE 95% CI BY INVERTING THE CIF CONFIDENCE LIMITS
  Self-contained, hand-checkable example. Input = a fixed CIF dataset
  (time / F / F_lcl / F_ucl), the normalised form of OUTCIF=.
=======================================================================*/
options nodate nonumber linesize=160 formdlim='-';

/*==== 1. INPUT: fixed CIF step function (single group "Demo") ==========*/
data cif_demo;
  input time F F_lcl F_ucl;
  length pop $8;
  pop = 'Demo';
  datalines;
1.0  0.08  0.02  0.22
2.0  0.15  0.06  0.34
3.0  0.21  0.10  0.41
4.0  0.27  0.14  0.48
5.0  0.32  0.18  0.55
6.0  0.36  0.21  0.60
7.0  0.39  0.24  0.64
8.0  0.41  0.26  0.67
;
run;

/*==== 2. THE INVERSION (non-decreasing F) ==============================*/
%let pslist = 0.25,0.50,0.75;
%let nq     = 3;

data cif_inverted;
  set cif_demo;
  by pop;

  array pctl[&nq.] _temporary_ (&pslist.);
  array tq[&nq.]   _temporary_;   /* q_hat : first t with F     >= p */
  array tl[&nq.]   _temporary_;   /* q_lcl : first t with F_ucl >= p */
  array tu[&nq.]   _temporary_;   /* q_ucl : first t with F_lcl >= p */
  array tp[&nq.]   _temporary_;   /* time just before the crossing   */
  array lv[&nq.]   _temporary_;   /* F(q_hat), the level at crossing */

  retain plateau plateau_lcl plateau_ucl maxtime;

  if first.pop then do;
    do i = 1 to &nq.;
      tq[i] = .; tl[i] = .; tu[i] = .; tp[i] = .; lv[i] = .;
    end;
    plateau = .; plateau_lcl = .; plateau_ucl = .; maxtime = .;
  end;

  maxtime     = max(maxtime,     time);
  plateau     = max(plateau,     F);
  plateau_lcl = max(plateau_lcl, F_lcl);
  plateau_ucl = max(plateau_ucl, F_ucl);

  do i = 1 to &nq.;
    if tq[i] = . then do;
      if F >= pctl[i] then do; tq[i] = time; lv[i] = F; end;
      else tp[i] = time;               /* still below p -- remember  */
    end;
    if tl[i] = . and F_ucl >= pctl[i] then tl[i] = time;
    if tu[i] = . and F_lcl >= pctl[i] then tu[i] = time;
  end;

  if last.pop then do i = 1 to &nq.;
    percentile = pctl[i];
    q_hat      = tq[i];
    q_lcl      = tl[i];
    q_ucl      = tu[i];
    q_prev     = tp[i];
    F_at_q     = lv[i];

    /* check 1: when all three exist, q_lcl <= q_hat <= q_ucl */
    chk_order = .;
    if n(q_hat, q_lcl, q_ucl) = 3 then chk_order = (q_lcl <= q_hat <= q_ucl);

    /* check 2: level at the crossing must be >= p (right-continuity) */
    chk_hit = .;
    if F_at_q ne . then chk_hit = (F_at_q >= percentile);

    output;
  end;

  keep pop percentile q_hat q_lcl q_ucl q_prev F_at_q
       plateau plateau_lcl plateau_ucl maxtime chk_order chk_hit;
  format percentile 4.2 F_at_q 5.3;
run;

/*==== 3. RESULTS =======================================================*/
title "C. Percentile inversion result (months)";
proc print data=cif_inverted noobs label;
  var pop percentile q_hat q_lcl q_ucl q_prev F_at_q
      plateau plateau_lcl plateau_ucl maxtime chk_order chk_hit;
run;
title;

手算核对与最终输出

以 p = 0.25 为例，逐行扫：CIF 在 t=4 才首次 ≥ 0.25（F(4)=0.27）→ q_hat=4；置信上限早在 t=2 就 ≥ 0.25（F_ucl(2)=0.34）→ q_lcl=2；置信下限要熬到 t=8 才 ≥ 0.25（F_lcl(8)=0.26）→ q_ucl=8。三端 2 ≤ 4 ≤ 8 正是「F 递增时置信限换边」的直观体现。完整三档结果如下（已用独立 Python 实现逐位复核，与 SAS 输出一致）：

7对 SAP 第 5.10.1 节的建议补充文字

以下文字可直接插入 SAP 5.10.1 节主分析段落之后，用于填补 CIF 分支的 CI 方法空白。写法上与 SAP 现有风格一致，仅补充方法、方差与「不可估」规则三项，不改变原有的估计目标与终点定义。

8版本对齐：SAS 9.4M7（TS1M7）与本文件的选项口径

本项目的基准环境为 SAS 9.4M7（TS1M7），随附 SAS/STAT 15.2。竞争风险相关功能并不受版本限制——它们自 SAS/STAT 14.1（随 SAS 9.4M3 发布）就已进入 PROC LIFETEST——因此本项目在功能层面没有障碍；真正需要对齐的是「文档口径」与「维护状态」两件事。

8.1 版本对应关系

8.2 本文件涉及的选项在 9.4M7 上的可用性

8.3 一条与交付有关的运维事实

9交付前核查清单

10参考文献

1. A　软件文档（SAS 9.4M7 / SAS/STAT 15.2）
2. SAS Institute Inc. The LIFETEST Procedure — PROC LIFETEST Statement（ALPHA=、ALPHAQT=、CONFBAND=、CONFTYPE=、ERROR=、METHOD=、NELSON | AALEN、CIFVAR、OUTCIF=、OUTSURV= 的取值与默认值）. SAS/STAT User's Guide. go.documentation.sas.com … /statug_lifetest_syntax01.htm
3. SAS Institute Inc. The LIFETEST Procedure — OUTCIF= Data Set（AtRisk、Event、AllEventTypes、CIF、CIF_STDERR、CIF_LCL、CIF_UCL、CONFTYPE、ALPHA）. … /statug_lifetest_details33.htm
4. SAS Institute Inc. The LIFETEST Procedure — OUTSURV= Data Set（SURVIVAL、SDF_STDERR、SDF_LCL、SDF_UCL、_CENSOR_、HW_*、EP_*；不含 at-risk 列）. … /statug_lifetest_details34.htm
5. SAS Institute Inc. The LIFETEST Procedure — Analysis of Competing-Risks Data / Estimation of the CIF（Aalen 方差、delta 方差、风险集 Yl 与全因 KM 的定义）. … /statug_lifetest_details25.htm
6. SAS Institute Inc. The LIFETEST Procedure — ODS Table Names（CIF、FailureSummary、GrayTest、Quartiles 的适用范围）. … /statug_lifetest_toc.htm
7. SAS Institute Inc. The LIFETEST Procedure — ODS Graphics（Table 7 图名清单与 cifPlot；Table 8 的 Median / LowerMedian / UpperMedian 仅属于生存曲线）. … /statug_lifetest_details86.htm
8. SAS Institute Inc. The LIFETEST Procedure — TIME Statement（FAILCODE / EVENTCODE）. SAS/STAT 15.1 User's Guide. support.sas.com/documentation/onlinedoc/stat/151/lifetest.pdf
9. SAS Institute Inc. SAS/STAT 与 SAS 9.4 维护版本的对应关系（15.2 (SAS 9.4M7)）与 SAS 9.4 各维护版本的支持级别（「SAS® 9.4M7 (TS1M7) … Limited Support began September 1, 2025」）. documentation.sas.com … Determining Your Update Path for SAS/STAT；support.sas.com … SAS 9.4 & Earlier Releases
10. B　临床应用与报告规范（优先采用较新文献）
11. Austin, P. C., Lee, D. S. & Fine, J. P. Introduction to the analysis of survival data in the presence of competing risks. Circulation 133, 601–609 (2016). doi:10.1161/CIRCULATIONAHA.115.017719 — 明确指出用 KM 的补集估计累积发生率会系统性偏高，且与竞争事件是否独立无关；附 R 与 SAS 代码
12. Austin, P. C. & Fine, J. P. Practical recommendations for reporting Fine–Gray model analyses for competing risk data. Stat. Med. 36, 4391–4400 (2017). doi:10.1002/sim.7501
13. Schuster, N. A., Hoogendijk, E. O., Kok, A. A. L., Twisk, J. W. R. & Heymans, M. W. Ignoring competing events in the analysis of survival data may lead to biased results: a nonmathematical illustration of competing risk analysis. J. Clin. Epidemiol. 122, 42–48 (2020). doi:10.1016/j.jclinepi.2020.03.004 — 以真实队列展示「KM 高估、CIF 应用」的全过程，适合作为对申办方的说明材料
14. Wolbers, M. et al. Competing risks analyses: objectives and approaches. Eur. Heart J. 35, 2936–2941 (2014). doi:10.1093/eurheartj/ehu131 — 其图注给出风险集的直白定义（「under follow-up and free of any event」）
15. Latouche, A., Allignol, A., Beyersmann, J., Labopin, M. & Fine, J. P. A competing risks analysis should report results on all cause-specific hazards and cumulative incidence functions. J. Clin. Epidemiol. 66, 648–653 (2013). doi:10.1016/j.jclinepi.2012.09.017
16. Andersen, P. K., Geskus, R. B., de Witte, T. & Putter, H. Competing risks in epidemiology: possibilities and pitfalls. Int. J. Epidemiol. 41, 861–870 (2012). doi:10.1093/ije/dyr213 — 给出「原因专属风险与累积发生率之间的一一对应关系丧失」这一核心论断
17. ICH. E9(R1) Addendum on Estimands and Sensitivity Analysis in Clinical Trials to the Guideline on Statistical Principles for Clinical Trials (International Council for Harmonisation, adopted 20 November 2019) — 竞争事件属于 intercurrent event，其处置方式必须在估计目标中预先声明
18. C　分位数与置信区间的构造（含方法原始出处）
19. Brookmeyer, R. & Crowley, J. A confidence interval for the median survival time. Biometrics 38, 29–41 (1982). — 分位数 CI 反演法的原始出处
20. Choudhury, J. B. Non-parametric confidence interval estimation for competing risks analysis: application to contraceptive data. Stat. Med. 21, 1129–1144 (2002). doi:10.1002/sim.1070
21. Braun, T. M. & Yuan, Z. Comparing the small sample performance of several variance estimators under competing risks. Stat. Med. 26, 1170–1180 (2007). doi:10.1002/sim.2661
22. Peng, L. & Fine, J. P. Nonparametric quantile inference with competing-risks data. Biometrika 94, 735–744 (2007). doi:10.1093/biomet/asm059
23. Lee, M. & Fine, J. P. Inference for cumulative incidence quantiles via parametric and nonparametric approaches. Stat. Med. 30, 3221–3235 (2011). doi:10.1002/sim.4349
24. Beyersmann, J., Di Termini, S. & Pauly, M. Weak convergence of the wild bootstrap for the Aalen–Johansen estimator of the cumulative incidence function of a competing risk. Scand. J. Stat. 40, 387–402 (2013)
25. Bluhmki, T. et al. A wild bootstrap approach for the Aalen–Johansen estimator. Biometrics 74, 977–985 (2018)
26. Aalen, O. O. & Johansen, S. An empirical transition matrix for non-homogeneous Markov chains based on censored observations. Scand. J. Stat. 5, 141–150 (1978). — Aalen–Johansen 估计量的原始出处
27. Aalen, O. Nonparametric estimation of partial transition probabilities in multiple decrement models. Ann. Stat. 6, 534–545 (1978). — SAS 文档将 ERROR=AALEN 的方差理论标注为「Aalen 1978」
28. Marubini, E. & Valsecchi, M. G. Analysing Survival Data from Clinical Trials and Observational Studies (Wiley, Chichester, 1995). — SAS 文档将 ERROR=DELTA 的方差标注为「Marubini and Valsecchi 1995」
29. Andersen, P. K., Borgan, Ø., Gill, R. D. & Keiding, N. Statistical Models Based on Counting Processes (Springer, New York, 1993). — 计数过程框架下风险集与可预测性的标准参考
30. Putter, H., Fiocco, M. & Geskus, R. B. Tutorial in biostatistics: competing risks and multi-state models. Stat. Med. 26, 2389–2430 (2007)
31. Gray, R. J. A class of K-sample tests for comparing the cumulative incidence of a competing risk. Ann. Stat. 16, 1141–1154 (1988). — SAS GrayTest 表的依据
32. Lin, D. Y. Non-parametric inference for cumulative incidence functions in competing risks studies. Stat. Med. 16, 901–910 (1997)
33. Allignol, A., Schumacher, M. & Beyersmann, J. Empirical transition matrix of multi-state models: the etm package. J. Stat. Softw. 38(4), 1–15 (2011)
34. Guo, C. & So, Y. Cause-specific analysis of competing risks using the PHREG procedure. In Proceedings of the SAS Global Forum 2018 Conference, Paper 2159-2018 (SAS Institute Inc., Cary, NC, 2018). support.sas.com/resources/papers/proceedings18/2159-2018.pdf