企业官网建设流程全解析

“生存分析”这个名字听起来很有距离感，实际上它在医学随访、客户流失、设备故障、保险精算里遍地都是，核心就一句话：我们要等多久，才会看到某个事件发生？R语言在这个领域的积累非常深，survival 包有三十多年历史，配套的 survminer 画图又漂亮又省事，所以一直到今天，学术界和工业界做生存分析，R 仍然是首选工具之一。

这篇文章我就从实际使用的角度，把生存分析从数据准备、KM 曲线、Cox 回归到常见坑位全部梳理一遍。适合刚接触生存分析的医学研究生、做用户留存分析的运营分析师，以及任何手里握着“时间 + 是否发生”这类数据的人。

1. 先搞懂生存分析在解决什么问题

1.1 生存时间和删失：这两个词绕不开

生存分析里最普通但最关键的数据结构只有两列：一列是时间，一列是事件是否发生。

举个例子，研究一款药物对晚期肺癌患者的效果，我们追踪每个患者从入组开始到死亡或者到研究结束。时间就是“存活了多少天”，事件就是“有没有死亡”。但问题来了：有人中途搬家联系不上，有人还在接受治疗但研究截止日期到了，这些人并没有发生死亡事件。这时候如果把他们的数据当作“没死”来处理，显然不对，因为人家可能出组第二天就去世了，只是我们没观察到。

这就是生存分析所谓的“删失”，专业一点的英文叫 censoring。删失不代表信息无效，它代表的是“我知道这个人至少活到了某个时点，之后我不确定”。R 的 survival 包里所有建模函数都要求你显式告诉它哪些是事件、哪些是删失，这一步做错了，后面全白搭。

删失还分几种：右删失最常见，意思是只知道事件发生在某个时间点之后；区间删失和左删失在医学、金融数据里也存在，但入门阶段右删失最需要理解。你只要记住：删失信息不是缺失值，千万不能粗暴删掉，也不能当成“没发生”，它是生存分析的核心信息之一。

1.2 风险函数：真正的主角

生存分析经常被误以为在分析“生存时间本身”，其实更严谨地说，我们关心的是风险函数 h(t)，它表示“活到时间 t 的那一刻，之后一瞬间发生事件的概率”。

这个概念像开车时的瞬时速度：生存曲线 S(t) 是累计的“到目前为止还活着多少人”，风险函数是“现在这一刻容不容易出事”。很多结论的直观解释都围绕风险展开，比如“吸烟组的死亡风险是不吸烟组的 2 倍”，这里的“风险”就是风险函数之比。

R 语言做生存分析时，你会在输出里看到 HR、exp(coef) 这类词，它们全都是在风险函数的基础上推出来的。理解“风险是一个瞬时概念”，比死记硬背公式更能帮助你读懂结果。

1.3 R 语言生态：survival 包为什么是首选

R 里的生存分析包不少，但我个人建议新手不要东张西望，直接围绕两个包展开就够了：survival 和 survminer。

• survival 是核心计算包，由 Terry Therneau 主导开发，KM 曲线、log-rank 检验、Cox 回归、竞争风险模型都靠它；
• survminer 是基于 ggplot2 的可视化包，专门用来画发表级别的生存曲线、森林图，不用自己手动调一堆坐标轴参数。

安装很简单，就两行：

install.packages("survival") install.packages("survminer")

有的环境里 survminer 依赖的包比较多，装的时候如果报错，多半是缺了 ggplot2、ggpubr 这类依赖包，顺手一起装上就行。

2. 用一张KM曲线快速上手

2.1 先找一份能跑通的示例数据

学生存分析最怕的是手边没有合适的数据。好在 survival 包自带一个非常经典的数据集：lung，收集的是晚期肺癌患者的生存数据。这是公众领域数据，拿来做案例非常合适，你不用自己去模拟一堆假数据。

看下数据结构：

library(survival) data(lung) str(lung)

这个数据有 228 行，每一行代表一位患者。常用的列有：

• time：从入组到死亡或删失的天数；
• status：1 表示删失，2 表示死亡；
• sex：1 表示男性，2 表示女性；
• age：年龄；
• ph.ecog：ECOG 体能评分，0 到 3，分数越高身体状态越差。

最关键的是 status 这一列，它不是常见的 0/1，而是 1/2。很多新手第一次跑就在这里翻车，后面我会专门讲。

