如何在R中用survival包绘制带p值标注的生存曲线？一步提升论文可视化水平

第一章：R中survival包生存分析概述

R语言中的`survival`包是进行生存分析的核心工具之一，广泛应用于医学研究、工程可靠性及社会科学等领域。该包提供了构建生存对象、拟合Kaplan-Meier曲线、Cox比例风险模型以及进行回归分析的完整方法体系。

核心功能与数据结构

`survival`包依赖于`Surv`对象来定义事件时间与状态变量。`Surv`函数可创建右删失、左删失或区间删失的生存对象，是后续建模的基础。

• time：表示事件发生的时间点
• event：指示事件是否发生（1=事件发生，0=删失）
• Surv(time, event)：构造基本的生存响应变量

安装与加载

在使用前需确保已安装并加载该包：

# 安装 survival 包
install.packages("survival")

# 加载包
library(survival)

典型应用场景示例

以肺癌数据集`lung`为例，展示如何快速生成Kaplan-Meier估计：

# 构建生存对象并拟合KM模型
fit <- survfit(Surv(time, status) ~ 1, data = lung)

# 查看结果摘要
summary(fit)

上述代码中，survfit函数用于估计生存曲线，~ 1表示整体人群的无分组模型。输出包含中位生存时间、风险表等关键统计量。

函数名	用途说明
Surv()	创建生存响应对象
survfit()	拟合生存曲线（如Kaplan-Meier）
coxph()	拟合Cox比例风险模型

graph TD A[原始数据] --> B[构建Surv对象] B --> C[选择分析方法] C --> D[survfit 或 coxph] D --> E[可视化与推断]

第二章：生存分析基础与数据准备

2.1 生存分析核心概念与Kaplan-Meier估计

生存分析用于研究事件发生时间的统计方法，广泛应用于医学、工程等领域。其核心概念包括**生存时间**、**删失数据**（censoring）和**生存函数** $S(t)$，即个体存活超过时间 $t$ 的概率。

Kaplan-Meier估计器

该非参数方法估算生存函数，适用于右删失数据。其公式为： $$ \hat{S}(t) = \prod_{t_i \leq t} \left(1 - \frac{d_i}{n_i}\right) $$ 其中 $d_i$ 为时间 $t_i$ 处的事件数，$n_i$ 为处于风险中的个体数。

library(survival)
fit <- Surv(time = lung$time, event = lung$status == 2)
km_fit <- survfit(fit ~ 1, data = lung)
plot(km_fit, xlab = "Time (days)", ylab = "Survival Probability")

上述R代码使用survival包构建生存对象并拟合Kaplan-Meier曲线。Surv函数定义时间与事件状态，survfit计算估计值，最终可视化生存概率随时间变化。

结果解读

• 每一步下降代表一个事件发生
• 删失个体用小竖线标记
• 曲线越平缓，生存率越高

2.2 使用survival包构建Surv对象

在R语言中，`survival`包是进行生存分析的核心工具。构建生存数据的第一步是创建`Surv`对象，它用于定义事件时间与事件状态。

Surv函数的基本语法

Surv(time, time2, event, type = "right", origin = 0)

其中，`time`表示右删失时间或区间起始时间；`event`指示事件是否发生（1=事件发生，0=删失）；`type`可设为"right"、"left"、"interval"等，用于指定删失类型。

常见使用示例

以肺癌数据为例：

library(survival)
Surv(lung$time, lung$status == 2)

此处`status == 2`将死亡事件编码为1，删失为0。该表达式生成一个右删失的Surv对象，作为后续`survfit()`或`coxph()`函数的输入。 通过正确构造`Surv`对象，能够准确表达生存数据结构，为模型分析奠定基础。

2.3 分组变量的处理与数据清洗技巧

在数据分析流程中，分组变量的正确处理是确保建模准确性的关键步骤。首先需识别分类变量中的异常值与不一致编码。

常见数据清洗策略

• 统一文本格式：如将“Male”和“male”标准化为“male”
• 处理缺失类别：使用众数或特定标签（如“Unknown”）填充
• 合并稀有类别：避免过细分组导致模型过拟合

代码示例：Pandas 中的分组清洗

import pandas as pd

# 示例数据
df = pd.DataFrame({'gender': ['Male', 'male', 'Female', None, 'FEMALE']})

# 清洗逻辑
df['gender_clean'] = df['gender'].str.lower().fillna('unknown')
df['gender_clean'] = df['gender_clean'].replace({'female': 'F', 'male': 'M'})

上述代码首先将所有值转为小写以实现格式统一，随后填充缺失值并映射为简洁编码，提升后续分组操作的效率与一致性。

2.4 数据结构检查与缺失值应对策略

在数据预处理阶段，首先需对数据结构进行系统性检查，确认字段类型、数据分布及完整性。可通过 pandas 的 info() 与 describe() 方法快速诊断。

常见缺失值识别方法

• isnull().sum()：统计各列缺失值数量
• dtypes：检查字段数据类型是否合理

