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

> 在医学研究中，我们常常构建回归模型来探寻风险因素与疾病结局的关联。一个根深蒂固的假设是：这种关联是线性的，如年龄每增加一岁，风险就固定增加多少。但人体非常复杂，真实世界往往并非如此。当连续变量与结局的关系是一条蜿蜒的曲线时，传统的线性模型便显得力不从心。此时，一种强大的可视化与分析方法——限制性立方样条图（RCS）正日益成为顶级期刊的“常客”，一张图就可以揭示隐藏的J型或U型关联，让医学研究从大概相关走向精确刻画。 ![image](https://upload-images.jianshu.io/upload_images/20783724-ee41a17398f677a6.png) --- # 一、何为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) ``` 本文由[mdnice](https://mdnice.com/?platform=6)多平台发布

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

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

相关阅读更多精彩内容

友情链接更多精彩内容