Stata实现RCS限制立方样条：节点选择、P值解读与绘图实战

前两天有个师弟拿着一张审稿意见来找我，说审稿人让他把BMI和代谢风险的关系用RCS限制立方样条画出来，他对着Stata界面半天不知道从哪下手。这种需求这几年越来越常见——RCS几乎成了剂量反应关系研究的标配工具，但Stata里跑起来并没有R语言那么顺手，软件自带帮助里的信息也偏零散。这篇文章我就从实战出发，把Stata做RCS的完整流程、节点怎么选、P值怎么读、图上那些坑怎么躲，一次性说透。

内容主要面向三类人：一是临床研究者，需要在论文里报告非线性关联却一直卡在操作上；二是流行病学方向的研究生，正在处理队列数据或剂量反应Meta分析；三是想用RCS做敏感性分析、应对审稿人质疑的同行。看完你至少能独立跑出一条规范的多变量RCS曲线，并且知道论文里该怎么写、P值该怎么解释。

1. RCS限制立方样条到底是什么——先说清楚“为什么用”

1.1 从线性假设的尴尬说起

回归模型的基本假设里，暴露X对结局Y的影响是一条直线，每增加一个单位的X，Y的改变量恒定。这个假设在绝大多数生物学指标面前站不住脚。拿BMI和全因死亡率的关系来说，临床上早已观察到典型的J型曲线：过瘦和过胖的人死亡风险都升高，某个中间范围风险最低。你把这种关系硬塞进线性模型，两端的升高会迫使拟合线整体倾斜，等你看到结果时，中段平坦区域的风险被严重高估，甚至可能得出“BMI越低越好”这种反直觉结论。

传统分段建模的思路是把X切成几段，每段单独拟合一条直线或低阶多项式。表面看解决了曲线拟合问题，实际操作中却让人头疼：切点位置是人为主观定的，切点不一样结果就跟着变，审稿人通常会追问“为什么选这个切点”；更麻烦的是分段模型在各个切点处往往不连续，曲线呈现出明显的折角甚至跳跃，违背了生物学上“剂量反应关系通常是渐变”的基本认识。

RCS的提出就是为了同时解决这两个问题，它用若干个三次多项式分段描述暴露与结局的关系，但在相邻曲线的连接点处强制满足函数值、一阶导数和二阶导数都连续，整条曲线看起来就像一条一次成型的光滑曲线，没有折角，也不需要你人工指定拐点。

1.2 RCS的核心优势在于“限制”二字

RCS里“限制”指的是在两个最外侧节点之外，曲线被强制约束为线性。为什么要加这个约束？因为纯三次多项式在数据范围之外极具“破坏力”，样本边界外的预测值会不受控制地上下翻滚，一条本来很平滑的曲线尾部会像鱼尾一样甩出去。加上两端线性约束后，曲线在边界外延续为一条稳定的直线，尾部的抖动被大幅抑制，这在样本边缘数据稀疏时尤为重要。

从实现角度看，RCS通常只需要三到五个基函数就能捕捉常见的J型、U型、倒U型、阈值效应等形态，自由度开销远低于loess这类非参数平滑方法。更关键的一点是，RCS本质上仍然是带参数的线性模型，系数可以直接用常规的Wald检验和似然比检验做推断。你在论文里常看到的“非线性P值”，就是通过对RCS中负责弯曲的那部分基函数做联合检验得到的。这一条让RCS在统计推断和可视化之间取得了很好的平衡。

实操提示：Stata的mkspline命令生成的是三次样条基变量，在实际建模使用中已经被广泛接受为RCS的Stata实现方案。如果审稿人要求严格意义上的限制立方样条，需关注两端线性约束的细节，但这并不影响常规剂量反应关系研究的结论。

2. 建模前的准备：数据、命令与三个关键决策

2.1 数据要求与常规预处理

RCS对数据类型的要求并不特殊。结局变量可以是连续值（用线性回归）、二分类变量（用logistic回归），也可以是生存时间数据（用Cox回归），暴露变量必须为连续型变量。第一步不是急着跑模型，而是先看暴露变量的分布。特别注意数据分布在两端是否过于极端：比如某个变量90%的人集中在5到10之间，剩下10%的人散布在10到100，这种数据会让最外侧节点的位置完全由少数极端个体决定，曲线尾部基本没有说服力。

常规做法是用centile命令查看暴露变量的第5、35、65、95百分位数，这几个位置通常是后续设置节点的依据，同时确认每个节点附近都有足够多的观测。缺失值处理同样不能马虎，如果直接使用complete case分析，要确认剔除比例不大，并且缺失机制没有明显偏倚；比例较高时建议先做多重插补再建模。

