生存分析Cox回归进阶¶

一句话概述：Cox比例风险模型是生存分析的核心方法，不假设基线风险函数的具体形式，通过风险比(HR)评估协变量对生存时间的影响，广泛用于预后因素分析和预后模型构建。

核心知识点速查表¶

概念	说明
Cox模型	半参数生存分析模型（白话：分析什么因素影响病人活多久）
HR(Hazard Ratio)	风险比，HR>1表示风险增加，HR<1表示保护因素
比例风险假设	Cox模型的核心假设：HR不随时间变化
C-index	一致性指数，衡量模型预测能力（0.5-1，越大越好）
Nomogram	列线图，将Cox模型可视化为临床评分工具
LASSO-Cox	LASSO正则化的Cox回归，用于高维特征选择
时间依赖ROC	评估预后模型在不同时间点的预测能力

一、基础Cox回归¶

# === Cox回归基础 ===
library(survival)
library(survminer)

# 准备数据
data(lung)                                  # 肺癌数据集
head(lung)                                  # 查看数据结构

# 单变量Cox回归
uni_cox <- coxph(                           # Cox回归
  Surv(time, status) ~ age,                 # 生存对象 ~ 协变量
  data = lung                                # 数据
)
summary(uni_cox)                             # 查看结果
# HR(exp(coef)), 95%CI, p-value

# 多变量Cox回归
multi_cox <- coxph(
  Surv(time, status) ~ age + sex + ph.ecog + ph.karno,  # 多个变量
  data = lung
)
summary(multi_cox)

# 森林图可视化
ggforest(multi_cox, data=lung)               # ★森林图（HR和95%CI）

二、比例风险假设检验¶

# === PH假设检验（非常重要！） ===
# Cox模型的前提条件：HR不随时间变化

# Schoenfeld残差检验
ph_test <- cox.zph(multi_cox)               # 检验PH假设
print(ph_test)                               # 查看结果
# p > 0.05 → 满足PH假设 ✓
# p < 0.05 → 违反PH假设 ✗ → 需要处理

# 可视化Schoenfeld残差
ggcoxzph(ph_test)                            # 图形检验（应为水平线）

# === 违反PH假设的处理方法 ===
# 方法1: 分层(Stratification)
strat_cox <- coxph(
  Surv(time, status) ~ age + strata(sex),   # 对sex分层
  data = lung
)

# 方法2: 时间依赖协变量
# 添加 变量*时间 的交互项
td_cox <- coxph(
  Surv(time, status) ~ age + sex + tt(sex),  # tt()=时间依赖
  data = lung,
  tt = function(x, t, ...) x * log(t)        # 时间变换函数
)

三、LASSO-Cox特征选择¶

# === LASSO-Cox：从高维数据中选择预后基因 ===
library(glmnet)

# 准备数据（假设有基因表达矩阵）
x <- as.matrix(expr_matrix)                  # 基因表达矩阵（n x p）
y <- Surv(clinical$time, clinical$status)     # 生存数据

# 交叉验证选择lambda
set.seed(42)
cv_fit <- cv.glmnet(
  x, y,
  family = "cox",                             # Cox模型
  alpha = 1,                                  # alpha=1是LASSO
  nfolds = 10                                 # 10折交叉验证
)

# 查看最优lambda
plot(cv_fit)                                   # lambda vs deviance图
lambda_min <- cv_fit$lambda.min                # 最小deviance的lambda
lambda_1se <- cv_fit$lambda.1se                # 1se规则的lambda（更简约）

# 提取选中的基因
coef_selected <- coef(cv_fit, s="lambda.1se")  # 获取系数
selected_genes <- rownames(coef_selected)[coef_selected[,1] != 0]  # 非零系数的基因
print(paste("选中", length(selected_genes), "个基因"))
print(selected_genes)

# 用选中的基因建立Cox模型
final_cox <- coxph(
  Surv(time, status) ~ .,                     # 所有选中基因
  data = data.frame(time=clinical$time,
                    status=clinical$status,
                    x[, selected_genes])
)
summary(final_cox)