缺失值处理策略

import pandas as pd
# 填充数值型变量中位数，分类变量众数
df['age'].fillna(df['age'].median(), inplace=True)
df['category'].fillna(df['category'].mode()[0], inplace=True)

上述代码通过中位数填充数值型缺失，避免均值受异常值干扰；分类变量采用众数填充，保持类别分布一致性。对于缺失比例超过60%的特征，建议结合业务逻辑评估是否剔除。

2.5 实战：从临床数据创建生存数据集

在生存分析中，构建结构化的生存数据集是关键步骤。临床数据通常分散在多个表中，需整合患者基本信息、随访记录和事件状态。

数据准备与字段定义

生存数据集核心包含三要素：随访时间（time）、事件状态（event）和协变量（如年龄、性别）。需从原始表中提取并清洗。

• time：最后一次随访或事件发生的时间（单位：月）
• event：1 表示事件发生（如死亡），0 表示删失
• 协变量：如肿瘤分期、治疗方式等

代码实现示例

# 合并患者基线与随访数据
surv_data <- merge(clinical, follow_up, by = "patient_id")

# 构建生存对象所需字段
surv_data$time <- pmin(surv_data$last_followup, surv_data$death_time, na.rm = TRUE)
surv_data$event <- ifelse(is.na(surv_data$death_time), 0, 1)

该代码段首先合并两个关键数据表，随后计算每个患者的观测时间与事件状态。pmin 函数确保取最早发生的终点或最后一次随访时间，逻辑判断生成删失标志，为后续使用 survival 包建模奠定基础。

第三章：绘制基础生存曲线

3.1 利用survfit函数拟合生存模型

在R语言中，`survfit` 函数是用于拟合生存分析模型的核心工具，常与 `Surv` 对象结合使用，以估计 Kaplan-Meier 生存曲线。

基本语法结构

library(survival)
fit <- survfit(Surv(time, status) ~ group, data = lung)
summary(fit)

该代码中，`Surv(time, status)` 创建生存对象，`time` 表示生存时间，`status` 指示事件是否发生（如死亡），`group` 为分组变量。`survfit` 根据分组计算生存率及其置信区间。

输出结果解析

• records：每组的总观测数
• events：发生的事件总数
• median：中位生存时间
• 0.95LCL/0.95UCL：95% 置信限

可视化生存曲线

可结合 `plot(fit)` 或 `ggsurvplot` 绘制直观的生存曲线，辅助判断组间差异。

3.2 使用plot()和ggsurvplot()绘制曲线

在生存分析中，可视化是理解模型结果的关键步骤。R语言提供了多种绘图工具，其中基础的plot()函数和survminer包中的ggsurvplot()函数尤为常用。

使用plot()绘制基础生存曲线

library(survival)
fit <- survfit(Surv(time, status) ~ sex, data = lung)
plot(fit, xlab = "时间 (天)", ylab = "生存概率", col = c("blue", "red"))
legend("topright", legend = c("性别1", "性别2"), col = c("blue", "red"), lty = 1)

该代码生成按性别分组的生存曲线。Surv()定义生存对象，survfit()拟合模型，plot()绘制曲线。参数xlab和ylab设置坐标轴标签，col指定线条颜色。

使用ggsurvplot()增强可视化效果

library(survminer)
ggsurvplot(fit, data = lung, pval = TRUE, risk.table = TRUE, conf.int = TRUE)

ggsurvplot()基于ggplot2，自动美化图形。参数pval添加显著性p值，risk.table显示风险表，conf.int展示置信区间，提升信息密度与可读性。

3.3 曲线样式定制与多组比较可视化

自定义曲线样式提升可读性

在多组数据对比中，合理设置线条颜色、宽度和样式能显著增强图表表现力。Matplotlib 支持通过 color、linewidth、linestyle 等参数精细控制曲线外观。

import matplotlib.pyplot as plt
plt.plot(x, y1, color='blue', linewidth=2, linestyle='-', label='Group A')
plt.plot(x, y2, color='red',  linewidth=1, linestyle='--', label='Group B')
plt.legend()

上述代码中，color 定义线条颜色，linewidth 控制粗细，linestyle 设置实线或虚线，便于区分不同数据组。

多组数据并列对比

使用图例（label 与 legend()）可清晰标识各曲线含义，确保多组数据在同一坐标系下有效比较，提升分析效率。

第四章：添加统计检验与P值标注

4.1 Log-rank检验原理及其在survival中的实现

Log-rank检验的基本思想

Log-rank检验是一种非参数假设检验方法，用于比较两组或多组生存曲线是否存在显著差异。其核心思想是基于“风险集”在每个事件时间点的期望与实际发生事件数之间的偏差进行加权求和，构造卡方统计量。

在R中使用survival包实现

library(survival)
# 构建生存对象并执行log-rank检验
surv_obj <- Surv(time = lung$time, event = lung$status == 2)
logrank_test <- survdiff(surv_obj ~ lung$sex)
print(logrank_test)