实操提示：我习惯在建模前把暴露变量的分布画出来，直方图和箱线图都看一眼。如果一个候选节点附近样本量只有几十个人，我会主动把节点往中间挪，这是RCS建模里非常容易被忽略的一步。

2.2 节点数量与位置的选择逻辑

RCS的结果受两个因素影响最大：节点数量和节点位置。节点数量方面，文献中最常见的是3到5个。节点越少曲线越平滑、越接近线性；节点越多拟合越灵活，但越容易把噪声当信号，造成过拟合，自由度增加还会削弱检验功效。Harrell的建议是大多数情况下4个节点足够捕捉常见的非线性形态；如果样本量特别大，比如上万例的队列研究，可以考虑5个节点去检验更复杂的形状。

节点位置一般按暴露变量的分位数确定，最常用的是第5、35、65、95百分位数对应4个节点。这样做的原因是保证各节点附近样本量相对均衡，避免节点落在数据稀疏区导致数值不稳定。也有研究者使用第10、50、90百分位数，目前没有统一标准。我的建议是不要纠结于“哪个分位数组合更标准”，而是做一次敏感性分析：分别用3节点、4节点、5节点跑一遍，看曲线形状和关键P值方向是否一致。如果一致，节点选择问题不大；如果不一致，说明数据可能撑不起一个稳定的非线性模型，这时候应优先考虑减少节点数而不是试图找到“最好的那一组”。

2.3 mkspline与外部命令的取舍

Stata里实现RCS主要有两条路线。第一条是用内置的mkspline命令生成样条基变量，然后放进任意回归命令。这样做的优势是兼容性最强、完全可控，基变量名由自己定义，后续的test、predict、esttab等操作都能顺畅衔接。第二条是使用SSC上的外部命令，比如postrcspline或rcspline，它们可以在模型拟合后自动画出RCS曲线并输出非线性P值，适合快速出图。但这些外部命令的封装程度较高，遇到复杂模型或多变量调整时反而不如手搓灵活。

我自己更倾向正式建模和写论文时走mkspline路线，把每个环节掌握在自己手里；画图阶段再结合外部命令来优化图形效果。两条路线并不冲突，关键在于至少先把第一条走通，后续扩展起来遇到问题才知道根源在哪。

3. 实战全流程：从拟合到可视化（含Stata代码）

3.1 最小可行示例：连续变量RCS建模

直接用Stata自带的auto数据集做一个连续结局的最小示例。假设我们要研究weight对price的影响，先按百分位数确定节点，再生成样条基变量。

* 1. 查看weight的百分位数分布，确定节点位置 sysuse auto, clear centile weight, centile(5 35 65 95) * 2. 生成RCS样条基变量：4个节点对应3个基变量 mkspline w1 w2 w3 = weight, cubic knots(1800 2550 3200 4200) * 3. 拟合线性回归 regress price w1 w2 w3 * 4. 整体关联检验：三个基变量做联合检验 test w1 w2 w3 * 5. 非线性检验：只检查偏离线性的部分 test w2 w3

这里有个关键点要说清楚。mkspline生成的三个基变量中，w1主要捕捉的是weight与price之间的线性趋势部分，w2和w3捕捉的是偏离线性的弯曲部分。因此：

• 对w1、w2、w3做联合检验，得到的是“weight与price是否存在关联”的整体P值；
• 只对w2和w3做联合检验，得到的是“关联形状是否显著偏离线性”的非线性P值。

如果你的数据不是auto，而是自己的研究数据，节点位置的数字不要照抄，必须根据第一步centile输出的实际百分位数来填。mkspline对knots列表的要求是严格从小到大排列，如果顺序不对会直接报错。

3.2 多变量调整与协变量处理

真实研究中极少只有暴露和结局两个变量，混了一切都要调整。多变量调整并不改变RCS的核心逻辑，只需要把协变量作为额外回归项放进模型即可：

* 二分类结局示例 logistic death w1 w2 w3 age sex * 生存结局示例 stset follow_time, failure(death) stcox w1 w2 w3 age sex

多变量场景下有一个容易踩的坑：样条基变量生成后，必须在同一个样本里完成后续的回归和检验。有些读者在清理缺失值后重新生成样条变量，导致基变量对应的样本变了，节点位置也随之改变，模型结果前后对不上。我建议在建模开始前先针对完整样本处理好缺失值，再一次性生成样条基变量，后面就不再动样本结构了。