四、列线图(Nomogram)¶

# === 构建列线图 ===
library(rms)

# 设置数据分布
dd <- datadist(lung)                          # 计算变量分布
options(datadist='dd')                         # 全局设置

# 拟合Cox模型（用rms的cph函数）
cox_rms <- cph(
  Surv(time, status) ~ age + sex + ph.ecog,
  data = lung,
  x = TRUE, y = TRUE, surv = TRUE
)

# 构建列线图
surv <- Survival(cox_rms)                     # 生存函数
nom <- nomogram(
  cox_rms,
  fun = list(function(x) surv(365, x),        # 1年生存率
             function(x) surv(730, x)),        # 2年生存率
  funlabel = c("1-Year Survival",
               "2-Year Survival"),
  lp = FALSE                                   # 不显示线性预测值
)
plot(nom)                                       # ★绘制列线图

# 校准曲线
cal <- calibrate(cox_rms, u=365, B=200)       # 1年生存率校准
plot(cal)                                       # 绘制校准曲线

五、时间依赖ROC¶

# === 时间依赖ROC曲线 ===
library(timeROC)

# 计算风险评分
risk_score <- predict(multi_cox, type="risk")  # 预测风险值

# 时间依赖ROC
roc <- timeROC(
  T = lung$time,                               # 生存时间
  delta = lung$status - 1,                     # 事件状态（0/1）
  marker = risk_score,                          # 风险评分
  cause = 1,                                    # 感兴趣的事件
  times = c(365, 730, 1095),                   # 1年、2年、3年
  iid = TRUE                                    # 计算置信区间
)

# 查看AUC
print(roc$AUC)                                 # 各时间点AUC

# 可视化
plot(roc, time=365, col="red", title=FALSE)     # 1年ROC
plot(roc, time=730, col="blue", add=TRUE)       # 2年ROC
legend("bottomright",
       legend=c(paste0("1-Year AUC=", round(roc$AUC[1],3)),
                paste0("2-Year AUC=", round(roc$AUC[2],3))),
       col=c("red","blue"), lty=1)

六、面试高频考点¶

Q1: Cox回归的HR怎么解读？¶

HR=1.5：风险增加50%（该因素每增加1单位，死亡风险增加50%）
HR=0.7：风险降低30%（保护因素）
HR=1.0：无影响
白话：HR就像"倍率"，>1是坏事，<1是好事

Q2: C-index和AUC的区别？¶

C-index衡量整体区分能力，AUC衡量特定时间点的预测能力
C-index=0.7-0.8算不错，>0.8算好
时间依赖ROC比单一C-index信息更丰富

Q3: 什么时候用单变量vs多变量Cox？¶

单变量：初步筛选，看每个变量单独的效果
多变量：控制混杂因素后的独立预后因子
标准流程：单变量筛选(p<0.1) → 多变量验证(p<0.05)

常见报错与解决¶

报错	原因	解决方案
`NA/NaN/Inf in foreign function`	数据有缺失值	先处理NA
`PH assumption violated`	违反比例风险假设	分层或加时间依赖项
`Convergence failure`	变量太多或共线性	减少变量或用LASSO筛选
`C-index太低`	模型预测力差	加入更多变量或换模型

速查表¶

# === Cox回归速查 ===
library(survival); library(survminer)

# 单变量
coxph(Surv(time, status) ~ var, data=df)

# 多变量
fit <- coxph(Surv(time, status) ~ var1+var2+var3, data=df)
ggforest(fit)                        # 森林图

# PH假设检验
cox.zph(fit)                         # p>0.05=满足

# LASSO-Cox
library(glmnet)
cv.glmnet(x, Surv(time,status), family="cox", alpha=1)

# 列线图
library(rms)
nomogram(cph_model, fun=list(...))

# 时间依赖ROC
library(timeROC)
timeROC(T, delta, marker, times=c(365,730,1095))