基于rms包的限制性立方样条回归（RCS）R代码实现-Toy模板网

这篇具有很好参考价值的文章主要介绍了基于rms包的限制性立方样条回归（RCS）R代码实现。希望对大家有所帮助。如果存在错误或未考虑完全的地方，请大家不吝赐教，您也可以点击"举报违法"按钮提交疑问。

1.原理

限制性立方样条(Restricted cubic spline，RCS)是分析非线性关系的最常见的方法之一。RCS用三次函数拟合不同节点之间的曲线并使其平滑连接，从而达到拟合整个曲线并检验其线性的过程。可以想见，RCS的节点数对拟合结果来说非常重要。通常，小于30个样本数的小样本取3个节点，大样本取5个节点。

2.R实现

1.cox回归

#Used for RCS(Restricted Cubic Spline)
#我们使用rms包

library(ggplot2)
library(rms)
library(survminer)
library(survival)

在这里我们使用survival包中的lung数据

#####基于cox回归
#这里用survival包里的lung数据集来做范例分析
head(lung)

# inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
#    3  306      2  74   1       1       90       100     1175      NA
#    3  455      2  68   1       0       90        90     1225      15
#    3 1010      1  56   1       0       90        90       NA      15
#    5  210      2  57   1       1       90        60     1150      11
#    1  883      2  60   1       0      100        90       NA       0
#   12 1022      1  74   1       1       50        80      513       0
#status:	censoring status 1=censored, 2=dead
#sex:	Male=1 Female=2
#ph.ecog：医生对患者的体能状态评级
#ph.karno：医生对患者的另一种体能状态评级karnofsky
#pat.karno:患者karnofsky自评
#meal.cal:摄入卡路里
#wt.loss:过去半年体重减轻



# 对数据进行打包，整理
dd <- datadist(lung) #为后续程序设定数据环境
options(datadist='dd') #为后续程序设定数据环境


#用AIC法计算不同节点数选择下的模型拟合度来决定最佳节点数
for (knot in 3:10) {
  fit <- cph(Surv(time,status) ~ rcs(meal.cal,knot) + sex+age , x=TRUE, y=TRUE,data=lung)
  tmp <- extractAIC(fit)
  if(knot==3){AIC=tmp[2];nk=3}
  if(tmp[2]<AIC){AIC=tmp[2];nk=knot}
}
nk  #3


#cox回归中自变量对HR的rcs
fit <- cph(Surv(time,status) ~ rcs(meal.cal,3) + sex+age , x=TRUE, y=TRUE,data=lung)#大样本5节点，小样本(<30)3节点
#比例风险PH假设检验,p>0.05满足假设检验
cox.zph(fit,"rank")
#非线性检验，p<0.05为有非线性关系
anova(fit)
#这里的结果是
#               Wald Statistics          Response: Surv(time, status) 
#
# Factor     Chi-Square d.f. P     
# meal.cal    0.42      2    0.8113
#  Nonlinear  0.09      1    0.7643，呈线性
# sex         6.61      1    0.0101
# age         1.99      1    0.1582
# TOTAL      10.29      4    0.0358
#查看各meal.cal对应的HR值
HR<-Predict(fit, meal.cal,fun=exp)
head(HR)
#画图
ggplot()+
  geom_line(data=HR, aes(meal.cal,yhat),
            linetype="solid",size=1,alpha = 0.7,colour="#0070b9")+
  geom_ribbon(data=HR, 
              aes(meal.cal,ymin = lower, ymax = upper),
              alpha = 0.1,fill="#0070b9")+
  theme_classic()+
  geom_hline(yintercept=1, linetype=2,size=1)+
  labs(title = "Lung Cancer Risk", x="Age", y="HR (95%CI)")

画出来的结果如图：

r 做rcs,r语言,回归

anova(fit)的结果以及可视化呈现显示在RCS下饮食能量摄入和死亡率之间的关系呈线性。

我们再找一个非线性的例子，这个例子使用survival包的colon数据：

