BMJ非线性指南：喊你用限制性立方样条rcs!推荐rcssci包
BMJ前2周发表了非线性数据的分析指南，网页https://www.bmj.com/content/390/bmj-2024-082440。文章串联介绍了目前市面上的非线性方法，详细阐述了正确处理连续变量的重要性，并描述了如何使用样条函数或分数多项式来模拟非线性关系，并说明了直接转化为分类入选模型的后果。下面让我们一起深入解析、复现这篇BMJ顶刊文章。
一、 前言
在医学研究中，年龄、生命体征或生物标志物浓度等连续变量非常常见。无论研究目的是描述（例如，年龄是否与确诊 COVID-19 后六个月的死亡率相关）、解释（例如，新抗癌药物的效果是否因连续生物标志物的值而异）还是预测（例如，将血压加入模型是否能提高心血管疾病风险预测的准确性），自变量与因变量之间的关联入选的形式十分重要。 在建模连续自变量与因变量之间的关联时，将x直接分类或直接假设存在线性关系，可能导致错误的结果！因为非线性关系在自然界中很常见，需要灵活的建模技术才能充分捕捉 样条函数和（分数）多项式是两种广泛使用的工具，用于允许自变量和因变量之间的非线性关系。
二、错误形态后果
线性化简化了真实的复杂性，分类化则粗暴切割了信息，两者都可能严重影响统计推断和预测质量。 线性建模的局限性：线性关系假设连续自变量每增加一个单位，因变量的效应是恒定的，例如年龄 30→31 岁与 90→91 岁的效应被认为相同。虽然这种方法建模简单、解释直观，也是默认选择，但在自然界真实的非线性关系中往往不成立。若强行使用线性假设，可能导致误导性的结论、预测性能下降，并降低临床研究中的统计功效。 分类建模的弊端：把连续变量分组看似直观，但会造成信息丢失、统计效能下降、预测性能受损，还引入不合理的阶跃函数。分组假设区间内效应恒定，且截点往往任意，导致结果难以比较和复现。尤其是二分法信息损失最大，而所谓“最佳截断点”方法更会增加假阳性和不可重复性。总体上，分类未能充分利用数据，常带来残余混杂，类别越少问题越严重，多分组只能部分缓解。
三、正确形态
数据从均匀分布中抽样，并拟合了四种不同的关系（常数、线性、二次和非单调）（黄线，即真实关系）。对于每种关系，都显示了一个线性函数（y～x）、一个使用受限三次样条的非线性函数（y～rcs(x,3)）和一个使用分数项多项式的非线性函数（y～fp(x)），以供比较。其他三条线越接近黄线，拟合度越好（即，用于拟合数据的函数更准确地描述了真实的潜在关系）。以非线性方式对常数或线性关系（前两幅图）进行建模只会导致由于过度拟合而产生的轻微误差。但是，当真实关系是非线性的（后两幅图）时拟合线性模型可能会导致严重错误。
3.1 不同领域线性误用
3.2 多个非线性x应同时入组模型
3.3 案例1：线性vs非线性
临床研究和预测建模中对连续变量的处理长期存在问题，主要表现为过度分类和不加检验的线性假设，导致信息损失、功效下降和结论误导。在预测建模、预后研究和随机对照试验中，连续变量经常被随意分类或假设为线性，这种做法造成信息丢失、混杂未控和功效下降，系统评价已反复证明其弊端，迫切需要用平滑的非线性方法更科学地处理连续变量。见图1，rcs的优点。
3.4 BMJ推荐
①预测建模：目标是追求最优预测精度，因此允许所有连续变量灵活建模（如样条函数、多项式）通常是最佳选择，只要样本量足够支撑。若自由度有限，应优先给予对结果最重要的变量更大的建模灵活性，对次要变量保持线性近似通常影响不大。 ②变量重要性评估：在缺乏学科知识时，可用统计量（如partial偏卡方、秩相关）衡量预测重要性，用于决定哪些变量需要更多自由度。但这种评估不能揭示函数形式，也不应仅因关联度低就剔除变量。 ③函数形式选择的陷阱：使用散点图或单独统计检验来决定变量形式会引入“幽灵自由度（ghost degrees of freedom）”，导致偏差（如 P 值过小、置信区间过窄）。模型选择步骤应纳入整体建模和验证过程，而不是事后修正。 ④解释性建模：对主要自变量必须灵活建模，以捕捉真实关系；混杂因素也应避免强制线性或分类，否则会残余混杂。 ⑤描述性建模：重点是准确反映数据结构，因真实关系通常未知，应优先选择样条或分数多项式等灵活方法，而非假定固定分布。 ⑥样本量与自由度限制：小样本下模型的灵活性受限，若效应大或关系已知，可用简约模型；否则需在研究设计阶段通过合理的样本量计算来保障建模能力。
四、非线性方法
4.1 分数多项式fp
历史上，当怀疑存在非线性关系时，会使用例如二次（ x 2 ）或三次（ x 3 ）多项式来变换连续变量。然而，这些低阶多项式通常不能很好地拟合数据，而高阶多项式往往与连续变量分布的尾部拟合不佳。 分数多项式（Fractional Polynomials, FP），通过在有限候选幂次集内灵活选择变换（包括分数和对数），为连续变量提供比传统多项式更稳健、灵活且可控的一阶或二阶非线性建模方式，并能有效控制 I 型错误。常用一阶（β?x?1）或二阶（β?x?1 + β?x?2）即可拟合大多数情况。幂指数从集合 S={?2, ?1, ?0.5, 0, 0.5, 1, 2, 3} 中选择（x? 表示 log(x)）。候选指数有限，但可产生多种函数形式：一阶 8 种，二阶 32 种。 先前的研究表明，在许多情况下，分数多项式和样条函数会产生相似的结果。多项式和样条函数在外推方面的表现都不太好。rcs样条函数（强制尾部线性）则表现更好！
4.2 限制性立方样条rcs
使用样条函数，将连续变量划分为不同的区间；在每个区间内，用同次（例如三次）的多项式函数拟合数据。这些单独的函数被约束为在定义区间的节点处平滑连接。样条函数与分类之间的重要区别在于，生成的函数保持连续，并且自变量变化对因变量的影响不仅可能在区间之间变化，也可能在区间内变化。分段拟合到每个区间的函数取决于所使用的样条函数。最简单的样条函数是线性样条函数。然而，由于它们通常不能很好地拟合非线性关系，因此可以使用非线性多项式代替。最常用的样条函数是三次样条函数，并且存在不同的变体，例如平滑样条函数、惩罚样条函数、B 样条函数或限制三次样条函数。所有这些样条函数都已被证明能够充分模拟非线性关系，因此它们中的任何一种都是灵活建模的良好选择。 限制性立方样条函数也称为自然样条函数，是在尾部（即第一个结点之前和最后一个结点之后）对线性进行约束的三次分段多项式，这使它们更加稳健。无约束三次样条函数与分数多项式一样，在尾部表现不佳，此时可用的数据点很少，而异常值的影响很大。要使用受限三次样条函数，需要指定结点的数量和位置。结点的确切位置通常不会对结果产生重大影响，因此在大多数情况下，将它们放在变量边缘分布的固定分位数（例如，三个结点的第 10、第 50 和第 90 分位数）处是合适的 。选择结点的数量时应考虑样本大小，通常在 3 到 5 个之间。随着节点数量的增加，模型关系的灵活性也会提升，但代价是模型复杂度的增加（使用的自由度增加），以及过拟合的风险也会增加。然而，与其他三次样条函数变体或分数阶多项式相比，限制型三次样条函数消耗的自由度（节点数量-1）更少，因此通常可以减少过拟合。模型质量估计量，例如赤池信息准则，可用于确定最优节点数量，因为在模型拟合度相似的情况下，该指标更倾向于选择更简单的模型。
4.3 智荟rcscsi
我在公众号推广了无数次一键rcs且超美的rcscsi包，之前的推文合集多篇阅读量都过万，大家可以回溯下，lancet顶刊也在用。最新方法学发在med research，中文翻译https://mp.weixin.qq.com/s/a-s7h72Anvo14xoJJTevJg，大家可引用https://doi.org/10.1002/mdr2.70015 。 实战医学统计rcs合集一次看够，点击地址：https://mp.weixin.qq.com/mp/appmsgalbum?__biz=MzAxNTUzNDMyNQ==&action=getalbum&album_id=3015203637993930755&scene=173&subscene=&sessionid=svr_9ed8d0340cf&enterid=1756087521&from_msgid=2650953533&from_itemidx=1&count=3&nolastread=1#wechat_redirect
4.4 案例2：rcs寻找真相- y轴预测概率的部分条件效应图
与模型 3 相比，模型 1 未能捕捉到脑脊液葡萄糖与急性细菌性脑膜炎之间的非单调关系，因此预测不准确。模型 1 未能捕捉到脑脊液葡萄糖水平约为 75 毫克/分升（4.16 毫摩尔/升）时急性细菌性脑膜炎的最低风险，以及随后随着血糖值升高而增加的风险。由于缺乏灵活性，模型 1 也无法准确预测低脑脊液葡萄糖患者的概率大小，并且对最低脑脊液葡萄糖值给出了约 65%-70%的最大概率点估计值。模型 2 的表现可能更差，因为它将脑脊液葡萄糖水平分别为 10 毫克/分升和 60 毫克/分升的患者归因于相同的急性细菌性脑膜炎概率。根据模型 3，即脑脊液葡萄糖和急性细菌性脑膜炎之间真实关系的更接近值，这两名患者的概率会有很大差异 - 即 10 mg/dL 的概率为 100%，而 60 mg/dL 的概率约为 20%。如果不能充分估计连续独立变量和因变量之间的真实关联，也会导致模型的诊断鉴别力较差。与模型 3（0.96（0.94 至 0.98））相比，模型 1（0.87（95% CI 0.83 至 0.91））和模型 2（0.87（0.83 至 0.90））之间的模型鉴别力（使用受试者工作特征曲线下面积（AUC）评估）存在显著差异。
4.5 案例3：rcs寻找真相- y轴OR的剂量反应图
剂量反应图基于以模型中剩余变量为条件的预测，并且根据所进行的分析，不仅可以检查自变量和因变量之间的非线性关系，还可以确定例如最高和最低风险水平、药物的最佳剂量或治疗效果。效应指标（例如比值比）是根据一个参考点计算的，在该参考点处置信区间会崩溃。该参考点默认为所研究连续变量的中位数。虽然这种方法可以作为合理的默认方法，但参考点的选择应以先验的临床知识为指导。在缺乏此类信息的情况下，可以使用具有临床意义的阈值（例如生物标志物的参考下限或上限）作为替代方案。
五、其他非线性方法
存在其他方法，例如广义加性模型gam和局部估计的散点图平滑（loess）。机器学习领域还可以找到其他替代方法（例如，k-最近邻算法和人工神经网络）。
六、R代码
原文代码有些bug，我一并修复并添加了备注。重新整理了逻辑，方便与本文一同阅读。
1数据清理
# ===================== 准备环境（依赖、全局设置） =====================
# R version 4.4.0
if (!requireNamespace("pacman", quietly = TRUE)) install.packages("pacman")
pacman::p_load(
survival, # 生存/回归底层工具
rms, # 回归建模策略（lrm, rcs, Predict 等）
data.table, # 高效数据框
ggplot2, # 作图
patchwork, # 拼图
gtsummary, # 描述性表格
gt, # 表格美化
Hmisc # label()/units() 等工具（rms 依赖）
)
# 统一主题（黑白底，字号 12）
theme_set(ggplot2::theme_bw(base_size = 12))
# ===================== 工作目录（如不存在则跳过） =====================
wd <- "D:/聂个人文件/我的公众号/RCS/BMJ非线性指南/BMJ代码"
if (dir.exists(wd)) setwd(wd)
# ===================== 数据导入与清洗 =====================
# 说明：官方数据 abm.xls 容易格式不标准，已转为 abm.csv 更稳妥
# 若将来需要自动识别多格式，可再封装 smart_import。
df <- data.table::setDT(read.csv("abm.csv"))
# 只保留分析所需变量，并把 ".z" 视为缺失 NA
df <- df[, .(abm, age, sex, whites, gl)
][, lapply(.SD, function(x) replace(x, which(x == ".z"), NA))]
# 排除结局缺失（abm）
df <- df[!is.na(abm)]
# 变量重命名/类型转换：
# - whites -> csf_leuk（脑脊液白细胞计数）
# - gl -> csf_gl（脑脊液葡萄糖，mg/dL）
df[, csf_leuk := as.numeric(whites)][, whites := NULL]
df[, csf_gl := as.numeric(gl)][, gl := NULL]
df[, abm := as.numeric(abm)] # 0/1：病毒/细菌性脑膜炎（假设原意）
# —— 构造按中位数分组的分类型变量（仅用于分类建模与展示）——
median_csf_gl <- median(df$csf_gl, na.rm = TRUE)
median_age <- median(df$age, na.rm = TRUE)
median_csf_leuk <- median(df$csf_leuk, na.rm = TRUE)
df[, csf_gl_cat := factor(
fifelse(csf_gl < median_csf_gl, paste("<", median_csf_gl, "mg/dL"),
paste("≥", median_csf_gl, "mg/dL"))
)]
df[, age_cat := factor(
fifelse(age < median_age, paste("<", median_age, "years"),
paste("≥", median_age, "years"))
)]
df[, csf_leuk_cat := factor(
fifelse(csf_leuk < median_csf_leuk, paste("<", median_csf_leuk, "count/mm3"),
paste("≥", median_csf_leuk, "count/mm3"))
)]
# —— 变量标签与单位（用于 rms/Hmisc 展示）——
label(df$csf_gl) <- "CSF Glucose"
units(df$csf_gl) <- "mg/dL"
label(df$csf_gl_cat) <- "CSF Glucose"
2 模型123汇总
# ===================== datadist（rms 所需的分布信息） =====================
dd <- rms::datadist(df)
options(datadist = "dd")
names(df)
# ===================== 描述性统计表（按 abm 分组） =====================
# 用 gtsummary，直接在 label= 里改中文名；add_overall() 增加“总体”列
table1 <- gtsummary::tbl_summary(
data = df,
by = abm,
label = list(
## abm ~ "Acute bacterial meningitis",
age ~ "Age (years)",
sex ~ "Sex",
csf_leuk ~ "CSF Leukocytes",
csf_gl_cat ~ "Glucose (categorised)",
age_cat ~ "Age (categorised)",
csf_leuk_cat ~ "CSF Leukocytes (categorised)"
),
missing = "ifany",
missing_text = "Missing"
) |>
gtsummary::add_overall() |>
gtsummary::as_gt()
# ===================== 模型构建：rms::lrm 二项逻辑回归 =====================
# 模型1：连续型（线性假定）
mod1 <- rms::lrm(
abm ~ csf_gl + age + sex + csf_leuk,
data = df, x = TRUE, y = TRUE
)
# 模型2：分类变量（按中位数分层）
mod2 <- rms::lrm(
abm ~ csf_gl_cat + age_cat + sex + csf_leuk_cat,
data = df, x = TRUE, y = TRUE
)
# 使用 rcs()（受限三次样条）为 csf_gl/age/csf_leuk 构造非线性项
# 先用 AIC 选择 csf_gl 的结点数（knots）= 3/4/5
nk <- c(3, 4, 5)
aic_tbl <- data.table::rbindlist(lapply(nk, function(k) {
fit <- rms::lrm(
abm ~ rms::rcs(csf_gl, k) + rms::rcs(age, 3) + sex + rms::rcs(csf_leuk, 3),
data = df, x = TRUE, y = TRUE
)
data.table::data.table(knots = k, AIC = round(AIC(fit), 3))
}))
print(aic_tbl)
……..
# 汇总模型与系数表（用于比较）
ls_mod <- list("Model 1" = mod1, "Model 2" = mod2, "Model 3" = mod3)
ls_tab_coef <- lapply(ls_mod, function(m)
data.table::data.table(
variables = names(coef(m)),
coefficients = round(coef(m), 3)
)
)
3 画图
# ===================== 画图（线性预测器、概率曲线、剂量-反应） =====================
## 1) 线性预测器：单变量（csf_gl）三种表示
mod_uni <- rms::lrm(abm ~ csf_gl, data = df)
mod_uni_cat <- rms::lrm(abm ~ csf_gl_cat, data = df)
mod_uni_rcs <- rms::lrm(abm ~ rms::rcs(csf_gl, 3), data = df)
lp <- as.numeric(predict(mod_uni, type = "lp"))
lp_cat <- as.numeric(predict(mod_uni_cat, type = "lp"))
lp_rcs <- as.numeric(predict(mod_uni_rcs, type = "lp"))
data_part <- data.table::data.table(
csf_gl = df$csf_gl, lp = lp, lp_cat = lp_cat, lp_rcs = lp_rcs
)[, melt(.SD, id.vars = "csf_gl")]
……..
ggplot2::ggsave("plot_lp.tiff", plot = plot_lp, dpi = 300,
width = 7, height = 4, units = "in")
## 2) 概率曲线：分别基于 mod1/mod2/mod3
plot_pred1 <- ggplot2::ggplot(
rms::Predict(mod1, csf_gl, fun = plogis),
adj.subtitle = FALSE
) +
ggplot2::scale_x_continuous(breaks = seq(0, 150, 25)) +
ggplot2::labs(y = "Probability")
plot_pred2 <- ggplot2::ggplot(
rms::Predict(mod2, csf_gl_cat, fun = plogis),
adj.subtitle = FALSE
) +
ggplot2::labs(y = "Probability")
plot_pred3 <- ggplot2::ggplot(
rms::Predict(mod3, csf_gl, fun = plogis),
adj.subtitle = FALSE
) +
ggplot2::scale_x_continuous(breaks = seq(0, 150, 25)) +
ggplot2::labs(y = "Probability")
patch_prob <- (plot_pred1 | plot_pred2) / plot_pred3 +
patchwork::plot_annotation(tag_levels = "A")
ggplot2::ggsave("patch_prob.tiff", plot = patch_prob, dpi = 300,
width = 7, height = 4, units = "in")
## 3) 剂量-反应曲线（OR），默认参考值为中位数
pred_mod <- rms::Predict(mod3, csf_gl, fun = exp, ref.zero = TRUE)
plot_dr <- ggplot2::ggplot(pred_mod) +
ggplot2::coord_cartesian(xlim = c(20, 150), ylim = c(0, 20)) +
ggplot2::scale_y_continuous(name = "Odds ratio", breaks = seq(0, 20, 5)) +
ggplot2::scale_x_continuous(breaks = seq(20, 140, 20)) +
ggplot2::geom_hline(yintercept = 1,
color = scales::alpha("black", 0.8), linetype = 2) +
ggplot2::geom_vline(xintercept = stats::quantile(df$csf_gl, na.rm = TRUE)[3],
color = scales::alpha("red", 0.6), linetype = 2) +
ggplot2::labs(caption = NULL)
# —— 修改参考值（例如 50 mg/dL），需更新 datadist 并重拟合 ——
ref_value <- 50
dd$limits["Adjust to", "csf_gl"] <- ref_value # 调整参考
mod3_upd <- update(mod3) # 基于新参考值更新模型
……..
patch_dr <- plot_dr + plot_dr_upd +
patchwork::plot_annotation(tag_levels = "A")
ggplot2::ggsave("patch_dr.tiff", plot = patch_dr, dpi = 300,
width = 7, height = 4, units = "in")
发散思维
1、限制性立方样条用rcssci包可以一键判断aic、knot、8款形态，十分方便。官网cran直接安装rcssci包就可用，欢迎大家引用，https://doi.org/10.1002/mdr2.70015
2、BMJ推荐rcs应用于任意x，尤其是重要性x，可更小bias。
小结
“rcs非线性拟合终极利器，BMJ让你用起来。rcssci包，一键精美出图。
PS: 福利网盘，有效期1周，到20250901。包含BMJ原文、方法附件、BMJ源代码以及实战医学备注版复现代码。
需要的朋友，公众号后台回复（不要空格）： 0825
本公众号建立了学术交流群(群)，仅供SCI学术交流，人数有限需要实名制。入群请加笔者微信 popnie，请备注说明：姓名-学校(单位)-专业。