R语言实战：RCS限制性立方条图分析

最新推荐文章于 2026-05-22 11:23:31 发布

原创最新推荐文章于 2026-05-22 11:23:31 发布 · 1k 阅读

8 ·

本内容遵循CC 4.0 BY-SA版权协议

GEO检测

在医学研究中，我们常常构建回归模型来探寻风险因素与疾病结局的关联。一个根深蒂固的假设是：这种关联是线性的，如年龄每增加一岁，风险就固定增加多少。

但人体非常复杂，真实世界往往并非如此。当连续变量与结局的关系是一条蜿蜒的曲线时，传统的线性模型便显得力不从心。

此时，一种强大的可视化与分析方法——限制性立方样条图（RCS）正日益成为顶级期刊的“常客”，一张图就可以揭示隐藏的J型或U型关联，让医学研究从大概相关走向精确刻画。

一、何为RCS分析

1.传统困境

在探究诸如年龄、血压、BMI等连续变量对疾病风险的影响时，研究者常面临两难选择。

常见做法是将连续变量人为分类，例如将年龄分为“青年、中年、老年”。然而，分类的切点选择往往带有主观性，且粗暴的分类会损失大量宝贵信息，可能导致错误结论。
另一种方法是强行拟合线性关系，但这可能完全扭曲了事实。例如，某种药物的疗效可能在某个剂量范围内最佳，过低或过高反而有害，这种“U型”或“J型”关系是线性模型无法捕捉的。

2.何为RCS

RCS是一种灵活的非参数拟合方法，用几段光滑的曲线（通常是3-5段）来拼接出自变量与因变量之间的真实关系。

核心思想是“分段拟合、平滑连接”。在自变量的不同区间内，用不同的多项式函数来拟合，并确保在连接点（称为“节点”）处曲线是平滑的。
所谓限制性，是指对曲线两端施加线性限制，避免在数据范围之外出现不合理的剧烈波动，使结果更稳健。

普通立方条样会让曲线在数据的两端（即 X 变量的最小值、最大值附近）出现无意义的剧烈波动（因为数据少，拟合容易失控）

而RCS强制要求在 X 变量的最小值点（左端点）和最大值点（右端点）处，曲线的二阶导数为 0。这个约束的效果是：让曲线在数据范围的两端保持线性趋势（不再剧烈弯曲），避免了无意义的波动，更贴合实际的生物学 / 统计学规律。

3.RCS的优势

RCS的核心优势在于其客观与精准。它不再依赖研究者的主观分组，而是让数据自己说话，描绘出最可能的关联形态。这使得它成为探索“剂量-反应关系”的利器。

例如，2018年发表在《柳叶刀-糖尿病与内分泌学》上的一项针对360万英国成年人的研究，正是运用RCS清晰揭示了BMI与全因死亡率之间显著的J型关系，找到了死亡风险最低的BMI拐点（约25 kg/m²）。

RCS的应用范围极广，不仅限于生存分析。无论是线性回归、Logistic回归，还是Meta分析中的剂量反应关系，只要想更精细地描述连续变量与结局的关系，都可以引入RCS。

4.R语言实现RCS分析

在R语言中，rms包是完成RCS分析的主力工具。分析流程通常包括几个关键步骤：

首先使用cph()函数（用于Cox回归）或lrm()函数（用于Logistic回归）拟合模型，在模型公式中通过rcs(自变量, 节点数)来指定对某个变量进行RCS拟合。
节点数的选择是关键：通常3-5个节点能在曲线平滑度和避免过拟合之间取得良好平衡。
模型拟合后，使用Predict()函数计算预测值及其置信区间，最后用ggplot2包进行可视化，得到一张带有平滑曲线和置信带的RCS图。

基于Cox回归的RCS代码如下：

#加载包
library(rms)
library(survival)
library(ggplot2)
#提取变量
mydata <- data[, c("stat", "CKM_stage_34","time", "X", "ragender", "age", "hrural", "raeducl", "marry", "Hibpe_treat", "Smoker", "Drinking", "Lip_treat", "Diabetes_treat")]
#转为因子
factor_vars <- c('ragender', 'hrural', 'raeducl', 'marry', 'Hibpe_treat', 'Smoker', 'Drinking', 'Lip_treat', 'Diabetes_treat',"age")
mydata[factor_vars] <- lapply(mydata[factor_vars], as.factor)
#设置datadist（rms包要求）
dd <- datadist(mydata)
options(datadist = 'dd')
#拟合COX回归RCS模型
fit <- cph(Surv(time, stat) ~ rcs(X, 3) + ragender + age + hrural + raeducl + marry + Hibpe_treat + Smoker + Drinking + Lip_treat + Diabetes_treat,data = mydata, x = TRUE, y = TRUE)
#方差分析提取P值,可以查看P for overall和P for nonlinear分别是多少，待会要标记进图里面
anova_result <- anova(fit)
print(anova_result)
#设置参考点（一般默认中位数）
ref_value <- median(mydata$X, na.rm = TRUE)
dd$limits$X[2] <- ref_value
options(datadist = 'dd')
fit <- update(fit)
#生成预测数据
pred_detailed <- Predict(fit, X = seq(min(mydata$X, na.rm = TRUE),max(mydata$X, na.rm = TRUE),length.out = 500),fun = exp, ref.zero = TRUE)
#查找HR=1对应的X值，找到最接近HR=1的点
pred_detailed$diff_to_one <- abs(pred_detailed$yhat - 1)
hr_eq_one_points <- pred_detailed[pred_detailed$diff_to_one < 0.01, ]
# 输出结果
cat("\n=== HR=1对应的X值 ===\n")
cat(sprintf("参考点（中位数）: %.3f (HR=1)\n", ref_value))
#生成RCS图
p <- ggplot(pred_detailed) +
  geom_line(aes(x = X, y = yhat), color = "red", linewidth = 0.8) +
  geom_ribbon(aes(x = X, ymin = lower, ymax = upper),alpha = 0.3, fill = "red") +
  geom_hline(yintercept = 1, linetype = "dashed", linewidth=0.6,color = "grey40") +
  geom_vline(xintercept = ref_value, linetype = "dashed", color = "blue") +
  theme_classic() +
#在图上标注P值,来源于之前方差分析的结果
  annotate("text", x = max(mydata$X, na.rm = TRUE) * 0.6, y = max(pred_detailed$upper, na.rm = TRUE) * 0.9,label = "P for overall < 0.001\nP for nonlinear = 0.375",size = 4, hjust = 0,fontface = "italic") +
# 标注参考点
  annotate("text", x = 7.9, y = 1.2,label = sprintf("%.2f",ref_value),size = 4, color = "blue") +
  labs(x = "X", y = "Hazard Ratio (95% CI)",subtitle = NULL,caption = NULL)+scale_y_continuous(limits = c(0, max(pred_detailed$upper, na.rm = TRUE) + 0.5),  # y轴范围从0开始expand = c(0, 0)   # 取消y轴两端的空白，让轴起点紧贴0) +
   theme(axis.text.x = element_text(size = 10), #X轴刻度字体大小axis.text.y = element_text(size = 10))   # Y轴刻度字体大小)
print(p)