关于亚组分析，还有一个常见误区。做性别亚组时，有些研究者会分别对男性和女性各建一个模型，然后看到男性P显著、女性P不显著，就得出结论说“效应存在性别差异”。这在统计上是不成立的，亚组间差异是否显著需要通过交互项检验，也就是在全样本模型中纳入width与sex的交互项，然后对交互项做联合检验，而不是比较两个独立模型的P值大小。

3.3 绘制RCS曲线图的细节调优

绘制RCS曲线，本质上是把样条基变量对应的预测值随暴露变量的变化趋势画出来。最直接的画法是先预测每个观测的风险值，再按暴露变量排序连线：

* 拟合后预测，并按weight排序画线 predict pprice sort weight twoway line pprice weight, sort

这在不含协变量的最简模型中是可用的。一旦模型里加了age、sex等协变量，每个观测的预测值里包含各自的协变量贡献，这时候画出来的曲线就不是“纯暴露效应”了。要展示纯效应，需要把协变量固定在某个水平，然后用外部命令来辅助绘图更省事：

ssc install postrcspline logistic death w1 w2 w3 age sex postrcspline weight, nk(4)

postrcspline会自动生成RCS曲线并在图上标注非线性P值，非常方便。不过不同版本的命令封装方式有差异，装好后建议先跑一下help postrcspline确认参数。如果这个命令在你当前版本不可用，手工绘图的替代方案是建立预测网格、固定协变量后自行计算预测值，工作量会大一些。

正式投稿的图还需要做一些美化：在曲线旁边或下方添加直方图展示暴露变量的分布，让读者直观看到数据支撑范围；在参考值位置画一条垂直参考线；如果Y轴是风险或HR，可以考虑把坐标轴标签设置为自然对数刻度，置信区间展示会更对称。最后导出时选择矢量格式如.eps或.pdf，方便期刊排版。

4. P值解读：非线性检验到底在测什么

4.1 整体关联P值与非线性P值的区别

RCS分析中不可能绕开两组P值，很多误解也出在它们身上。我把它们的区别拆成一张表：

整体P值的零假设是“所有样条基变量的系数都为0”，即暴露与结局完全没有关联。非线性P值的零假设是“负责弯曲的基变量系数都为0”，即数据可以用一条直线来描述而不损失太多拟合优度。

很多人把这两个P值用反了。整体P值不显著，不能断言暴露与结局“没关系”，只能说在当前样本量下没有足够证据拒绝“无关联”的零假设，真实情况可能是效应量太小或者样本量不足。非线性P值不显著，也不等于关系就是线性的，只是在线性模型下拟合尚可，除非你有明确理由强调弯曲，否则按线性趋势报告是合理的。

更具体地说，组合情况有四种。整体P显著、非线性P不显著时，报告为“存在显著关联，未观察到显著的非线性偏离，趋势呈线性”是最佳表达。两者都显著时，可以报告非线性关联存在，并描述曲线形态，比如J型、倒U型或阈值效应。整体P不显著、非线性P显著的情况比较罕见，通常解释为数据整体不支持关联，但某个局部区间可能有风险变化，需要结合临床意义谨慎讨论。两者都不显著就老老实实报告未观察到显著关联。

4.2 三个P值如何对应到论文报告

论文里常出现的趋势P值、整体关联P值、非线性P值三者各有定位。趋势P值在很多文献里直接用暴露变量进模型后的系数P值，有时也用序数化后的P值；整体关联P值和非线性P值在RCS语境下更规范，部分期刊还要求报告自由度，例如“P for overall association = 0.002 (df = 3)”，这样做的好处是读者能直接看出检验的复杂程度。

我自己的报告模板大致是：

• 如果非线性P值小于0.05：报告RCS曲线，附整体P值和非线性P值，明确描述曲线形态，可以标注参考值对应的HR或OR；
• 如果非线性P值大于等于0.05但整体P值小于0.05：报告为线性趋势，附趋势P值，不再强调曲线弯曲；
• 如果两者都不显著：报告未观察到显著关联，除非研究本身是验证性的假说检验，否则不适合过度解读。

需要注意一点，很多期刊只要求报告“非线性P值”，审稿人默认你用的是RCS或者类似平滑方法。如果审稿人追问自由度，就如实报告检验使用的基变量数量，这一般等于节点数减1。

4.3 常见错误解读与大忌

最典型的错误是把非线性P值当成“整体关联是否存在”的指标。有人看到非线性P值等于0.08，就写成“变量与结局之间没有显著的非线性关系，但近似线性且显著”，这句话本身没问题；但换成“变量与结局无关联”就完全错了。一个变量完全可以是显著线性的，非线性P值同样不显著，两者并不矛盾。

