全书目录 / 实战案例

生存分析：KM、log-rank 与 Cox

PROC LIFETEST + PROC PHREG 是肿瘤项目 SAS 程序员的日常。这个案例用 R 的 survival 包在同一份经典公开数据上复刻全流程：KM 估计 → log-rank 检验 → Cox 多变量模型 → KM 图落盘，一路给你 SAS 对照坐标。

背景与学习目标

数据是 NCCTG 肺癌试验（survival::lung），228 例晚期肺癌患者的随访记录，是生存分析教学的事实标准数据集。脚本为 cases/c2_survival.R，仓库根目录下 Rscript cases/c2_survival.R 即可复现。学完你应当能够：

目标 1：用 Surv() 构造生存对象、survfit() 估计分组 KM 曲线并读懂中位生存期与置信区间（对照 PROC LIFETEST）。
目标 2：用 survdiff() 做 log-rank 检验比较组间生存分布。
目标 3：用 coxph() 拟合多变量 Cox 模型（对照 PROC PHREG），正确解释 HR 的方向与参考水平，并用 survminer 输出带 risk table 的 KM 图。

数据源：survival::lung — NCCTG（North Central Cancer Treatment Group）肺癌试验数据，228 例晚期肺癌患者，含生存时间、删失状态、年龄、性别、ECOG 评分（ph.ecog）等变量 · 获取方式：R 包 survival 内置（data(lung) 或直接引用），随 R 安装即得、无需联网 · 许可：公开教学数据，广泛用于教材与课程。

步骤 1：数据与清洗

# C2：公开肿瘤试验生存分析（survival::lung，NCCTG 肺癌试验）
# 运行：Rscript cases/c2_survival.R
suppressPackageStartupMessages({
  library(survival); library(survminer); library(ggplot2)
})

cat("== 步骤 1：数据与清洗 ==\n")
dat <- lung
dat$status <- dat$status - 1   # survival::lung 的 status 为 1=删失 2=死亡 → 转 0/1
cat("n =", nrow(dat), "；事件数 =", sum(dat$status), "；删失 =", sum(dat$status == 0), "\n")

# 实跑捕获（R 4.5.0）
== 步骤 1：数据与清洗 ==
n = 228 ；事件数 = 165 ；删失 = 63

逐行解读：

dat$status <- dat$status - 1：全案例最关键的一行。lung 的原始 status 编码是 1=删失、2=死亡，而 Surv() 的惯例是 0=删失、1=事件——减 1 完成转换。SAS 等价：death = (status = 2); 或 status01 = status - 1;。
拿到任何公开数据先做"事件/删失"清点：165 + 63 = 228，对得上——相当于先看 PROC LIFETEST 日志首页的 censored 汇总，再决定往下走。
不清洗直接建模不会报错，但结果全错且无警告——这类"静默错误"比 SAS 的 log ERROR 更危险，务必养成核对编码的习惯。

SAS 误区：PROC LIFETEST 里你写 time*status(1)，括号里声明的是删失值；而 R 的 Surv(time, status) 默认 status 非零即事件。两边"谁是 1"的语义相反，照搬编码必翻车——先统一成 0/1，再进任何模型。

步骤 2：KM 按性别分层（survfit ↔ PROC LIFETEST）

cat("\n== 步骤 2：KM 按性别分层 ==\n")
fit <- survfit(Surv(time, status) ~ sex, data = dat)
print(fit)

# 实跑捕获（R 4.5.0）
== 步骤 2：KM 按性别分层 ==
Call: survfit(formula = Surv(time, status) ~ sex, data = dat)

        n events median 0.95LCL 0.95UCL
sex=1 138    112    270     212     310
sex=2  90     53    426     348     550

逐行解读：

Surv(time, status) 构造生存对象（时间 + 事件指示），~ sex 声明分层变量——整句对应 proc lifetest data=dat; time time*status(0); strata sex;，公式接口一行搞定。
survfit() 默认 Kaplan-Meier 估计；print 给出每组 n、事件数、中位生存期及 95%CI（默认补数法），即 LIFETEST 的 "Median Survival" 表。
读结果：sex=1（男）中位 270 天（95%CI 212–310），sex=2（女）426 天（348–550）——女性中位生存几乎多出半年，两组 CI 不重叠，差异肉眼可判。

