RCS限制性立方样条-高分套路全解析
本章思维导图：
近两年RCS和Nomogram一样炙手可热(泛滥成灾)，有必要做一个更专业更SCI化的专业梳理。出于此目的，汇总心得如下。
前言
临床连续数据与结局的线性关联方式与模型优化是一个重要问题。医学常见非线性曲线形态如下图，有S L log -log J threshold U ∩ 形态，数据拟合越接近真实关联曲线形态才越利于提高预测性能。统计学家George Box 的名言说到 “all models are wrong, but some are useful”。意思是，不能说某个构建的模型100%正确，只可以说它看起来似乎有点用。这几年热门的RCS(限制性立方样条)就是这样一种方法，它可以为模型拟合提供线索, 让模型看起来更好用，让数据分析更亮眼。
我们经常通过线性回归、Logistic回归、Cox回归探讨连续自变量和因变量的关系，但上述模型要求合适的链接函数g(y)，保证自变量和因变量服从一般/广义线性前提条件。其中线性前提条件转换可参考后续转换章节，例如正态转换、logitP转换。这里需要注意，上述假设不成立时，我们需要采用什么样的方法探讨自变量和因变量是线性还是非线性关联呢？平滑函数 非线性关系可以平滑工具来实现。平滑工具主要包括核平滑、多项式平滑和样条平滑。核平滑是使用一组局部权重来生成平滑的估计，在实际应用中很难实现；多项式是最简单的函数平滑结构，但只是一种粗略的估计方法，存在过拟合、共线性、全局性等诸多问题。样条(spline)本质上是一个分段多项式, 但它一般要求每个分段点上连续并且二阶可导，这样可以保证曲线的平滑性，如下图鸡蛋曲面knot可导。Spline适用条件为①数据x与y关系不符合直线或广义线性前提 ②数据多项式回归R2低③某个节点前后趋势发生明显改变。
RCS
RCS(Restricted Cubic Spline或Natural cubic splines )。立方样条本质上是连续平滑的分段三次多项式。其中“分段pieces”的数量由使用的“节点knots”数量决定。每个节点内实际上是一个三次多项式。经典的分段回归使每个段的内部效应被强制统一, 在节点的位置跳跃,“瞬时变化”不合理， 这不但不符合很多实际情况, 而且不能发现最大值和最小值的点。RCS的数学原理实际上比每个节点拟合三次多项式要复杂一些，需要施加进一步限制以便spline是连续平滑没有间隙。简而言之，RCS实质上是通过选择节点的位置和数量，拟合样条函数RCS（X），使得连续变量X在整个取值范围内呈现光滑的曲线，如图。Restricted RCS在回归样条的基础上附加Restricted：样条函数在自变量数据范围左右两端的两个区间内为线性函数。因为左右两端预测的区间非常宽加一个线性函数作为约束条件，这样使得两边的预测更为准确。由于加了约束条件，正常情况下会导致各段区间内RCS拟合效果稍差于分段回归，但是总体趋势上会更加合理。
RCS通过样条函数RCS（X）转换自变量X后，然后根据因变量的分布类型选择合适的链接函数，进而拟合模型g（Y）=常数项+RCS（X）+其他自变量，其中g为链接函数。样条函数RCS（X）包括一个线性项X以及K-2个立方项（S），即RCS（X）=β0X+β1S1+…+β（k-2）S（k-2），完整公式见下图。Ci(x)是落在第i个节点中的三次分量，g是所谓的链接函数。
其中，节点的位置对样条函数的拟合影响不大，一般根据连续变量的百分位数选择。节点数量对样条函数的影响较大，节点的数量决定曲线的形状，当节点的数量为2时，得到的拟合曲线就是一条直线。研究显示，节点数量为3~5时样条函数拟合较好，一般推荐knot=4。
表1 harrell 推荐knot 3-7
RCS形状受 knots影响很大，knot越多曲线约复杂。由此可见，RCS 只能用来判断非线性趋势，为切点选择提供线索，而不能精确提供切点！
SCI实战中，BMJ、EI文章对于RCS非线性关系描述的非常好，值得参考。例如，BMJ文章截取了其中的脂肪量和死亡风险的曲线图，作者采用中位数设置为拐点，从图中可以看出脂肪量在<21Kg，死亡风险基本保持不变，而>21Kg以后，死亡风险随着脂肪量增大而逐渐增高。通过切点前后的HR95%CI,侧面说明研究结论。通过绘制RCS图，既直观又能输出p for non-linearity值用来全面描述x与y之间的非线性关系，还可发现潜在的有趣切点，可以为文章增色不少。RCS对切点的线索，多是采用knot 3-8个且视觉上95%CI最狭窄处，但是很多文章都没有交代，似乎很多都默认了。R包会自动给出对应的切点位置。通过SCI分析，我们也发现knot对应的分位数也并不是一成不变的，如BMJ文献knot=4 (reference value= median，percentage= P5 P35 P65 P95)，EI文献knot=4 (reference value= P10，percentage= P20 P40 P60 P80)。主流RCS表达套路如下：
① RCS+ HR左右开弓，切点median
② RCS+ 直方图，切点自定义
③ RCS+ 密度图，切点自定义
发散思维1：选择分位数还是RCS
连续数据，我们往往第一选择会依据临床共识进行分类，常见的如二分类法（如高血压 SBP≥140mmhg。有时关注变量并没有公认的分类切点，此时我们会选择中位数或多分类/等级分类(三分位、四分位)。这些分类虽然严格遵循临床意义或分位数原则，但是可能削弱模型预测能力。忽略连续数据非线性形态，而直接引入分类进行模型预测，显然容易丢失许多信息，如下图，具体测试可在https://drjgauthier.shinyapps.io/spliny/ 中自行感受。简而言之，能RCS就先RCS, RCS后再根据图形定性形态，选择对应的knots点、临床界值、ROC界值、中位数、分位数、segmented包分段对应切点，进行分析。
发散思维2：其他非线性拟合
RCS很强，但是非线性拟合过程也有其他方法，https://www.r-bloggers.com/2014/09/an-exercise-in-non-linear-modeling/ 对主流的几类非线性方法进行了比较。为了避免过度拟合，根据AIC最小值标准选择模型。结论是RCS的AIC最小，拟合效果最佳。所以一句话，看不明白RCS准没错。B样条，别名基本样条，是RCS的常见替代品，后者也使用节点来控制灵活性。与RCS区别在于左右末端不受限制，因此它们比受限制的三次样条具有更灵活的尾部。广义加性模型 (GAM) 是广义线性模型的扩展，专门研究非线性关系，GAM用途很广。所用惩罚参数，薄板回归样条、三次回归样条、P 样条(惩罚B样条)。GAM不光用来为非线性关系提供线索，更重要的是直接参与GAM模型构建，尤其是环境研究领域常用，详细查看相关章节。
表 spline 平滑相关R包
R软件代码
SAS在unconditional logistic, pooled logistic, conditional logistic, 甚至是 proportional hazards regression models里均可实现RCS，但是绘图结果较丑陋，不推荐！这里如果感兴趣的同学深入研究，也会发现SAS 提供的 p for non-linearity 值定义是 P >0.05 为非线性关联，与R提供的P<0.05 为非线性关联相反，一定要注意这个坑。STATA也可以跑RCS，也很便捷，但是依旧存在图形需要转换再编辑问题，也就不推荐了。由于日程使用，临床数据大多结局为Y=01, 考虑OR或HR为效应指标。本程序仅提供R语言logistic 与 cox主流代码。
rm(list=ls()) #清除内存
pacman::p_load(rms, survival)
setwd("D:/合作者/黄素丽/高血压重金属")
library(haven)
#dt <- read_sav("D:/合作者/黄素丽/高血压重金属/metal-hyper.sav")
load(file="rt.Rdata")
y <- "hypertension"
x <- "LnAl"
covs <- c("age", "BMI", "gender", "smoking", "drinking", "UA","dm","hyperlipidemia","eGFR")
#indf 新增 xy
indf <- rt
indf[, "y"] <- indf[, y]
indf[, "x"] <- indf[, x]
# 删除自变量x covs中，有缺失值的95行， 统计RCS用数据集 indf
sum(!complete.cases(indf[, c(y, x, covs)]))
indf <- indf[complete.cases(indf[, c(y, x, covs)]),]
#设定数据环境
dd <- datadist(indf)
options(datadist="dd")
#========测试二分类 logistic为主分析，AIC最小 选择knot=3;lrm 测试logistic
## 获取 knot minAIC 方法1,=3
aics <- NULL
for (knot in 3:10) {
formula <- paste0(
"y ~ rcs(x, ", knot, ")", " + ", paste0(covs, collapse=" + ")
)
model <- lrm(as.formula(formula), data=indf, x=TRUE)
summary(model)
aics <- c(aics, AIC(model))
}
aics
#1521.189 1523.673 1521.889 1523.185 1524.816 1525.510 1524.746 1525.997；knot=3 最小
seq(3, 10)[which.min(aics)] #3
# 方法2, 需要输入covs,nk=3
for (knot in 3:10) {
fit <- glm(y~rcs(x,knot)+age+BMI+gender+smoking+drinking+UA+dm+hyperlipidemia+eGFR,data=indf,family=binomial() )
tmp <- extractAIC(fit)
if(knot==3){AIC=tmp[2];nk=3}
if(tmp[2]<AIC){AIC=tmp[2];nk=knot}
}
nk #3
#获取非线性检验p， x=LnAl,knot=3
#y为二分类，lrm函数拟合：fit<- lrm(y ~ rcs(x1),data=data)
#y为连续变量，ols函数拟合模型：fit<-ols(y~ rcs(x1)
#y为 time相关随访结局，cph函数拟合
#采用默认 knot=3，P10 P50 P90为位置参数
fit <- lrm(y~rcs(LnAl,3)+age+BMI+gender+smoking+drinking+UA+dm+hyperlipidemia+eGFR,data=indf,x=T,y=T)
anova(fit)
## 获得x的5%分位间隔的分位数，以此获得P5，P35，P50，P65，P95等等..分位
quantile(indf$LnAl,probs = seq(0,1,0.05))
#采用自定义 knot=3，P20 P50 P80为位置参数
fit.LnAl <- lrm(y~rcs(LnAl,c(2.74,3.17,3.71))+age+BMI+gender+smoking+drinking+UA+dm+hyperlipidemia+eGFR,data=indf,x=T,y=T)
anova(fit.LnAl)
# p-non-linear 0.0107 或 0.0079 p<0.05, 非线性趋势
p <-round(anova(fit.LnAl)[,3],3)
Pre_OR.LnAl <-rms::Predict(fit.LnAl,LnAl,fun=exp,type="predictions",ref.zero=T,conf.int = 0.95,digits=2)
ggplot(Pre_OR.LnAl)
# y-hat =3.23 是OR=1切点
View(Pre_OR.LnAl)
#对数据进行打包，并指定参考值默认median value=3.17; 本例value=3.23，OR=1，
#refvalue <- quantile(indf[, "x"], prob=0.5)
refvalue <- 3.23
ddist <- datadist(indf)
##这里是设置参考点，也就是HR为1的点，常见的为中位数或者临床有意义的点
ddist$limits$LnAl[2] <- refvalue
options(datadist="ddist")
# OR计算，fun是转化函数
Pre_OR.LnAl <-rms::Predict(fit.LnAl,LnAl,fun=exp,type="predictions",ref.zero=T,conf.int = 0.95,digits=2)
ggplot(Pre_OR.LnAl)
#========测试time数据 cox为主分析，AIC最小 选择knot=3; cph 测试cox
pacman::p_load(haven,smoothHR,survival,rms,Hmisc)
#测试cox
#indfcox <- read_sav("D:/合作者/黄素丽/高血压重金属/rcscox.sav")
load(file="indf.Rdata")
y <- "hypertension"
x <- "LnAl"
covs <- c("age", "BMI", "gender", "smoking", "drinking", "UA","dm","hyperlipidemia","eGFR")
#indf 新增 xy
indf[, "y"] <- indf[, y]
indf[, "x"] <- indf[, x]
# 删除自变量x covs中，有缺失值的95行， 统计RCS用数据集 indf
sum(!complete.cases(indf[, c(y, x, covs)]))
indf <- indf[complete.cases(indf[, c(y, x, covs)]),]
#test 数据集 indf需带time status
quantile(indf$LnAl,probs = seq(0,1,0.05))
# 获得x的5%分位间隔的分位数，以此获得5，35，50，65，95等等分位
#对数据进行打包，寻找y-hat=1
ddist <- datadist(indf)
options(datadist="ddist")
S <- Surv(indf$time,indf$status==1)
#获取多因素 cox分析 AIC 最小 knot, 单变量LnAl
for (knot in 3:10) {
fit <- cph(S~rcs(LnAl,knot)+gender,data=indf,x= TRUE, y= TRUE, surv = TRUE)
tmp <- extractAIC(fit)
if(knot==3){AIC=tmp[2];nk=3}
if(tmp[2]<AIC){AIC=tmp[2];nk=knot}
}
nk #4
#---------基于cox绘制 单因素LnAl 的 RCS
pacman::p_load(rms,survminer,ggplot2,ggsci)
#构建cph 函数获取rcs, 单指标LnAl；也可校正其他因素 S ~ rcs(LnAl,4)+ BMI 等
fit.LnAl <- cph(S ~ rcs(LnAl,4)+gender, x=TRUE, y=TRUE,data=indf)
# PH 检验,P>0.05 符合PH假设
cox.zph(fit.LnAl, "rank")
#残差图的横轴是时间，纵轴是残差，残差均匀分布则表示残差与时间相互独立，满足ph假设
ggcoxzph(cox.zph(fit.LnAl, "rank"))
#非线性检验p-non-linear，0.9289
anova(fit.LnAl)
p <-round(anova(fit.LnAl)[,3],3)
## HR计算，fun是转化函数
Pre_HR.LnAl <-rms::Predict(fit.LnAl,LnAl,fun=exp,type="predictions",ref.zero=T,conf.int = 0.95,digits=2)
ggplot(Pre_HR.LnAl)
# y-hat=1.00000, 对应LnAl= 3.16
View(Pre_HR.LnAl)
# 根据y-hat =1 调整ref，指定参考值对数据再次进行打包，OR=1
#refvalue <- 3.16
#ddist <- datadist(indf)
#ddist$limits$LnAl[2] <- refvalue
#options(datadist="ddist")
#S <- Surv(indf$time,indf$status==1)
# 全人群的COX模型 RCS绘图
rcs_all <- ggplot()+
geom_line(data=Pre_HR.LnAl,
aes(LnAl,yhat),# x y轴距
linetype="solid",#曲线加粗
size=1,
alpha=0.7,
colour="purple")+
scale_color_nejm()+ ##采用ggsci包中英格兰调色，也可以其他
geom_ribbon(data=Pre_HR.LnAl,
aes(LnAl, ymin=lower,ymax=upper,fill="purple"),alpha=0.1)+
theme_classic()+
scale_fill_nejm()+
geom_hline(yintercept=1,linetype=2,size=0.75) #y=1水平线
labs(title ="风险随LnAl变化曲RCS",
x="LnAl",
y="HR (95%CI)"
)
rcs_all
# 分性别的，分层人群的COX模型 RCS绘图；gender=1 男性，2女性
# 男性无拐点，女性拐点3.16
Pre1 <- rms::Predict(fit.LnAl, LnAl, gender=c('1','2'),fun=exp,type="predictions",ref.zero=TRUE,conf.int = 0.95,digits=2)
par(mfrow=c(1,2))
ggplot(Pre1)
View(Pre1)
###---直方图+平滑曲线----
## HR计算，fun是转化函数；承接上述fit.LnAl
Pre0 <-rms::Predict(fit.LnAl,LnAl,fun=exp,type="predictions",ref.zero=T,conf.int = 0.95,digits=2)
par(mar=c(3,4,1,5),new=T)
#mar是设置图形边界空白宽度，四个数值分别代表，下左上右的边界宽度，取值单位是线条宽度lwd 。new=T表示可以添加新的图层。
#####-----------------------------
col=c('#5BA5DA','#F06955','#FDC059','skyblue','blue') ## 蓝 红 黄
hist(indf$LnAl,axes=F,xlab="",ylab="",
xlim=c(1.5,4.5),
col=col[4],main='',freq=T)
axis(4) ##加右边Y轴,1，2，3，4分别代表下左上右
par(new=T)## 创建新的画布，类似于ggplot2中的图层叠加
plot(Pre0[,1],axes=T,Pre0$yhat,type='l',lty=1,lwd=2,
#Plot中type绘制散点图的类型常用的6种类型，其中“p”表示绘制散点；“b”表示点连线；
#“l”表示绘制线图；“s”表示绘制阶梯；“o”表示绘制点连线，此处线条穿过每一个点，与“b”类型稍有不同
#“n”表示不绘制任何图形元素。其中s就是绘制生存曲线KM的阶梯线图
col=col[2],xlim = c(1.5,4.5),
ylim=c(0,3),xlab='LnAl',ylab='HR(95%CI)')
lines(Pre0[,1],Pre0$lower,type='l',lty=2,lwd=2,col=col[2])
lines(Pre0[,1],Pre0$upper,type='l',lty=2,lwd=2,col=col[2])
mtext('freq',side=4,line = 3,outer = F)
##图形添加散点
points(x=ddist$limits$LnAl[2], ##散点X轴坐标
y=1, ####散点y轴坐标
pch=19,##pch是点的形状，数字从1-25代表不同的点形状。
col='red')
Text <- paste('Ref=',round(ddist$limits$LnAl[2],3),sep='')
##图像添加文字，这里添加参考点
text(Text,x=ddist$limits$LnAl[2]+5,y=0.5,cex = 1,font = 3)
#abline(h=1,col='black')
outfile <- "fig_rcs.pdf"
pdf(outfile, onefile=FALSE, width=15/2.54, height=15/2.54)
plot
#####直方图+平滑曲线+置信区间填充
col=c('#E0E7F2','#EE947B')
par(mar=c(4,4,1,5))
#####-----------------------------
hist(indf$LnAl,axes=F,xlab="",ylab="",
#ylim=c(0,0.035),
xlim=c(1.5,4.5),
#两个plot的X轴范围要保持一致，保持对齐
col=col[1],main='',freq=T,border = 'gray')
axis(4) ##加右边Y轴
par(new=T)## 创建新的画布，类似于ggplot2中的图层叠加
plot(Pre0[,1],axes=T,Pre0$yhat,type='l',
lty=1,lwd=2,
col=col[2],
xlim = c(1.5,4.5),
ylim=c(0,5),xlab='LnAl',ylab='HR(95%CI)')
c= rgb(1,0,0,alpha=0.15)
polygon(c(Pre0[,1],rev(Pre0[,1])),
c(Pre0$lower,rev(Pre0$upper)),col = c,border = c)
mtext('freq',side=4,line = 3,outer = F)
##图形添加散点
points(x=ddist$limits$LnAl[2], ##散点X轴坐标
y=1, ####散点y轴坐标
pch=19,##pch是点的形状，数字从1-25代表不同的点形状。
col='red')
abline(h=1,col='black',lty=2)
outfile <- "fig_rcs.pdf"
pdf(outfile, onefile=FALSE, width=15/2.54, height=15/2.54)
plot
dev.off()
outfile <- "fig_rcs.tiff"
tiff(outfile, width=15, height=15, unit="cm", res=600, compression="lzw+p")
plot
dev.off()
1. https://www.nature.com/articles/s41409-019-0679-x.
2. 王晓晓, 陶立元, 李楠, 等. 限制性立方样条在非线性关联分析中的应用 [J] . 中华儿科杂志,2020,58 (08): 652-652. DOI: 10.3760/cma.j.cn112140-20200615-00621
3. Lee DH, Keum N, Hu FB, et al. Predicted lean body mass, fat mass, and all cause and cause specific mortality in men: prospective US cohort study. BMJ. 2018;362:k2575. Published 2018 Jul 3. doi:10.1136/bmj.k2575
4. Wen Y, Huang S, Zhang Y, et al. Associations of multiple plasma metals with the risk of ischemic stroke: A case-control study. Environ Int. 2019;125:125-134. doi:10.1016/j.envint.2018.12.037
5. R语言绘图|3.基于Cox+限制性立方样条：直方图+平滑曲线+双坐标轴绘制， 易学统计。
原创不易，转载请说明来自本公众号。