限制性立方样条图,一种美的不行的趋势性分析方法(附R语言详细教程)

2021年4月专题:趋势性分析方法系列1.最基本的趋势性检验方法2.利用回归进行趋势性分析3.限制性立方样条图,一种美的不行的趋势性分析方法4.Loess回归,非线性趋势性分析,怎么能少了它呢?

一般医学研究中,我们经常会通过构建回归模型来分析自变量和因变量之间的关系。但是,大多数回归模型有一个重要的假设条件:自变量和因变量呈线性关联。这个条件一般难以满足。当条件不满足时,可以①将连续型变量转化为分类变量,但是分类变量的类别数目以及节点位置的选择一般会带有主观性并且分类变量会损失部分信息;②也可以直接拟合自变量和因变量之间的非线性关系。

若自变量x与因变量y之间存在非线性关系时,常用的方法是绘制限制性立方样条图(Restricted cubic spline,RCS)。

这图太美了!

一、非线性关系与限制性立方样条图

非线性关系可以构建多项式回归或者样条回归来进行说明,但是直接构建多项式回归存在以下问题:①过度拟合②共线性③全局性(全局性是针对所有数据讲的,也就是说所有用来拟合的数据都需要符合多项式的规律),但当所有数据不能用一种关系来表示时(即可能有的数据在小于某个值之前是直线关系,在到了这个值以后是二次项关系),这种情况下需要构建样条回归,才可更准确的描述连续变量与结局之间的关系。

样条回归是由于数据在自变量取不同范围时有不同变化趋势,需要将数据分开,分别拟合模型,可以拟合直线、二次项、三次项式回归,其中拟合回归的类型根据实际情况而定。样条回归本质上其实是一个分段多项式,但它一般要求每个节点上连续且二阶可导,这样是为了保证曲线的平滑性。适用条件为①数据无法用一条直线描述②数据多项式回归效果一般好③本身想要了解某个事件前后变化的趋势④发现数据在某个节点前后趋势发生了改变。

其中如果数据全部看起来是一个变化趋势,那么没必要进行样条回归。样条回归和分段回归区别在于前者是加了约束条件的后者,因为有的变量是一个缓慢变化的过程,分段回归使每个段的内部效应被强制统一, 在节点的位置跳跃,“瞬时变化”不合理, 这不但不符合很多实际情况, 而且不能发现最大值和最小值的点。由于加了约束条件,正常情况下会导致其拟合效果稍差于分段回归,但是分析自变量和因变量之间的关系会更加合理。

样条回归往往会在曲线的两头,预测的区间会非常宽,因此需要再加一个边界限制条件,即限制性立方样条图。限制性立方样条是在回归样条的基础上再加一个约束条件,即样条函数在自变量数据范围两端的两个区间内为线性函数,这样使得两边的预测更为准确一些。

使用限制性立方样条图绘制非线性关系时,即将连续变量分为几段, 南京四小凤进行分段回归,通常需要设置样条函数截断值的个数以及位置。大多数情况下截断值的位置基本不会影响限制性立方样条的拟合效果,截断值的个数会决定曲线的形状。当节点个数等于样本量时, 相当于将各个点用线段相连, 得到的是完全拟合但是不平滑的折线。由于节点个数的选择和自由度有关, 所以当样本量比较大的时候可以取较多的节点。但是节点越多, 自由度越大, 模型越复杂, 越难解释。

在«Regression Modeling Strategies»这本书中,Harrell建议节点数为4时,模型的拟合效果较好,即同时可以兼顾曲线的平滑程度以及避免过拟合造成的精确度降低。当样本量较大时,5个节点是更好的选择。小样本(n<30)可以选择3个节点。当节点的个数为2时,得到的拟合曲线就是一条直线。大多数研究者推荐的节点为3-5个。

二、文献对RCS的描述

最近几年来经常见到在Jama、BMJ、Lancet 等杂志中见到别人利用限制性立方样条图对非线性关系的描述。当连续变量与结局为非线性关系时,使用RCS曲线可以更准确的描述连续型变量与结局发生的风险关系。

01 案例

图片

图片

图片

02 案例

图片

图片

图片

这些案例,大家都可以发生关键词“RCS”到公众号获取代码哦!

