Stata实现RCS限制立方样条：探索非线性剂量-反应关系

三年前我拿到一份随访数据，BMI和死亡风险，Cox模型里只放BMI线性项，P=0.31，完全不显著。但把BMI按五分位分组画K-M曲线，趋势清晰得吓人——最低组和最高组的死亡风险都比中间组高出一截。这个矛盾让我意识到，连续性暴露变量的效应，往往就藏在"模型的线性假设"里。后来改用RCS（限制立方样条，Restricted Cubic Spline）重新建模，U型关系一目了然，非线性P值0.008。这篇攻略就是把从那次之后踩过的坑、积累的Stata实操经验完整写出来，适合要做剂量-反应关系分析、准备SCI论文投稿或者复核审稿人意见的医学统计从业者。

1. 为什么RCS是探索连续变量关联形态的首选

1.1 线性模型看不清的"隐藏形态"

线性回归或者Cox回归里放一个连续变量X，本质上是假设X每增加一个单位，结局风险就固定变化一个常数倍。这个假设在真实医学数据里经常不成立。BMI和死亡风险是教科书级的反例：过瘦和过胖死亡风险都升高，中间是谷底，关系呈U型。如果只用线性项拟合，得到的HR其实是被"平均化"的斜率——它既低估了高BMI的风险，也掩盖了低BMI的风险，甚至可能因为两端同时升高而相互抵消，得出"不显著"的结论。

有人会想：那我用四分位数分组，把BMI变成分类变量不就行了？分组确实能看出趋势，但是代价很大。第一，组别划分方式会影响结果，五分位和四分位的结论可能不一样；第二，组内信息全部丢失，一组几万人被压缩成一个哑变量；第三，分组比较得到的是"组间差异是否显著"，而不是"剂量-反应关系形态如何"，审稿人一句"分组切点是否有生物学依据"就能问到哑口无言。

1.2 多项式拟合的三大硬伤

既然线性不行，那直接放二次项BMI^2、三次项BMI^3行不行？这是很多人的第一反应，我也走过这条路。多项式回归在x取值范围内拟合效果尚可，但有三个硬伤。

第一是全局性。多项式每个系数都作用于整个定义域，局部某个区间的弯曲会牵动整条曲线变形，尤其当数据在某个区间比较稀疏时，曲线会被强行拽过去。第二是边界震荡。这就是数值分析里著名的Runge现象——次数越高，边界处振荡越剧烈，BMI取到35以上时曲线可能突然飙上去再跌下来，毫无医学解释。第三是外推灾难。三次以上的多项式在数据范围外会急速发散，你拿到一个BMI=52的新样本，模型预测的风险可能是负的。RCS正是针对这三个痛点设计的：局部拟合、边界线性、限制外推。

RCS的思路很直接：用几个节点把X的取值范围切成若干段，每一段分别用三次多项式拟合，同时保证连接点处曲线光滑连续，这就是"样条（Spline）"；所谓"限制立方"，是指首尾两段强制变成直线，从而避免边界处不可控的抖动。用大白话说，RCS就是"分段三次多项式，节点处平滑衔接，两头拉直"。

2. RCS的数学直觉，不需要懂公式也能理解

2.1 分段三次多项式与节点的作用

假设我们选了4个节点，记为k1、k2、k3、k4，把BMI的取值范围切成3段。每一段里用一个三次多项式去拟合数据，同时在k2、k3这两个内部节点处，要求左右两边的函数值、一阶导数、二阶导数都相等。这个"相等"的约束保证了线段之间没有折角，曲线看起来是浑然一体的光滑弯曲。

节点数量决定了模型的灵活度。节点越多，曲线能捕捉的细节越多，但也越容易跟着噪音走。Harrell在《Regression Modeling Strategies》里给出过一个广为接受的经验：4到5个节点（也就是df=3到4）对绝大多数医学数据都够用了。节点位置一般放在预测变量的分位数上，比如5%、35%、65%、95%百分位，而不是均匀等距。因为样条在数据稀疏的地方本来就不稳定，放在分位数上能确保每段都有足够样本支撑。

2.2 "限制"在哪里：首尾段线性

"限制立方样条"里的"限制"，学术上叫linear tail restriction，即第一个节点之前和最后一个节点之后的曲线被强制为线性。这个设计极其重要。如果不加限制，样条在两端外推时，三次项会让曲线像脱缰野马一样乱跑，样本量一少，95%置信区间就张成一个大喇叭口。限制之后，两段变成直线，外推行为可控得多。

有人会问：那既然两头是直线，为什么不用纯粹的分段线性（linear spline）？分段线性也就是折线，每条线段是直线，节点处是折角，曲线不光滑，且局部拟合能力弱，对真实曲线的逼近效果远不如三次样条。所以RCS是"光滑"和"稳定"之间最好的折中。

