生存分析Cox回归¶

一句话概述：生存分析研究"多久会发生某事件"，Kaplan-Meier画生存曲线看差异，Cox回归量化各因素的风险比（HR），是临床生信分析的核心技能。

核心知识点速览¶

概念	白话解释
生存分析	研究"从开始到事件发生经过了多久"的统计方法
生存时间	从起点（如确诊）到终点（如死亡）的时间
删失（Censoring）	观察期内没发生事件，你只知道"至少活了这么久"
Kaplan-Meier曲线	用阶梯形图展示不同组的"存活概率随时间变化"
Log-rank检验	比较两条KM曲线是否有统计学差异
Cox回归	多因素分析，量化每个变量对死亡风险的影响大小
HR（风险比）	HR=2表示风险是对照组的2倍，HR=0.5表示风险减半
比例风险假设	Cox回归的前提：两组的风险比在任何时间点都恒定
森林图	把多个变量的HR和置信区间画在一起的可视化
C-index	模型区分能力指标，0.5=随机，>0.7=好，>0.8=优秀

一、基本概念¶

1.1 生存数据的特殊性¶

普通数据：身高 = 175cm（精确值）
生存数据：存活 36个月后失访（删失值）

删失的意思是：
  患者A：观察了36个月，还没死，但联系不上了
  → 你只知道他至少活了36个月
  → 不能扔掉这个数据（信息浪费）
  → 也不能当成36个月就死了（不准确）
  → 生存分析专门处理这种"不完整"数据

1.2 Surv对象¶

library(survival)  # 加载survival包（核心）

# 创建生存对象（Surv对象）
# time = 观察时间（月/年/天）
# status = 事件是否发生（1=死亡/复发，0=删失/存活）
surv_obj <- Surv(time = patient_data$time, event = patient_data$status)

# 查看Surv对象
head(surv_obj)
# [1]  36+  24  48   12+  60   30
# 数字后面带 "+" 的表示删失（没发生事件）
# 没有 "+" 的表示观察到了事件（如死亡）

二、Kaplan-Meier生存分析¶

2.1 基本KM分析¶

library(survival)   # 生存分析核心包
library(survminer)  # 生存曲线美化绘图

# 读取示例数据
data(lung)  # R自带的肺癌数据集
head(lung)  # 查看数据结构
# inst: 机构 | time: 生存时间 | status: 1=删失,2=死亡
# age: 年龄 | sex: 1=男,2=女 | ph.ecog: ECOG评分

# 修正status：survival包需要 0=删失, 1=事件
lung$status <- lung$status - 1  # 1→0(删失), 2→1(死亡)

# 拟合KM生存曲线（按性别分组）
km_fit <- survfit(
  Surv(time, status) ~ sex,  # 公式：Surv(时间, 状态) ~ 分组变量
  data = lung                 # 数据框
)

# 查看摘要
summary(km_fit)                # 每个时间点的存活率
print(km_fit)                  # 简要信息：中位生存时间等

2.2 绘制KM曲线（发文章质量）¶

# 使用survminer画高质量KM曲线
ggsurvplot(
  km_fit,                           # KM拟合结果
  data = lung,                       # 原始数据
  pval = TRUE,                       # 显示log-rank p值
  pval.coord = c(100, 0.9),          # p值位置
  conf.int = TRUE,                   # 显示95%置信区间
  risk.table = TRUE,                 # 显示风险表（每个时间点的人数）
  risk.table.col = "strata",         # 风险表按组着色
  risk.table.height = 0.25,          # 风险表高度
  linetype = "strata",               # 不同组用不同线型
  palette = c("#E7B800", "#2E9FDF"), # 自定义颜色
  ggtheme = theme_bw(),              # 主题
  legend.title = "Sex",              # 图例标题
  legend.labs = c("Male", "Female"), # 图例标签
  xlab = "Time (days)",              # X轴标签
  ylab = "Survival Probability",     # Y轴标签
  title = "Kaplan-Meier Survival Curve",  # 图标题
  surv.median.line = "hv",           # 中位生存时间辅助线
  break.x.by = 100                   # X轴刻度间隔
)

2.3 Log-rank检验¶

# Log-rank检验：比较两组生存曲线是否有差异
logrank_test <- survdiff(
  Surv(time, status) ~ sex,  # 公式
  data = lung                 # 数据
)
print(logrank_test)  # 查看结果
# Chi-squared p值 < 0.05 说明两组有显著差异

2.4 按基因表达分组¶

# 生信分析中常见：按基因表达高/低分组做KM分析
# 方法1：用中位数分组
lung$gene_group <- ifelse(lung$gene_expr > median(lung$gene_expr),
                           "High", "Low")  # 高于中位数=High