#结肠癌病人数据
Colon <- colon
Colon$sex <- as.factor(Colon$sex)#1 for male,0 for female
Colon$etype <- as.factor(Colon$etype-1)
dd <- datadist(Colon)
options(datadist='dd') 
for (knot in 3:10) {
  fit <- cph(Surv(time,etype==1) ~ rcs(age,knot) +sex , x=TRUE, y=TRUE,data=Colon)
  tmp <- extractAIC(fit)
  if(knot==3){AIC=tmp[2];nk1=3}
  if(tmp[2]<AIC){AIC=tmp[2];nk1=knot}
}
nk1 #3
fit <- cph(Surv(time,etype==1) ~ rcs(age,3) +sex , x=TRUE, y=TRUE,data=Colon)
cox.zph(fit,"rank")
anova(fit)
#               Wald Statistics          Response: Surv(time, etype == 1) 
#
# Factor     Chi-Square d.f. P     
# age        6.90       2    0.0317
#  Nonlinear 6.32       1    0.0120，非线性
# sex        1.20       1    0.2732
# TOTAL      7.54       3    0.0565
HR<-Predict(fit, age,fun=exp)
head(HR)
ggplot()+
  geom_line(data=HR, aes(age,yhat),
            linetype="solid",size=1,alpha = 0.7,colour="#0070b9")+
  geom_ribbon(data=HR, 
              aes(age,ymin = lower, ymax = upper),
              alpha = 0.1,fill="#0070b9")+
  theme_classic()+
  geom_hline(yintercept=1, linetype=2,size=1)+
  geom_vline(xintercept=47.35176,size=1,color = '#d40e8c')+#查表HR=1对应的age
  geom_vline(xintercept=65.26131,size=1,color = '#d40e8c')+
  labs(title = "Colon Cancer Risk", x="Age", y="HR (95%CI)")

结果如下： r 做rcs,r语言,回归

anova(fit)的结果以及可视化呈现显示该示例呈现非线性关系

我们还可以进行分组的研究和可视化呈现，只需修改Predict()函数中的参数

HR1 <- Predict(fit, age, sex=c('0','1'),
               fun=exp,type="predictions",
               conf.int = 0.95,digits =2)
HR1
ggplot()+
  geom_line(data=HR1, aes(age,yhat, color = sex),
            linetype="solid",size=1,alpha = 0.7)+
  geom_ribbon(data=HR1, 
              aes(age,ymin = lower, ymax = upper,fill = sex),
              alpha = 0.1)+
  scale_color_manual(values = c('#0070b9','#d40e8c'))+
  scale_fill_manual(values = c('#0070b9','#d40e8c'))+
  theme_classic()+
  geom_hline(yintercept=1, linetype=2,size=1)+
  labs(title = "Colon Cancer Risk", x="Age", y="HR (95%CI)")

结果：

r 做rcs,r语言,回归

2.logistic回归

如果不是cox回归，用的是logistic回归，总体差不多，只需要替代模型就行

#####基于logistic回归的rcs
#建模型
fit <-lrm(status ~ rcs(age, 3)+sex,data=lung)  
OR <- Predict(fit, age,fun=exp)
#画图
ggplot()+
  geom_line(data=OR, aes(age,yhat),
            linetype="solid",size=1,alpha = 0.7,colour="#0070b9")+
  geom_ribbon(data=OR, 
              aes(age,ymin = lower, ymax = upper),
              alpha = 0.1,fill="#0070b9")+
  theme_classic()+
  geom_hline(yintercept=1, linetype=2,size=1)+
  geom_vline(xintercept=38.93970,size=1,color = '#d40e8c')+ #查表OR=1对应的age
  labs(title = "Lung Cancer Risk", x="Age", y="OR (95%CI)")

3.线性回归

线性回归也一样文章来源地址https://www.toymoban.com/news/detail-773379.html

#####基于线性回归的rcs
fit <- ols(meal.cal ~rcs(age,3)+sex,data=lung)
Kcal <- Predict(fit,age)
#画图
ggplot()+
  geom_line(data=Kcal, aes(age,yhat),
            linetype="solid",size=1,alpha = 0.7,colour="#0070b9")+
  geom_ribbon(data=Kcal, 
              aes(age,ymin = lower, ymax = upper),
              alpha = 0.1,fill="#0070b9")+
  theme_classic()+
    labs(title = "RCS", x="Age", y="Kcal")

到了这里，关于基于rms包的限制性立方样条回归（RCS）R代码实现的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！