步骤 3：log-rank 检验（survdiff）

cat("\n== 步骤 3：log-rank 检验 ==\n")
print(survdiff(Surv(time, status) ~ sex, data = dat))

# 实跑捕获（R 4.5.0）
== 步骤 3：log-rank 检验 ==
Call:
survdiff(formula = Surv(time, status) ~ sex, data = dat)

        N Observed Expected (O-E)^2/E (O-E)^2/V
sex=1 138      112     91.6      4.55      10.3
sex=2  90       53     73.4      5.68      10.3

 Chisq= 10.3  on 1 degrees of freedom, p= 0.001

逐行解读：

survdiff() 默认 log-rank 检验，对应 PROC LIFETEST "Equality of Survival Curves" 表里的 Log-Rank 行：Chisq=10.3、df=1、p=0.001，组间生存分布差异显著。
Observed/Expected 列可以直接读方向：男性组观测死亡 112 > 期望 91.6，女性组观测 53 < 期望 73.4——男性死亡"超发"。
两组时 (O-E)^2/V 每行都是同一个统计量（10.3），多组时才需要看整体卡方 + 两两比较。
想换检验口径（如 Wilcoxon/Gehan）加参数 rho=1；SAS 侧对应 test w 选项——概念一一对应。

步骤 4：Cox 多变量模型（coxph ↔ PROC PHREG）

cat("\n== 步骤 4：Cox 多变量模型 ==\n")
cox <- coxph(Surv(time, status) ~ age + sex + ph.ecog, data = dat)
print(summary(cox)$coefficients)
cat("HR(sex 每+1单位, 女 vs 男) =", round(exp(coef(cox)["sex"]), 3), "\n")

# 实跑捕获（R 4.5.0）
== 步骤 4：Cox 多变量模型 ==
               coef exp(coef)    se(coef)         z     Pr(>|z|)
age      0.01106676 1.0111282 0.009267411  1.194159 2.324157e-01
sex     -0.55261240 0.5754446 0.167739054 -3.294477 9.860514e-04
ph.ecog  0.46372848 1.5899912 0.113577266  4.082934 4.447067e-05
HR(sex 每+1单位, 女 vs 男) = 0.575

逐行解读：

coxph() 即 PROC PHREG：同一份系数表（coef / exp(coef) / se / z / p），Efron 法处理并列（SAS 默认也是 Efron）。
ph.ecog 的 HR=1.59（p<0.001）：ECOG 体力状态每恶化 1 分，死亡风险约升 59%——临床上最强的预后因子，符合预期。
age 的 HR=1.011、p=0.23：校正性别与 ECOG 后年龄不显著。
方向陷阱：lung 的 sex 是数值编码 1=男、2=女。R 把它当连续变量，exp(coef)=0.575 的真实含义是"sex 每 +1 单位（男→女）"的 HR，即女性死亡风险约为男性的 57.5%。本脚本的 cat 标签已修正为方向自明版本；但审旧代码时常见 "HR(男 vs 女)" 这类含糊标签——SAS 里你会写 class sex; ref='1' 让标签自明，R 里则应 factor(sex, levels=c(1,2), labels=c("Male","Female")) 或至少把方向写进注释。

SAS 误区：在 PHREG 里你习惯了 class 语句自动做哑变量并指定参考水平；coxph 里数值型变量一律按"每增加 1 单位"进模型，多水平分类变量若不 factor 化，得到的 HR 可能根本不是你想要的对比。递交级程序里，分类协变量必须先显式 factor 并声明水平顺序。

步骤 5：KM 图落盘（ggsurvplot）

cat("\n== 步骤 5：KM 图落盘 ==\n")
out_png <- file.path(if (grepl("cases$", getwd())) ".." else ".",
                     "docs", "run-outputs", "c2_km.png")