# 方法2：用survminer找最佳切点
res_cut <- surv_cutpoint(
  lung,                    # 数据
  time = "time",           # 时间列
  event = "status",        # 事件列
  variables = "gene_expr"  # 要分组的连续变量
)
summary(res_cut)  # 查看最佳切点

# 按最佳切点分组
res_cat <- surv_categorize(res_cut)  # 自动分组
km_gene <- survfit(Surv(time, status) ~ gene_expr, data = res_cat)

# 画图
ggsurvplot(km_gene, data = res_cat, pval = TRUE, risk.table = TRUE)

三、Cox比例风险回归¶

3.1 单因素Cox回归¶

# 单因素Cox回归（逐个变量分析）
# 分析年龄对生存的影响
cox_age <- coxph(Surv(time, status) ~ age, data = lung)
summary(cox_age)

# 批量做单因素Cox回归
variables <- c("age", "sex", "ph.ecog", "ph.karno", "pat.karno", "wt.loss")

# 循环每个变量
univ_results <- lapply(variables, function(var) {
  formula <- as.formula(paste("Surv(time, status) ~", var))  # 构建公式
  cox_model <- coxph(formula, data = lung)                    # 拟合模型
  s <- summary(cox_model)                                      # 获取摘要

  # 提取关键信息
  data.frame(
    Variable = var,                                  # 变量名
    HR = round(s$conf.int[, "exp(coef)"], 3),        # 风险比
    HR_lower = round(s$conf.int[, "lower .95"], 3),  # HR 95%CI下限
    HR_upper = round(s$conf.int[, "upper .95"], 3),  # HR 95%CI上限
    P_value = round(s$coefficients[, "Pr(>|z|)"], 4) # p值
  )
})

univ_table <- do.call(rbind, univ_results)  # 合并结果
print(univ_table)  # 查看单因素分析表

3.2 多因素Cox回归¶

# 多因素Cox回归（所有显著变量一起分析）
# 通常把单因素p<0.05或p<0.1的变量纳入多因素
cox_multi <- coxph(
  Surv(time, status) ~ age + sex + ph.ecog,  # 多个自变量
  data = lung
)
summary(cox_multi)

# 关键输出解读：
# coef       — 回归系数
# exp(coef)  — HR（风险比），最重要的指标
# se(coef)   — 标准误
# z          — z统计量
# Pr(>|z|)   — p值
# Concordance — C-index（模型区分能力）

3.3 HR的解读¶

HR（Hazard Ratio，风险比）的含义：
  HR = 1.0  → 没有影响
  HR = 2.0  → 风险是对照的2倍（危险因素）
  HR = 0.5  → 风险是对照的一半（保护因素）
  HR = 1.5  → 风险增加50%

95% CI的解读：
  HR = 1.5, 95%CI (1.1-2.0), p=0.01
  → HR=1.5且CI不跨过1.0，说明显著

  HR = 1.3, 95%CI (0.8-2.1), p=0.25
  → CI跨过了1.0，说明不显著

四、森林图（Forest Plot）¶

# 方法1：用survminer画森林图
ggforest(
  cox_multi,                    # 多因素Cox模型
  data = lung,                  # 原始数据
  main = "Hazard Ratio (95% CI)",  # 标题
  fontsize = 1.0                # 字体大小
)

# 方法2：自定义森林图（用forestplot包）
library(forestplot)

# 准备数据
forest_data <- data.frame(
  Variable = c("Age", "Sex (Female)", "ECOG Score"),
  HR = c(1.02, 0.59, 1.59),
  Lower = c(1.00, 0.42, 1.23),
  Upper = c(1.04, 0.83, 2.05),
  P = c(0.042, 0.003, 0.001)
)

# 画图
forestplot(
  labeltext = cbind(
    forest_data$Variable,
    paste0(forest_data$HR, " (", forest_data$Lower, "-", forest_data$Upper, ")"),
    forest_data$P
  ),
  mean = forest_data$HR,         # HR值
  lower = forest_data$Lower,     # CI下限
  upper = forest_data$Upper,     # CI上限
  zero = 1,                       # 参考线位置（HR=1）
  is.summary = c(TRUE, rep(FALSE, nrow(forest_data))),
  col = fpColors(box = "royalblue", line = "darkblue"),
  xlab = "Hazard Ratio"
)

五、比例风险假设检验¶

# Cox回归的前提假设：比例风险（PH假设）
# 即两组的风险比HR在任何时间点都恒定

# 1. Schoenfeld残差检验（标准方法）
ph_test <- cox.zph(cox_multi)  # 检验PH假设
print(ph_test)
# 如果 p > 0.05，说明满足PH假设（好事）
# 如果 p < 0.05，说明违反PH假设（需要处理）