上述代码中，Surv() 创建右删失生存对象，survdiff() 执行Log-rank检验。survdiff 返回卡方统计量与p值，用于判断不同分组间生存分布是否显著不同。

结果解读

输出包含各组观察事件数、期望事件数及统计量，其自由度为组数减一，适用于大样本场景。

4.2 提取并格式化假设检验P值

在统计分析中，p值是判断假设检验显著性的关键指标。准确提取并规范呈现p值，有助于提升结果的可读性与专业性。

从统计模型中提取p值

以Python中的`statsmodels`库为例，线性回归模型返回的结果对象包含详细的统计信息，可通过属性访问p值：

import statsmodels.api as sm
X = sm.add_constant(data['x'])
y = data['y']
model = sm.OLS(y, X).fit()
p_values = model.pvalues

上述代码中，model.pvalues 返回一个包含每个回归系数p值的Series，索引对应变量名。

格式化输出p值

为增强可读性，常将p值四舍五入至小数点后三位，并标记显著性水平：

• < 0.001: ***
• < 0.01: **
• < 0.05: *

变量	P值	显著性
截距	0.002	**
斜率	0.000	***

4.3 在图形中自动标注P值（base R与ggplot2方案）

在统计可视化中，标注显著性P值是展示组间差异的关键步骤。Base R 和 ggplot2 提供了灵活的实现方式。

Base R 中添加P值

使用 text() 或 mtext() 可将P值手动添加到绘图区域：

# 示例：t检验后标注P值
test_result <- t.test(len ~ supp, data = ToothGrowth)
p_value <- round(test_result$p.value, 3)
boxplot(len ~ supp, data = ToothGrowth)
text(1.5, max(ToothGrowth$len), labels = paste("p =", p_value), cex = 0.8)

该方法直接控制文本位置，适合简单图形，但需手动调整坐标以避免重叠。

ggplot2 自动化标注

结合 ggpubr 包可自动完成统计检验与标注：

• stat_compare_means() 支持多种检验方法
• 可自定义标签位置与格式
• 无缝集成于 ggplot 流程

library(ggpubr)
ggboxplot(ToothGrowth, x = "supp", y = "len") +
  stat_compare_means(method = "t.test")

此方案提升可复用性，适用于复杂图表的批量生成。

4.4 综合示例：发表级生存图的完整绘制流程

数据准备与生存对象构建

在绘制生存曲线前，需确保数据包含时间（time）和事件状态（status）变量。使用R语言的survival包构建生存对象。

library(survival)
surv_obj <- Surv(time = lung$time, event = lung$status == 2)

该代码创建一个Surv对象，其中time表示生存时间，status == 2将死亡事件编码为1，其他为0，是生存分析的标准输入格式。

拟合Kaplan-Meier模型并绘图

通过survfit()拟合非参数生存曲线，并按性别分组：

km_fit <- survfit(surv_obj ~ sex, data = lung)

此模型将生成分组的生存估计，用于后续可视化。

发表级图形绘制

使用ggsurvplot来自survminer包，生成高分辨率、符合期刊要求的图形：

library(survminer)
ggsurvplot(km_fit, data = lung, pval = TRUE, risk.table = TRUE, 
           xlab = "Time (days)", ylab = "Survival Probability")

参数pval = TRUE添加对数秩检验p值，risk.table展示各时间点的危险人数，提升图表信息密度。

第五章：提升科研图表表现力的进阶建议

优化视觉层次与色彩对比

在科研图表中，合理的色彩搭配能显著提升信息传达效率。避免使用高饱和度颜色组合，推荐采用 ColorBrewer 提供的科学配色方案。例如，在绘制多类别柱状图时，使用发散型调色板突出对照组与实验组差异。

增强数据可读性

当图表包含多个数据系列时，应通过图例位置、线型与标记形状进行区分。以下代码展示了 Matplotlib 中如何自定义图例并提升字体可读性：

import matplotlib.pyplot as plt

plt.figure(figsize=(8, 5))
plt.plot(x, y1, label='Control', color='#1f77b4', linewidth=2)
plt.plot(x, y2, label='Treatment', color='#d62728', linestyle='--', marker='o', markersize=4)
plt.legend(loc='upper left', frameon=False, fontsize=12)
plt.xlabel('Time (hours)', fontsize=11)
plt.ylabel('Expression Level', fontsize=11)
plt.tick_params(axis='both', which='major', labelsize=10)
plt.show()

选择合适的图表类型

不同类型的数据应匹配最能体现其特征的可视化形式。下表列举了常见科研数据场景与推荐图表：

数据类型	推荐图表	适用场景
时间序列	折线图	基因表达动态变化
分类比较	箱形图	多组间统计分布对比
相关性分析	热图	基因共表达网络

嵌入交互式元素（适用于电子出版）

现代科研论文常以 PDF 或网页形式发布，可嵌入可缩放矢量图形（SVG）或使用 Plotly 生成交互式图表。用户可通过悬停查看具体数值，提升数据探索体验。