另一个高频错误是对曲线尾部的小波动过度解读。RCS曲线两端的置信区间通常很宽，因为极端分位数附近的样本量稀少。如果尾部出现一个看似陡峭的下降或上升，需要先看该位置的置信区间是否跨越了无效应线。如果跨越，大概率是噪声，不能作为拐点证据，更不能因此得出“低于某个值开始风险升高”这种精确结论。

还有一个大忌：直接从图上“看”出拐点位置，然后把暴露变量在拐点处二分或三分，再做分组比较。这种做法在统计上极为脆弱，因为拐点位置本身带有很大的不确定性，完全忽略这种不确定性会导致多重比较和过度拟合问题。审稿人对此非常敏感。更合理的做法是保留RCS的连续结果来描述趋势，或者使用两段线性样条做显式的阈值检验。

避坑提醒：如果审稿人要求你报告“拐点及95%置信区间”，这就意味着你应该使用分段回归或阈值分析的方法，而不是从RCS图上目测一个位置。两者的统计推断逻辑不同，务必分清楚。

5. 常见问题与排查技巧实录

5.1 节点数不同结果差异巨大怎么办

如果你跑完3节点、4节点、5节点三种设置，曲线形态基本一致、P值方向不变，那结论很稳健，可以在论文里报告一种设置，同时在补充材料中展示敏感性分析结果。如果节点一变结论就翻转，优先检查两件事：一是样本量，尤其要怀疑非线性项是否被少数极端个体驱动；二是节点位置是否落在数据稀疏区。

排查方法很直接：在数据集中找出各节点附近的观测数。如果某个节点两侧加起来不到几十个观测，换节点位置或者减少节点数量通常能解决问题。还有一个细节，不同回归类型对节点个数的敏感程度不同，Cox模型往往比线性回归更容易出现数值波动，生存数据里事件数不足时尤其明显。

5.2 曲线端头“发飘”的应对

RCS虽然加了端部线性约束，但两端仍然可能因为数据稀疏出现宽置信区间，甚至整条曲线在尾部大幅度摆动。我遇到这种情况时一般做三个检查：第一，看数据中是否存在极端离群值，考虑是否为记录错误，比如前文提到的BMI等于80这种明显不合理数值；第二，尝试把最外侧节点向中间移动，比如从5%和95%百分位数改成10%和90%百分位数，这样会牺牲一部分曲线形状，但能显著提高尾部稳定性；第三，如果前两步都无效，可能真的需要更多样本，这只能在研究设计阶段解决，分析阶段能做的有限。

5.3 参考值（Reference）如何设定

论文里的HR或OR必须相对于某个参考值才有意义。RCS模型的输出默认以模型截距对应的预测值为基线，但这个值不一定有临床意义。比如研究BMI与死亡风险，临床上普遍习惯以BMI等于25作为参照，就需要在建模后显式设定参考值。

Stata中实现这一点的基本思路是：在预测网格上生成对应样条基变量的值，计算出参考值处和每个网格点处的线性预测值，两者相减后取指数，就得到以参考值为基准的OR或HR曲线。具体可以用lincom逐步计算，也可以在预测数据集中用generate手动构造。这一步没有捷径，耐心处理好会让论文结果的可解释性提升一个档次。

5.4 结果导出与报告模板

模型结果建议用esttab导出到Word或Excel，论文正文不需要列出每个基变量的系数，只需要报告非线性P值和曲线图。但如果期刊要求补充材料，把关键模型的结果表放进去也很常见。esttab的基本用法：

eststo model1: logistic death w1 w2 w3 age sex esttab model1 using rcs_results.rtf, replace

如果只需要导出P值，可以用putexcel手动构建结果表。曲线图则务必输出矢量格式，避免用截图直接放入论文。论文里的标准描述句大概是：“采用限制立方样条拟合BMI与全因死亡风险的剂量反应关系，节点设置于第5、35、65、95百分位数，以BMI=25为参考，结果显示两者呈J型关联（整体P<0.001，非线性P=0.003）。”这样一句就足以概括核心结果。

最后说句掏心窝的话。我见过太多人把RCS当成一个出图工具，跑通就完事，完全不理解背后的假设和检验逻辑。等审稿人问一句“你的非线性P值为什么用Wald检验，自由度是多少”，就卡壳了。我的建议是：正式投稿前，把节点数从3换到5跑一遍，把参考值换一个再跑一遍，把协变量调整方案换一套再跑一遍。如果结论始终稳定，你再把图放进论文，心里是有底的；如果不稳定，正好趁早发现，总比审稿时被打回来强。