2.2 三步画出 KM 生存曲线

KM 曲线，全称 Kaplan-Meier 生存曲线，是最基础的生存率可视化方法。它的思路很朴素：在每个事件发生的时间点上，先计算“活过这个点的概率”，然后把所有点上的条件概率连乘起来，得到一个阶梯状的生存率曲线。

在 R 里画这条曲线分三步。

第一步，把“时间 + 事件状态”包装成 survival 包认识的 Surv 对象：

lung$status2 <- ifelse(lung$status == 2, 1, 0) Surv(lung$time, lung$status2)

Surv 的第一个参数是时间，第二个参数是事件指示变量，一般 0 表示删失，1 表示事件发生。如果你手头数据的 status 正好是 0/1 编码，可以不用转换，直接写 Surv(time, status)。

第二步，用 survfit 拟合 KM 估计。这里加了一个分组变量 sex，看看男性和女性的生存曲线有没有差异：

km_fit <- survfit(Surv(time, status2) ~ sex, data = lung)

第三步，用 survminer 画图：

library(survminer) ggsurvplot( km_fit, data = lung, conf.int = TRUE, risk.table = TRUE, palette = c("#E74C3C", "#3498DB"), xlab = "Time (days)", ylab = "Overall Survival Probability" )

这一步出来的图就是你在论文里常见的那种带置信区间、带风险人数表的生存曲线。

2.3 KM 曲线怎么看

曲线本身的信息量很大，但新手最喜欢盯着 p 值看，反而忽略了几个关键细节。

第一，曲线的台阶处代表有事件发生。每次有人死亡，曲线就往下掉一截；如果曲线上只有竖线没有台阶，说明那个时间点发生了删失，患者失去随访，但曲线本身不下降。

第二，中位生存时间。很多人以为中位生存时间就是“一半患者死亡时对应的时间”，这句话基本对，但因为删失的存在，严格来说它是“生存率降到 50% 时对应的时间点”。R 里可以直接算：

km_fit summary(km_fit)

输出结果里的 median 就是中位生存时间。如果曲线的最后一部分一直没降到 50%，R 会显示 NA，这说明随访时间不够长，没法估计中位生存期。

第三，置信区间。置信区间越宽，说明该时间点上可用的样本越少，估计越不稳定。曲线尾部通常都有一个“喇叭口”，那是样本量越来越少导致的，看到这种图别慌。

2.4 log-rank 检验：两组差异到底有没有意义

看图只能说“男的曲线好像比女的高”，到底有没有统计学差异，需要做检验。最常用的是 log-rank 检验，它的原理是在每个事件时间点上，比较“观察到的死亡人数”和“如果两组没有差异时预期的死亡人数”，最后汇总成一个卡方统计量。

R 里一行代码：

survdiff(Surv(time, status2) ~ sex, data = lung)

输出里的 p 值如果小于 0.05，就认为两组生存曲线有显著差异。对 lung 这个数据来说，通常能观察到性别之间存在显著差异，女性的生存情况好于男性，这也符合很多肺癌临床研究的结论。

需要注意的是，log-rank 检验对两组生存曲线“后期交叉”的情况不敏感。如果两条曲线先在一起、后来分开、最后又交叉，log-rank 检验可能给不出显著结果，但那不代表没有差异。遇到这种情况，可以考虑更灵活的检验方法，比如 Peto 检验或限制平均生存时间，不过这是进阶话题了。

3. Cox 回归：从单因素到多因素

3.1 有了KM曲线为什么还要做Cox

KM 曲线和 log-rank 检验解决的是“分组比较”的问题，但它们有两个明显短板。

第一个短板是连续变量。你想看年龄对生存有没有影响，总不能把年龄切成若干个年龄段画一堆曲线，那样既损失信息又麻烦。第二个短板是混杂因素。两组患者的年龄、体能状态可能本来就不均衡，这时候观察到的差异，到底是性别带来的，还是其他因素带来的？KM 曲线回答不了这种“多因素校正”的问题。

Cox 比例风险回归就是来解决这两个短板的。它能把年龄、性别、体能评分、治疗方式等多个变量同时放进模型，给出每个变量对风险的独立贡献，这也是临床论文中多因素分析最常用的方法。

Cox 模型的核心公式长这样：

h(t) = h0(t) * exp(b1x1 + b2x2 + ... + bp*xp)

