R语言小白绘图系列

第 35 弹|连续变量生存曲线(最优cutoff)

目录

1. 连续变量生存曲线的背景和应用
2. 选择最优cutoff的理论基础
3. 常见的最优cutoff选择方法解析
4. 使用 R 语言进行最优cutoff的实践操作
5. 如何解释结果
6. 高级话题：考虑其他生存分析方法和模型优化

1 连续变量生存曲线的背景和应用

在生存分析中，连续变量（如基因表达量、临床生物标志物等）通常用于预测预后。然而，直接将连续变量应用于生存分析可能面临挑战，因为常规的Kaplan-Meier生存曲线需要二分类变量。为此，通常我们需要确定一个最佳cutoff，将连续变量分为“高”和“低”两组。这个cutoff值可以最大化两组之间的生存差异，从而揭示该变量对预后的影响。

应用场景
最优cutoff方法广泛应用于肿瘤生物标志物的发现，心血管疾病的预后预测，基因组学和转录组学数据分析等。
通过将复杂的连续变量转化为易于解读的二分类变量，临床医生和研究人员能够更直观地评估特定因子对患者生存的影响。

2. 选择最优cutoff的理论基础

在实际应用中，选择cutoff并不简单。假设我们有一个连续变量，它可以预测患者的生存期，如何选择一个cutoff值，将样本分为两组，以便在生存曲线中显示显著差异？

关键思想：我们需要通过一种策略来找到一个cutoff点，使得分组后的生存差异（如通过log-rank检验的p值或卡方检验统计量）最大化。这里存在以下几个关键挑战：

• 统计学与生物学意义：cutoff选择不仅要具有统计显著性，还应考虑生物学合理性。
• 过拟合：在寻找cutoff时，如果不加以限制，容易选择一个p值最小的cutoff，导致过拟合。
• 多重检验校正：遍历多个cutoff值的过程中，需要控 制多重检验问题，避免出现假阳性结果。

3. 常见的最优cutoff选择方法解析

1. 基于中位数和四分位数的cutoff：

• 中位数：中位数法将连续变量按50%分位数切分为两组，通常是最直接且常用的方法。适用于样本量较小且变量分布均匀的场景。
• 四分位数：可进一步精细化，选择上、下四分位数或其它百分位数来进行多组分类（如低、中、高表达）。
• 优缺点：优点：简单，易解释；缺点：在数据分布非对称或有显著偏态时，可能会产生不准确的分组。

2. 最大χ²法（Maximal Chi-square Method）：

• 通过遍历所有可能的cutoff值，计算每个cutoff值对应的两组之间的生存曲线的差异（通常使用log-rank检验），选择使得检验统计量（如χ²值）最大的那个cutoff。这样，生存曲线之间的差异最大。
• 优缺点：优点：能够从多个潜在cutoff中选择最优值，保证统计差异最大；缺点：易引入多重检验问题，需考虑适当的校正。

3. 最小p值法（Minimum P-value Method）：

• 与最大χ²法类似，此方法选择使得log-rank检验p值最小的cutoff。通常与最大χ²法配合使用。

4. 时间依赖性ROC曲线法（Time-dependent ROC Curve）：

• 基于时间依赖性ROC曲线（如AUC）的cutoff选择方法，通过最大化灵敏度和特异性来确定cutoff。此法常用于生存时间预测模型的评估。

实践中的应用

在应用这些方法时，需要注意以下几点：

• 过拟合：在选择cutoff时，如果只关注p值最小或χ²值最大，可能会导致过拟合。建议在独立验证集上进行验证。
• 多重检验问题：如果通过遍历多个cutoff值来确定最佳cutoff，需考虑多重检验校正，避免假阳性结果。
• 生物学意义：除了统计学上的显著性，还应结合临床或生物学意义选择合理的cutoff。

应用实例

在下图中：

• 横坐标代表生存时间，单位是年；
• 纵坐标代表生存率；
• 根据目标基因的表达水平可以寻找一个最优cutoff，再根据cutoff将样本分为高低表达两组，然后可以比较两组之间生存是否具有差异。
• 通过P值可以了解两组之间生存率是否存在差异，若p<0.05则说明高低表达两组之间的生存率存在差异；
• 如下图，p<0.001,两组之间生存存在差异，说明该基因与生存相关。生存曲线下方表格代表不同时间点高低表达组存活的患者数目。

4. 使用 R 语言进行最优cutoff的实践操作

现在进入R语言中的实际操作部分。

4.1 环境准备

# 安装并加载所需的包
install.packages("survival")
install.packages("survminer")

这两行代码用于安装survival和survminer包。这两个包在R中用于生存分析和生存曲线的可视化。

library(survival)
library(survminer)

这两行代码用于加载之前安装的survival和survminer包。加载这些包后，可以使用其中的函数进行生存分析和绘图。

4.2 数据准备与处理

假设我们有一个包含生存数据的文件 input.txt，格式如下：