2.3 rcsgen生成的变量与自由度df

Stata和R中，RCS的实现逻辑同源。Stata里最常用的命令是rcsgen，它也是直接借鉴了R的rms包里的rcs()函数。运行rcsgen x, df(4)会生成4个衍生变量，对应4个自由度，也就意味着在模型里要放4个系数。

这里有一个理解上极为关键的点：生成的这些基函数变量里，第一个分量承载的是X的线性趋势，剩下几个分量专门负责"偏离线性"的部分。所以后面做非线性检验时，只需要联合检验第2到第4个变量的系数是否同时为0，如果它们显著不等于0，就说明数据中存在线性模型无法解释的弯曲。理解了这一点，P值解读就有了理论基础。

3. Stata全流程实操：生成变量、拟合模型、出P值

3.1 准备数据与安装rcsgen

先用Stata自带的NHANES II随访数据做演示，这份数据包含人群的BMI、年龄、性别和死亡结局，非常适合复现RCS分析。

webuse nhanes2f, clear stset t2death, failure(death) des bmi death age female

安装rcsgen：

ssc install rcsgen

3.2 生成RCS基函数并拟合Cox模型

rcsgen bmi, df(4) gen(bmi_rcs) list bmi bmi_rcs1 bmi_rcs2 bmi_rcs3 bmi_rcs4 in 1/5, sep(0)

默认情况下，df(4)会在BMI的5%、35%、65%、95%分位数放置4个节点，生成4个基函数变量bmi_rcs1~bmi_rcs4。你也可以手动指定位置，用knots(22 27 32 38)这样的方式，但我建议在多数场景下直接用默认分位数节点，避免人为干预带来的选择性偏差。

拟合Cox模型：

stcox bmi_rcs* age female

输出的LR chi2会显著高于只放bmi的线性模型。最关键的是看下面两条检验命令。

3.3 核心检验：整体关联性与非线性检验

* 整体关联性检验：4个基函数系数是否同时为0 testparm bmi_rcs* * 非线性检验：剔除线性分量后，其余基函数是否同时为0 testparm bmi_rcs2 bmi_rcs3 bmi_rcs4

第一条命令的P值回答的问题是：BMI和死亡风险到底有没有关系？这个P值相当于全局Wald检验，对"任何形式的关联"都敏感。第二条命令的P值回答的问题是：这种关系能不能用一条直线描述？如果它显著，说明曲线的弯曲是真实的统计学信号，而不是随机波动。

如果你做的是二分类结局，把stcox换成logistic或logit即可，检验命令完全一样。线性回归结局则用regress。

4. 结果解读：两个P值怎么看、怎么写进论文

4.1 整体关联P值回答"有没有关系"

整体关联P值检验的是四个基函数系数同时为0的原假设。如果P<0.05，说明在调整协变量之后，BMI与死亡风险存在统计学关联，至于是直线还是曲线、是正向还是负向，这个检验不负责回答。它的意义类似于回归模型整体的F检验，是我们对外报告"这个变量有效"的底气。

实际操作里有一个细节：整体关联检验的是"在节点位置既定条件下"的关联。节点位置变了，基函数矩阵就变了，检验结果也会微调。所以规范做法是在方法部分写明节点的数量和位置。

4.2 非线性P值回答"是不是直线"

非线性P值检验的是bmi_rcs2、bmi_rcs3、bmi_rcs4这三个系数的联合显著性。如果这三个系数都是0，那曲线就退化成一条直线，说明BMI每增加一个单位风险对数值固定变化，此时没必要用RCS这种复杂模型，直接报告线性HR即可。如果它们不全为0，曲线就是弯的。

注意一个常见误区：有人把非线性P值当作"RCS模型的整体P值"来报告，或者在非线性P值不显著时依然声称发现了U型关系，这会被审稿人一眼看穿。两个P值各司其职，不能混用。

4.3 四种组合判断表与论文报告句式

下面这个判断表是我在实际项目里反复使用的一个参考框架：

论文中可以直接使用的句式，以BMI和死亡风险为例：

"采用限制立方样条Cox回归模型探索BMI与全因死亡风险的剂量-反应关系，节点置于BMI分布的5%、35%、65%、95%分位数。结果显示BMI与死亡风险显著相关（全局Wald检验P=0.007），非线性检验提示二者呈非线性关系（非线性P=0.012）。"

这样写，审稿人要的信息——模型类型、节点位置、两个检验结果——都齐了，基本不会再追问方法学细节。

5. 画图与汇报：别让审稿人挑出刺

5.1 用marginsplot快速看趋势