h0(t) 是基线风险函数，它随时间是变化的，但 Cox 模型一个巧妙的地方在于：我们根本不用估计 h0(t) 长什么样，只需要估计变量前面的系数 b。这也就是为什么 Cox 回归属“半参数模型”——它不假设基线风险的具体分布，只对变量的效应做参数假设。

3.2 Cox 模型结果怎么读：HR 和置信区间

用 lung 数据跑一个包含年龄、性别、体能评分的 Cox 模型：

cox_model <- coxph( Surv(time, status2) ~ age + factor(sex) + ph.ecog, data = lung ) summary(cox_model)

输出结果会有一大堆，你先把注意力放在这几列上：

• coef：回归系数，也就是公式里的 b；
• exp(coef)：风险比 HR（Hazard Ratio）；
• Pr(>|z|)：p 值；
• 后面的 95% 置信区间。

HR 怎么解释？如果某个变量的 HR 是 1.30，意思是这个变量每增加一个单位，事件发生风险提高 30%。如果 HR 小于 1，说明是保护因素，事件风险下降。

举个例子，我用这个模型跑过的结果里，sex 的 HR 大概在 0.58 左右，这意味着女性患者的死亡风险大约是男性的 58%，也就是低了 42%。这里要注意，因素的水平设置很重要，R 默认会把第一个出现的水平当成参考组。sex 这个变量如果用数值直接丢进模型，参考组就是 1（男性），HR 的符号完全取决于数据的排序，所以最好转成因子并显式指定水平。

lung$sex <- factor(lung$sex, levels = c(1, 2), labels = c("Male", "Female"))

这样在输出里可以明确看到 Male 是参考组，结果解释不会拧巴。

ph.ecog 这个变量也要留意，它默认按连续变量处理，HR 大概是 1.3 左右，意思是 ECOG 评分每升高 1 分，死亡风险增加约 30%。如果你觉得线性效应太强，也可以把它转成因子，看不同评分之间的阶梯效应，这在临床分析中很常见。

3.3 PH 假设检验到底在查什么

Cox 模型有个关键前提，叫“比例风险假设”，英文是 Proportional Hazards Assumption，简称 PH 假设。意思是：不同组别的风险函数比值应该不随时间变化。换句话说，如果性别带来的死亡风险始终是男性比女性高 1.7 倍，那么它就满足 PH 假设；如果第一年男性风险很高，两年后男性和女性风险没差别，那就不满足。

R 里检验 PH 假设非常方便：

cox_zph <- cox.zph(cox_model) print(cox_zph)

输出结果里每个变量有一行，最后还有一个 GLOBAL 综合检验。如果 p 值小于 0.05，通常认为该变量不满足 PH 假设。

遇到不满足 PH 假设的情况，处理方式有几种。

第一，最常用的是把不满足的变量作为分层变量，用 strata() 放进模型，比如：

coxph(Surv(time, status2) ~ age + strata(sex) + ph.ecog, data = lung)

这样每个性别有自己单独的基线风险函数，但其他变量的效应仍然是同一个，模型照样能解释。

第二，如果问题出在某个连续变量上，可以考虑引入时间交互项，用 tt() 处理。但这个方式解释起来比较绕，实际临床文章中并不常用。

第三，如果两组曲线明确交叉，那意味着更复杂的情况，可能要考虑分段模型或其他方法，这时候建议咨询统计师。

3.4 用森林图把结果摆出来

模型跑完，结果怎么展示？survminer 里自带一个 ggforest 函数，可以直接把 Cox 模型的变量、HR、置信区间、p 值画成森林图：

ggforest(cox_model)

森林图的中间部分是每个变量的 HR 点和置信区间横线，横线横跨 1 就说明效应不显著。右边的 p 值也一目了然，临床审稿人非常喜欢这种图。

如果后续要出论文级的图，建议再微调一下字体、配色和标签名称,不过这些细节不影响分析逻辑，先把模型结果读对再说。

4. 实战中踩过的坑

4.1 status 编码颠倒：最隐蔽的坑

我在帮人复现分析时，见过太多人把事件状态搞反。有的数据集里 0 表示“删失”，1 表示“事件”，有的反过来 1 表示活着、2 表示死亡，还有用字符串 "Yes"/"No" 的。