三、如何利用R语言绘制RCS图

目录(1)导入包和数据

(2)判断是否为非线性关系

(3)选择节点个数与位置

(4)作图

1、导入包和数据

library(rms) #RCS

library(survminer)#曲线

library(ggplot2) #作图

library(haven) #读取spss数据集

#清理运行环境

rm(list = ls())

#读入数据集(某病理类型胰腺癌数据集)

aa<-read_sav("pancer.sav")

#为后续程序设定数据环境(rms包函数很多都要做这两步)

dd<-datadist(aa)

#为后续程序设定数据环境

options(datadist='dd')

2、survminer包中的ggcoxfunctional()函数可直接绘图

#例如,我们检验数据集中age与风险的非线性关系

ggcoxfunctional(Surv(time, censor)~age+log(age)+sqrt(age),加盟服务data=aa)

这个研究可以进行age和生存的线性关系诊断,采用的方法是鞅残差,从下图age的1次、2次和3次方来看,可以看出年龄和死亡之间存在非线性关系。

图片

3、如何利用RCS曲线确认cut off值

节点个数:分段回归,使每段都符合线性假设。Harrel建议节点数,大多数研究者推荐的节点为3-5个。

#1.构建模型。【节点数为3.4.5均进行了构建和比较】

fit1<- cph(Surv(time,censor)~rcs(age,3),data=aa);anova(fit1);fit1

fit2<- cph(Surv(time,censor)~rcs(age,4),data=aa);anova(fit2);fit2

fit3<- cph(Surv(time,censor)~rcs(age,5),data=aa);anova(fit3);fit3

R²和Dxy越大,拟合的模型越优。

(1)3节点 R^2=0.193

图片

(2)4节点 R^2=0.245

图片

(3)5节点 R^2=0.248

图片

因此,节点取5,R²和Dxy值最大,拟合的模型更优。

2.计算风险比HR值与自变量(age)变化关系

HR<-Predict(fit3,age,fun=exp,ref.zero=TRUE);HR

图片

图片

3.可视化上述关系

P1<-ggplot(HR);P

图片

4.图形美化

P2<-ggplot()+ geom_line(data=HR, #数据来源           aes(age,yhat),#xy轴的数据           linetype="solid",#曲线加粗           size=1,           alpha=0.7,           colour="red")+ geom_ribbon(data=HR, #加入置信区间             aes(age,             ymin=lower,              ymax=upper),             alpha=0.1,             fill="red")+ theme_classic()+ geom_hline(yintercept=1,linetype=2,size=0.75)+ #y=1的水平线  labs(title="死亡风险随年龄的变化曲线:基于限制性立方样条(RCS)",       x="年龄 (即连续变量)",      y="HR (95%CI)"  );P2

图片

年龄£64岁,死亡风险随年龄变化不是很明显;年龄>64岁后,死亡风险随年龄的增加而增加。

四、R语言代码获取

过程有点复杂,有兴趣的朋友可以跟着实操一下,文中涉及的代码、英文文献,均可发送关键词“RCS”到公众号直接获取

本文就到此结束!

更多信息

本公众号作为医学数据分析公众号,提供一些免费医学统计学学习资源下载,欢迎点击下载。1.免费下载!统计初学者的福音!《妙趣横生统计学》视频,生动有趣的统计学!2.医学研究样本量如何计算?原创高清教程视频来了,完全免费下载!3.绝对值得收藏!原创高清SPSS 操作视频免费下载4.推荐!这个流行病大神制作的公共卫生研究小工具,可以计算标准化率及置信区间5.2006-2020中国卫生统计年鉴完整合集免费下载6.全网最简单的SPSS教程,160页PPT学会SPSS统计分析!免费下载!7.“如何在90分钟学会统计分析?” 来下载PPT学习吧!还有免费直播视频特别提醒:上述资源每天限分享和下载一个。

培训通告

2021年,我们召集了一批富有经验的高校专业队伍,着手举行短期统计课程培训班。如果您有需求,不妨点击查看:来参加吧,通俗易懂的统计培训课:R、SPSS、Meta、重复测量以及量表分析

图片

如果您觉得文章不错,为我们打“call”,点击“分享”吧