p <- ggsurvplot(fit, data = dat, pval = TRUE, conf.int = FALSE,
                risk.table = TRUE, palette = c("#7d8f4e", "#b0413e"),
                xlab = "Days", ylab = "Survival probability")
ggsave(out_png, print(p), width = 8, height = 6, dpi = 150)
cat("KM 图已保存:", out_png, "；文件大小 =", file.size(out_png), "bytes\n")
cat("== C2 完成 ==\n")

# 实跑捕获（R 4.5.0）
== 步骤 5：KM 图落盘 ==
KM 图已保存: ./docs/run-outputs/c2_km.png ；文件大小 = 5187 bytes
== C2 完成 ==

实跑产物（KM 曲线 + 底部 risk table，即"在险人数表"）：

按性别分层的 Kaplan-Meier 生存曲线（含 log-rank p 值与 risk table）

逐行解读：

ggsurvplot() 一步拼齐递交级 KM 图：曲线 + log-rank p 值（pval=TRUE）+ 底部在险人数表（risk.table=TRUE）——SAS 里这要用 LIFETEST 的 ODS 图形再加一段 atrisk 数据步拼装。
palette 显式指定两色，图表用色进代码而非手工调 GUI——可复现、可审阅。
ggsave(..., width=8, height=6, dpi=150)：尺寸与分辨率也是规格的一部分，类比 ODS 的 options(orientation/size) 与 RTF 模板设置；产物落 docs/run-outputs/ 便于 QC 归档。

测验：C2 开头为什么必须执行 dat$status <- dat$status - 1？

选 A。Surv(time, status) 把非零 status 一律当"事件"；若不转换，原编码 1=删失、2=死亡两个值都非零，会被解读为"228 例全部发生事件"，删失信息完全丢失，KM 与 Cox 全错且不报任何警告。这是公开数据实战的第一课：先核对事件编码，再进模型。

解读要点

一个包顶两个 PROC：survival = PROC LIFETEST + PROC PHREG 的合体，公式接口（Surv(time, status) ~ covariates）统一贯穿估计、检验与建模。
事件编码是第一道 QC：0=删失、1=事件的惯例与 lung 原始编码相反；任何公开数据/外部数据先 table(status) 清点再动手。
HR 方向必须说清参考：数值型分类变量按"每 +1 单位"进模型，sex 的 HR=0.575 实为女 vs 男；递交程序里先 factor 化并声明水平。
证据三件套：中位生存 270 vs 426、log-rank p=0.001、校正后 Cox 仍显著（HR=0.575，p<0.001）——单变量到多变量结论一致，这正是审评问答里"稳健性"的最小示范。

练习

（易）KM 按 ph.ecog 分层：取 dat2 <- subset(dat, ph.ecog %in% 0:1)，survfit(Surv(time, status) ~ factor(ph.ecog), data = dat2) 并跑 survdiff，比较 ECOG 0 与 1 两组的生存差异，与按性别的结论对照。
（中）Cox 加 age×sex 交互：coxph(Surv(time, status) ~ age * sex + ph.ecog, data = dat)，看交互项 p 值是否显著，并解释交互的含义（age 的效应是否随性别不同而不同）；若显著，写出男/女各自的 age 斜率。
（难）用 flexsurv 做参数生存模型：flexsurvreg(Surv(time, status) ~ sex + age + ph.ecog, data = dat, dist = "weibull")，再换 dist="exp"、"llogis"、"gompertz"，比较各模型 AIC 选优，并与 Cox 半参数结论对照——体会"参数假设换效率"的取舍。

扩展阅读

Project Data Sphere — 注册制公开肿瘤试验数据平台，可申请真实试验级数据练手数据
hbiostat Data Archive — Frank Harrell 维护的教学数据仓库（含 GUSTO、SUPPORT 等经典生存数据）数据

下一案例：离开试验数据，直接连 FDA 的真实世界安全数据库——用 openFDA API 做一次诚实的信号检测。