限制性立方样条图 R语言实现

限制性立方样条
多次回归样条往往会在曲线的两头，预测的区间会非常宽，因此需要再加一个边界限制条件，即限制性立方样条图。

使用RCS曲线可以更准确的描述连续型变量与结局发生的风险关系。 

节点的选择
在利用限制性立方样条绘制曲线关系时，通常需要设置样条函数节点的个数和位置。绝大多数情况下, 节点的位置对限制性立方样条的拟合影响不大, 而节点的个数则决定曲线的形状, 或者说平滑程度。当节点的个数为2时, 得到的拟合曲线就是一条直线。

节点的个数选择
大多数研究者推荐的节点为3-5个。在《Regression Modeling Strategies》这本书中，Harrell建议节点数为4时，模型的拟合较好，同时可以兼顾曲线的平滑程度和避免过拟合造成的精度降低。而当样本量较大时，例如因变量为未删失的连续变量并且大于100时，5个节点是更好的选择。小样本（如n<30）可以选择3个节点。以下是Harrell推荐的节点数和相应的节点位置，大家可以参考。 

节点的位置选择
平均分配

选取剧烈变化的区域

数据变化复杂的地方多设置节点，更稳定的地方少设置节点

交叉验证 

实战及代码复现
1.导入包和数据
n <- 1000

set.seed(731)

age <- 50 + 12*rnorm(n)

label(age) <- "Age"

sex <- factor(sample(c('Male','Female'), n, rep=TRUE, prob=c(.6, .4)))

cens <- 15*runif(n)

h <- .02*exp(.04*(age-50)+.8*(sex=='Female'))

time<- -log(runif(n))/h

label(time) <- 'Follow-up Time'

death<- ifelse(time <= cens,1,0)

time <- pmin(time, cens)

units(time) <- "Year"

data<-data.frame(age,sex,time,death) #以上代码都是随机生成一份数据，无需过分关注，只需将自己的数据导入即可

#加载所需要的R包

library(rms) #回归/RCS的包

library(ggplot2) #作图

2.做回归、拟合曲线
View(data) #查看数据信息

dd <- datadist(data) #为后续程序设定数据环境

options(datadist='dd') #为后续程序设定数据环境

fit<- cph(Surv(time,death) ~ rcs(age,4) + sex,data=data) #限制性样条回归

HR<-Predict(fit, age,fun=exp,ref.zero = TRUE) #拟合曲线

3.作图
P1<-ggplot(HR)#####用ggplot2直接画图

P1

P2<-ggplot()+geom_line(data=HR, aes(age,yhat),linetype="solid",size=1,alpha = 0.7,colour="red")+
 
 geom_ribbon(data=HR, aes(age,ymin = lower, ymax = upper),alpha = 0.1,fill="red") #优化曲线

P2

P2<-P2+theme_classic()+geom_hline(yintercept=1, linetype=2,size=1)+
 
 labs(title = "RCS", x="age", y="HR (95%CI)") #优化曲线

P2

4.线性或逻辑回归参考
#线性回归示例参考
ggplot(data=data, aes(x=age, y=time)) +  
 geom_point()+
 stat_smooth(method = lm, formula = y ~ rcs(x,4))  

#逻辑回归示例参考
fit3 <- lrm(death ~ rcs(age,4)+sex,data=data)
OR <- Predict(fit3, age,fun=exp,ref.zero = TRUE)
ggplot(OR) 

原文链接：限制性立方样条（Restricted cubic spline，RCS） - 哔哩哔哩