输入文件有四列信息

• id: 样品名称
• futime:生存时间，单位为年。
• fustat:生存状态，0代表存活，1代表死亡。
• 基因表达量

setwd("D:\\biowolf\\bioR\\36.survivalCutoff")

这行代码将当前工作目录设置为"D:\\biowolf\\bioR\\36.survivalCutoff"。这样，所有的文件读取和写入操作都会在这个目录下进行。

# 输入文件有四列信息
# id: 样品名称
# futime: 生存时间，单位为年
# fustat: 生存状态，0代表存活，1代表死亡
# 基因表达量: 连续变量 (如MATN3)

inputFile = "input.txt"  # 文件路径
rt = read.table(inputFile, header = TRUE, sep = "\t", check.names = FALSE)

# 选择生存分析所需的列
rt = rt[, c("futime", "fustat", "MATN3")]
colnames(rt) = c("futime", "fustat", "var")

• inputFile 变量设置了输入数据文件的名称为"input.txt"。
• outFile 变量设置了输出文件的名称为"survival.pdf"，这是将保存的生存曲线图。
• var 变量设置为"MATN3"，表示将使用输入数据中的MATN3列进行生存分析。

4.3 计算最优cutoff

我们将使用 surv_cutpoint 函数来计算cutoff。

# 计算最优cutoff
res.cut = surv_cutpoint(rt, time = "futime", event = "fustat", variables = c("var"))
res.cut

输出的res.cut对象包含最优cutoff值及其相关信息。

4.4 数据分类并拟合生存曲线

根据cutoff值将样本分为高、低表达组，进行生存分析：

# 根据cutoff值分类
res.cat = surv_categorize(res.cut)
# 拟合生存曲线
fit = survfit(Surv(futime, fustat) ~ var, data = res.cat)

4.5 计算生存差异并绘制生存曲线

# 计算生存差异的p值
diff = survdiff(Surv(futime, fustat) ~ var, data = res.cat)
pValue = 1 - pchisq(diff$chisq, df = 1)
if (pValue < 0.001) {
  pValue = "p<0.001"
} else {
  pValue = paste0("p=", sprintf("%.3f", pValue))
}

# 绘制生存曲线
surPlot = ggsurvplot(fit, 
                     data = res.cat,
                     conf.int = TRUE,  
                     pval = pValue, 
                     legend.labs = c("High", "Low"),
                     xlab = "Time (years)", 
                     palette = c("red", "blue"),  
                     risk.table = TRUE)
# 保存图像
pdf("survival.pdf")
print(surPlot)
dev.off()

这段代码的目的是使用ggsurvplot函数生成一个包含生存曲线、置信区间、p值和风险表的图，并将其保存为PDF文件。这种可视化方式清晰展示了高低表达组之间的生存差异。

结果解释：

生成的生存曲线将显示高低表达组之间的生存差异，p值反映了组间差异的显著性。如果p值小于0.05，则说明组间存在显著差异。

5. 如何解释结果

生存曲线的可视化展示了变量对患者预后的影响。关键点如下：

• 高/低表达组的生存曲线差异：生存曲线之间的差异越大，说明该变量的预后价值越高。
• p值的意义：p值小于0.05表示高、低表达组在统计学上有显著的生存差异。
• 风险表：风险表展示了在不同时间点上，各组的风险人数，便于解读生存趋势。

6. 高级话题：其他生存分析方法

6.1 Cox回归模型

对于多变量分析，Cox比例风险回归模型是常用方法，可以同时考虑多个因素的影响。

coxph.fit = coxph(Surv(futime, fustat) ~ var + 其它变量, data = rt)
summary(coxph.fit)

6.2 LASSO-Cox回归

如果有多个变量，可以使用LASSO-Cox回归进行变量筛选，确保模型的简化和优化。

6.3 多重检验校正

为了避免在多个cutoff值上选择时出现假阳性结果，使用多重检验校正如Bonferroni校正或FDR校正是必要的。

总结

通过这篇深入的 R 语言教学贴，我们从连续变量的生存分析出发，详细讨论了选择最优 cutoff 的各种方法，并通过代码示例展示了如何在 R 语言中实际操作，从计算 cutoff 到拟合生存曲线。我们探讨了常见的中位数法、最大 χ² 法、最小 p 值法等不同的选择方法，并结合实际分析给出了直观的图形和 p 值解释。同时，我们也讨论了如何使用 Cox 回归等高级生存分析方法来进一步优化和解释数据。

这篇教学贴不仅涵盖了基础操作，还通过高级话题延伸了对生存分析的理解，希望可以帮助读者深入掌握生存分析中的最优 cutoff 选择，并灵活运用于实际研究中。为了避免过拟合和多重检验带来的统计陷阱，文中还提醒了多重检验校正和生物学意义的重要性。这些步骤能够帮助研究者从统计和生物学的角度全面评估变量对生存预后的影响。