回归样条本质上是一个分段多项式，但它一般要求每个分段点上连续并且二阶可导。在建立限制性立方样条的时候，节点（knots）数量的选择决定曲线的形状, 或者说平滑程度，大多数研究者推荐的节点为3-5个。在《Regression Modeling Strategies》这本书中，Harrell建议节点数为4时，模型的拟合较好，同时可以兼顾曲线的平滑程度和避免过拟合造成的精度降低。
但是我们作为医学统计学的运用者，就和之前郑老师说的那样，我们只要学会把一头猪放进机器等着香肠制作完成就好，至于步骤是怎么制作的，交给统计学家去完成就好了。
适用条件
首先我们要知道为什么去使用限制性立方样条。
1、将连续型变量转化为分类变量，但是分类变量的类别数目以及节点位置的选择一般会带有主观性并且分类变量会损失部分信息；
2、可以直接拟合自变量和因变量之间的非线性关系，于是我们常用的方法就是绘制限制性立方样条图。
具体条件如下：
①数据无法用一条直线描述；
②数据多项式回归效果一般好；
③本身想要了解某个事件前后变化的趋势；
④发现数据在某个节点前后趋势发生了改变。
下面我将结合案例给大家运行一下如何应用R语言进行RCS的绘制
本例子以睡眠时间为例，当你想进一步了解具体睡眠时间对心血管之间的变化趋势，但又不能用线性关系去解释的时候，这就可以用到限制性立方样条。
首先安装加载一些包
#包的安装
install.package(“foreign”)
install.package(“ggplot2”)
install.package(“rms”)
install.package(“survival”)
install.package(“Hmisc”)
install.package(“splines”)
加载相关的包
#包的加载，安装成功的包才能加载
library(“foreign”)
library(“ggplot2”)
library(“rms”)
library(“survival”)
library(“Hmisc”)
library(“splines”)
然后导入数据集
操作所需数据请大家在公众号后台回复“样图”获取
#设置工作空间，需要在D盘下建一个R文件夹，后把cc数据放进去
setwd("D:/R")
ae <- read.csv("cc.csv")
names(ae)
attach(ae)
变量因子化
ae$sex<-as.factor(ae$sex)
ae$bmi1<-as.factor(ae$bmi1)
接着为后续程序设定数据环境，也就是打包数据，这一步在预测模型中也常做
dd <- datadist(ae)
options(datadist='dd')
然后我们拟合logistic回归的限制性立方样条
【心血管与高血压之间的关系（xinxueguan是结局指标） （da049是睡眠时间，3是拟合曲线的时候采用三个节点）后面的+都是一些协变量，当然我没有完全放进去】
#拟合回归方程
fit<-lrm(xinxueguan~rcs(da049,3)+sex+bmi1+age+wenhua+xiyan,data=ae)
an<-anova(fit)
以上的3这个节点数量是自己选的，严谨起见还是得让系统选择一个节点数
#通过以下循环语句选择拟合最好的结点
for (i in 3:7) {
fit <- lrm(xinxueguan~rcs(da049,3)+sex+bmi1+age+wenhua+xiyan,data=ae)
tmp <- extractAIC(fit)
if(i == 3) {AIC = tmp[2]; nk = 3}
if(tmp[2] < AIC) {AIC = tmp[2]; nk = i}
}
*AIC实际上是评估拟合效果的，选AIC绝对值最小的就行，涉及到复杂的统计原理，我也不是很清楚
生成预测值，并作图 （exp相当于把预测值转换成了OR）
plot(Predict(fit, da049,fun=exp), anova=an, pval=T)
OR<-Predict(fit, da049,fun=exp,ref.zero = TRUE)
剩下来就是画图
这么简单语句就可以了
ggplot(OR)
#后面的一串代码就是修饰
ggplot()+geom_line(data=OR, aes(da049,yhat),linetype=1,size=1,alpha = 0.9,colour="red")+
geom_ribbon(data=HR, aes(da049,ymin = lower, ymax = upper),alpha = 0.3,fill="red")+
geom_hline(yintercept=1, linetype=2,size=1)+theme_classic()+
labs(title = "RCS", x="Sleep duration", y="OR (95%CI)")
从图里可以看到睡眠时间和冠心病发生的风险呈现显著的非线性关系，睡眠时间小于6小时会显著增加心血管疾病的发病率，并且在睡眠时间高于9小时，心血管疾病的发病率也有增加的趋势，整体呈现U型关系。
当然我们也可以对性别进行分层，代码如下
#计算不同性别的OR值
OR1 <- Predict(fit, da049, sex=c('1','2'),
fun=exp,type="predictions",
ref.zero=TRUE,conf.int = 0.95,digits=2)
剩下就是画图，以下一大串的就是进行图行绘制，基本上不用进行更改
ggplot()+
geom_line(data=OR1, aes(da049,yhat, color = sex),
linetype="solid",size=1,alpha = 0.9)+
geom_ribbon(data=OR1,
aes(da049,ymin = lower, ymax = upper,fill = sex),
alpha = 0.2)+
scale_color_manual(values = c('red','blue'))+
scale_fill_manual(values = c("red","blue"))+
theme_classic()+
geom_hline(yintercept=1, linetype=2,size=1)+
labs(title = "CVDs Risk", x="Sleep duration", y="OR(95%CI)")
HR1 <- Predict(fit, da049, sex=c('1','2'),
fun=exp,type="predictions",
ref.zero=TRUE,conf.int = 0.95,digits=2)
结果图
这样就可以看出男性1和女性2之间差别了，结果显示同样呈现U型趋势的情况下，睡眠时间对女性心血管疾病的发病更为易感。
另外，如果是生存分析的话，拟合回归样条的曲线的代码也很相似，区别在于一个因变量是二分类的，一个是生存资料，所以在拟合模型时，会有一定差异
fit<- cph(Surv(livetime,xinxueguan) ~ rcs(da049,3)+sex+age+ bmi1+ wenhua+xiyan,data=be)
其他的代码都是类似
详情请点击下方：
https://mp.weixin.qq.com/s/bELmEO513R0VbBoe55OYcw
返回搜狐，查看更多