模型跑完后，第一件该做的事是画图确认曲线形态，光看系数是看不出U型还是倒U型的。最快的方式是配合margins和marginsplot：

stcox bmi_rcs* age female margins, at(bmi=(15(1)45)) atmeans marginsplot

这段代码的含义是：把BMI从15到45每隔1取一个点，其他协变量固定在样本均值，计算每个点的预测相对风险并绘制曲线。得到的图会清晰展示风险随BMI变化的走势，是否有谷底、是否有平台期，一目了然。如果你想展示的是预测概率而非相对风险，可以先跑一个logistic模型再用相同方式画，趋势形态一般差别不大。

5.2 标准剂量反应曲线与参考值设定

如果你要投稿，审稿人更希望看到的是以某个参考值为基准、HR=1的剂量反应曲线。在Stata里可以直接用社区命令postgrsp绘制：

ssc install postgrsp stcox bmi_rcs* age female postgrsp bmi, test1(25)

test1(25)表示以BMI=25作为参考值，该点的HR=1，曲线在这个点必然经过1。实际使用中我会建议把参考值设在中位数或者临床公认的正常切点，比如BMI取25，而不是设在曲线的某个极端位置，否则整条曲线的形状会被强行改变，误导读者。

5.3 节点、CI、极端值：图表中的三个细节

第一，图中要标注节点位置，至少在图注里写清楚"节点置于5%、35%、65%、95%分位数"，这是可复现性的基本要求。第二，置信区间是必须带的，没有CI的剂量反应曲线和折线图没有区别，审稿人看到光秃秃一条线通常会直接打回。第三，注意BMI极端值的处理。nhanes2f里BMI有超过50的个例，样条尾部会出现一个大喇叭口，因为极端值样本量极少，95%CI会急剧变宽。此时常见做法是把作图范围限定在2.5%~97.5%分位数，或者做缩尾处理，并在方法部分说明。

6. 实战中容易翻车的五个细节

6.1 节点数不是越多越好

节点数过多是新手最容易犯的错误。df设到6甚至8，曲线会开始捕捉个体噪音，表现为波浪形抖动，局部出现毫无医学意义的"驼峰"，而非线性P值因为自由度膨胀变得极小，看起来无比显著，实则全是伪信号。我的经验是：单变量RCS先试df=3和df=4，对比AIC或BIC，如果两者差异不大，取df=3更保守稳妥。

6.2 样本量与自由度怎么权衡

小样本下样条模型极易过拟合。粗略经验是每个节点段至少要有20~30个事件，如果你做的是Cox模型，要看的是"死亡例数"而不是总样本量。亚组分析时尤其要警惕：总样本5000但某个亚组只有200例，还硬上df=4，结果必然不稳定。此时建议df=3，甚至2都可以考虑——自由度的减少换来的模型稳定性是值得的。

6.3 亚组分析与交互项的处理

很多人拿到整体结果后喜欢分亚组分别跑RCS，比如分男女各跑一遍。这种做法本身没错，但要注意两点。第一，各组样本量差异会导致自由度上的可比性变差，最好在全样本模型里放交互项来正式检验效应修饰作用。一个可行的做法是生成RCS基函数后，分别生成与分组变量的交互项，再做似然比检验。第二，亚组RCS的节点位置尽量沿用主分析的位置，不要每组重新优化，否则组间比较就会混入"节点选择差异"这一额外变量。

6.4 稳健标准误对检验的影响

如果你在模型里加了vce(robust)或者vce(cluster id)，testparm会自动基于稳健方差矩阵做Wald检验，数学上没有问题。但要留意：cluster数量太少时Wald检验会偏激进，此时更稳妥的办法是用似然比检验。Stata里可以先估计无约束模型和约束模型，再用lrtest，不过样条的系数约束不是单一系数等于某个值，手动操作比较麻烦，所以多数时候我们仍然依赖testparm，但要在心里清楚它的Wald性质。

6.5 R与Stata结果的对应关系

如果你在R里用的是rms包的rcs()函数，Stata的rcsgen与其同源，默认节点位置也是一致的，因此两个软件算出的P值、曲线形状通常高度接近。但画图的便捷度上R更强，rms的Predict()配合ggplot可以轻松输出带CI且以参考值为基线的标准RCS图。有些团队的习惯是Stata算数、R出图，这种跨软件工作流也是可行的，只要在方法部分统一写清楚节点位置和参考值即可。

我在实际使用中还有一个习惯：所有RCS分析结果一定回查原始数据，尤其是曲线拐点附近，看看有没有个别极端样本把曲线"拽"出奇怪的形状。任何统计方法都不能替代对数据的仔细审视，RCS作为一个灵活的非线性工具，如果用它的人对底层数据没概念，灵活性反过来就会变成危险的自由度。