# 2. 可视化Schoenfeld残差
ggcoxzph(ph_test)  # 残差应该是水平线（无时间趋势）

# 3. 如果违反PH假设的处理方法
# 方法a：分时间段分析
# 方法b：使用时间依赖系数
# 方法c：使用分层Cox（strata）
cox_strata <- coxph(
  Surv(time, status) ~ age + strata(sex),  # 按sex分层
  data = lung
)

六、TCGA生存分析实例¶

# 从TCGA数据库下载临床数据并做生存分析
library(TCGAbiolinks)   # TCGA数据下载
library(SummarizedExperiment)

# 下载临床数据
clinical <- GDCquery_clinic(project = "TCGA-LUAD", type = "clinical")

# 准备生存数据
surv_data <- data.frame(
  patient = clinical$submitter_id,
  time = clinical$days_to_last_follow_up,        # 随访时间
  status = ifelse(clinical$vital_status == "Dead", 1, 0),  # 生存状态
  stage = clinical$ajcc_pathologic_stage          # 临床分期
)

# KM分析（按分期分组）
km_stage <- survfit(Surv(time, status) ~ stage, data = surv_data)
ggsurvplot(km_stage, data = surv_data, pval = TRUE, risk.table = TRUE,
           palette = "jco", xlab = "Days")

常见报错与解决¶

报错信息	原因	解决方案
`Error in Surv(): invalid status`	status列不是0/1	确保0=删失，1=事件
`Loglik converged before variable`	完全分离	该变量预测太完美，去掉或合并组别
`Error in coxph.fit: NA/NaN/Inf`	数据有缺失值	用na.omit()或complete.cases()过滤
`cox.zph: p < 0.05`	违反PH假设	用分层Cox或时间依赖变量
`ggsurvplot: strata not found`	公式里没有分组变量	添加 ~group 或 ~1（无分组）
`risk.table显示不正常`	survminer版本问题	更新survminer或检查ggplot版本

速查表¶

# KM分析三步走
survfit(Surv(time, status) ~ group, data) → ggsurvplot() → survdiff()

# Cox回归三步走
coxph(Surv(time, status) ~ var1 + var2, data) → summary() → ggforest()

# PH假设检验
cox.zph(model) → ggcoxzph()

# 关键包
survival    — 核心函数（Surv, survfit, coxph）
survminer   — 美化绑图（ggsurvplot, ggforest, ggcoxzph）
forestplot  — 自定义森林图

# HR解读速记
HR > 1: 危险因素（数值越大风险越高）
HR < 1: 保护因素（数值越小保护越强）
HR = 1: 无影响
95%CI不跨1: 显著 | 跨1: 不显著

面试高频问题¶

Q1：什么是删失（Censoring）？为什么重要？¶

答：删失是指在观察期内没有观察到事件发生。比如研究5年生存率，但某患者只随访了3年就失访了——你只知道他至少活了3年，不知道后来怎样。删失数据很重要，因为：①如果直接丢弃，会浪费大量数据；②如果当成"死亡"处理，会高估死亡率；③KM和Cox回归能正确利用删失信息，给出无偏估计。

Q2：Log-rank检验和Cox回归有什么区别？¶

答：Log-rank检验是非参数方法，只能比较两组的整体生存差异（给一个p值），不能调整混杂因素。Cox回归是半参数方法，可以同时分析多个变量的影响，给出每个变量的HR和p值，还能调整混杂因素。通常先用KM+Log-rank做初步探索，再用Cox回归做多因素分析。

Q3：什么是比例风险假设？违反了怎么办？¶

答：比例风险假设（PH假设）是Cox回归的核心前提，意思是两组的风险比HR在任何时间点都是恒定的（曲线不交叉）。用cox.zph()检验，p>0.05表示满足。如果违反（p<0.05），可以：①用分层Cox（strata）——对违反PH的变量分层；②引入时间依赖系数（time-varying coefficients）；③分时间段分析。

Q4：HR=2.0和OR=2.0有什么区别？¶

答：HR（Hazard Ratio）来自Cox回归，考虑了时间因素，表示"在任何时间点，某组的瞬时风险是对照组的2倍"。OR（Odds Ratio）来自logistic回归，不考虑时间，只比较两组的事件发生比例。对于罕见事件（发生率<10%），HR≈OR；但对于常见事件，两者差异很大。生存分析必须用HR。

Q5：如何从TCGA数据做基因的生存分析？¶

答：①用TCGAbiolinks下载临床数据和表达数据；②提取感兴趣基因的表达值和患者生存信息（time + status）；③按基因表达高低分组（中位数或surv_cutpoint找最佳切点）；④做KM分析+Log-rank检验看两组差异；⑤做多因素Cox回归调整年龄、分期等混杂因素。关键是确保样本ID能正确匹配临床和表达数据。