如果你直接把 status 喂给 Surv()，事件指示符不是 0/1，R 大多数情况下会报错或者给出奇怪的警告，但有时候它也会默默把非 0 值当作事件处理，这个更危险。比如 lung 数据里 status 是 1 和 2，如果你直接写 Surv(time, status)，你会发现所有患者都被当成了死亡，log-rank 检验和 Cox 模型的样本量看起来都对，但所有结果就都错了。

我的习惯是，任何数据集拿到手，先做三件事：

table(lung$status) class(lung$time) head(lung)

先看事件状态的频数分布，再确认时间变量是不是数值型，最后扫一眼数据结构。这三步非常快，但能拦下 90% 的低级错误。

4.2 时间变量的格式问题

时间变量是另一个重灾区。很多原始数据表里的时间是日期格式，比如 “2021-03-15”，而生存分析需要的是“从起点到终点的时间间隔”，不是具体日期。

有人直接把日期放进 Surv()，结果系统报错；有人把日期转成字符串再转成数值，那算出来的“时间”就完全不可解释。正确做法是先算两个日期之间的差值，再统一单位。

start_date <- as.Date("2021-03-15") end_date <- as.Date("2022-09-10") follow_up_days <- as.numeric(difftime(end_date, start_date, units = "days"))

换算成月还是天，要结合领域习惯。临床生存分析常用“月”或“天”，用户流失分析常用“天”或“周”。单位本身不影响 Cox 模型的有效性，但会影响 HR 和生存时间的数值大小，写文章时一定要统一并说明。

还有一个小坑：如果患者入组当天就发生了事件，生存时间是 0 天。部分方法对 time=0 的数据非常敏感，处理不当会直接报错。建议看看有多少这样的样本，如果很少可以和统计专家讨论是否排除，或者在敏感性分析里单独验证。

4.3 不满足 PH 假设怎么办

前面已经简单提过分层，这里再补充一个实际经验：很多刚接触的人一看到 cox.zph 的 p 值小于 0.05，就急着把整个模型推翻，其实不需要。

PH 假设是针对某个变量的，不是针对整个模型。检验结果里通常有两类情况：一类是个别变量有问题，另一类是整体有问题。如果只是某个分类变量不满足，分层是最直接的解决办法。如果是连续变量有问题，可以考虑把它离散化后分层，或者换成时间依赖模型。

另外，大样本很容易让 PH 检验变得敏感，微小偏离也会出显著 p 值。这时候不要机械地只看 p 值，可以画一个 scale Schoenfeld 残差图，看残差是不是随机分布在水平线附近，如果整体趋势平稳，就算 p 值临界，也通常可以接受。

4.4 竞争风险：多结局数据要注意

生存分析里还有一个常见情况：你想研究的结局是“因特定疾病死亡”，但有些患者在观察期内可能因为车祸、其他疾病等原因死亡。这些“其他原因死亡”会阻碍目标事件的发生，就叫竞争风险事件。

比如研究肿瘤复发，如果患者因为其他原因死亡了，那就观察不到复发了。这时候如果直接用普通 Cox 模型，会高估累积发生率。稍微专业一点的回答是换用 Fine-Gray 模型，R 里可以用 cmprsk 包的 crr() 函数实现。

但竞争风险模型不是所有场景都要用。如果你的结局是“全因死亡”，也就是不管什么原因死亡都算事件，那就不存在竞争风险问题。所以先想清楚结局定义，再决定要不要换模型。

4.5 给新手的核对清单

把思路整理成一套可以直接对照的流程，分享给大家：

• 确认事件状态编码：0/1 还是 1/2，明确哪个是事件，哪个是删失；
• 确认时间变量类型：是数值型的时间间隔，还是需要换算的日期；
• 确认单位统一：天、月、周，和论文其他部分一致；
• 确认分类变量的因子水平：谁是参考组，解释 HR 时心里要有数；
• 先画 KM 曲线，看数据基本情况，再跑 Cox 回归；
• 跑完模型用 cox.zph 检查 PH 假设，画森林图看结果；
• 如果结局是多原因的，先判断是否需要竞争风险模型。

这套流程看着简单，但每一步都踩过不少人的坑。尤其是事件状态编码，几乎每个拿着自定义数据来找我的人，都要在这上面卡一会儿。我个人实操的体会是：生存分析真正难的不是 R 代码，而是搞清楚每一个数据列在你研究问题里到底是什么意思。代码只是把你想清楚的事情严谨地表达出来而已。先把数据定义弄明白，再动手跑模型，你会少走很多弯路。