一、引言
在医学和统计学领域，生存分析是一种分析个体生命长度和生存时间的重要方法。了解人们生存的期限和影响因素，对于制定健康政策、优化医疗资源的分配以及个体护理方案的制定都至关重要。传统的生存分析方法如Kaplan-Meier曲线和Cox比例风险模型已经广泛应用，但它们对数据的假设有一定的限制。
为了克服传统方法的限制，限制立方样条（Restricted Cubic Splines，简称RCS）成为一种强大而灵活的工具，可以更准确地估计个体生存函数。RCS是一种非参数的灵活拟合方法，通过将生存时间转化为各个节点上的分段函数来建模生存曲线。它可以适应各种生存时间分布，而不需要过多的假设。
本文旨在探讨RCS在生存分析中的应用和优势。首先，我们将介绍RCS的基本原理和建模方法。然后，我们将讨论RCS在生存分析中的应用案例，并与传统方法进行比较。接下来，我们将探讨RCS的局限性和可能的改进方向。最后，我们将总结本文的主要内容。
二、文献分享
-
「研究对象」：360万名成年人。
-
「研究目的」：探究BMI对死亡率的影响，并了解不同BMI范围内的死亡原因分布。研究采用了大规模队列研究的方法，跟踪参与者的BMI和死亡情况，通过对数据进行回归分析和风险评估，来评估BMI与总体和特定原因死亡率之间的关联。
-
「研究方法」：在这项基于人口的队列研究中，我们使用了与国家死亡登记数据相连的英国一级护理数据，数据来自临床实践研究数据联机(CPRD)。我们使用调整后的Cox回归模型，检查了BMI与全因死亡率的关联，以及BMI与广泛范围的特定死因（记录为国际疾病分类第10版[ICD-10]代码）之间的关联。我们包括所有在16岁及以上收集BMI数据并进行后续随访的个体。随访开始于以下情况中的最后一个：CPRD研究标准随访开始时间，第一个BMI记录的5周年纪念日，或1998年1月1日（死亡登记数据的起始日期）；随访终止于死亡或2016年3月8日。完全调整的模型按性别分层，并根据基线年龄、吸烟、酒精使用、糖尿病、多重剥夺指数和日历期调整。模型既适用于从未吸烟者，也适用于整个研究人群。我们还进行了广泛的敏感性分析。使用Poisson模型，包括BMI、年龄和性别，估计了40岁时基线的男性和女性在不同BMI范畴下的预期寿命。
-
「研究结果」：BMI与更具体的死亡结果之间的关联显示在非线性样条模型（这里是部分结果）从图中我们可以看到，不管是在全人群还是非吸烟人群中，BMI与all-cause, communicable, and non-communicable disease mortality之间呈U型关系，死亡率最低点的BMI大概是25Kg/m2。死亡风险均随着BMI先降低后增加。
-
三、RCS简介
研究数据通常来自于观察某个特定群体的个体，并收集其相关信息。这些数据可以是从临床试验、队列研究、调查问卷等途径获得的。数据的特征通常包括个体的生存时间、生存状态（例如是否死亡）、危险因素（例如年龄、性别、治疗方案等）以及其他可能影响生存的变量。
-
「限制立方样条的原理和优势」
限制立方样条是一种非参数的拟合方法，通过将连续变量（如生存时间）转化为多个分段函数来建模生存曲线。RCS在建模过程中不对数据分布作出假设，因此适用于各种类型的生存数据。与传统方法相比，RCS具有更高的灵活性，可以更准确地拟合生存曲线的形状。通过使用适当的节点数和位置，RCS可以探测到生存风险的非线性关系，并提供更精细的生存预测。
-
「生存时间和时间截断指示符的定义」
生存时间是指从某个初始事件（如诊断、手术等）到目标事件（如死亡、复发等）之间的时间间隔。在生存分析中，我们通常使用时间截断指示符来处理那些在观察期限结束时没有遇到目标事件的个体。即使没有遇到目标事件，我们仍然希望利用这些数据对生存曲线进行建模和估计。
四、RSC建模和结果
5.1 数据集载入
library(survival)
str(gbsg)
结果展示：
> str(gbsg)
'data.frame': 686 obs. of 10 variables:
$ age : int 49 55 56 45 65 48 48 37 67 45 ...
$ meno : int 0 1 1 0 1 0 0 0 1 0 ...
$ size : int 18 20 40 25 30 52 21 20 20 30 ...
$ grade : int 2 3 3 3 2 2 3 2 2 2 ...
$ nodes : int 2 16 3 1 5 11 8 9 1 1 ...
$ pgr : int 0 0 0 0 0 0 0 0 0 0 ...
$ er : int 0 0 0 4 36 0 0 0 0 0 ...
$ hormon : int 0 0 0 0 1 0 0 1 1 0 ...
$ rfstime: int 1838 403 1603 177 1855 842 293 42 564 1093 ...
$ status : Factor w/ 2 levels "0","1": 1 2 1 1 1 2 2 1 2 2 ...
age：患者年龄
meno：更年期状态（0表示未更年期，1表示已更年期）
size：肿瘤大小
grade：肿瘤分级
nodes：受累淋巴结数量
pgr：孕激素受体表达水平
er：雌激素受体表达水平
hormon：激素治疗（0表示否，1表示是）
rfstime：复发或死亡时间（以天为单位）
status：事件状态（0表示被截尾，1表示事件发生）
5.2 安装和加载依赖包
install.packages("remotes")
library(remotes)
remotes::install_github("liuqiang070488/ggrcs")
library(ggrcs)
library(survival)
library(rms)
library(ggplot2)
library(scales)
library(survminer)
5.3 检验是否符合线性关系
#例如，我们检验数据集中age与风险的非线性关系
ggcoxfunctional(Surv(rfstime, status)~age+log(age)+sqrt(age),data=gbsg)
这个函数可以进行age和生存的线性关系诊断，采用的方法是鞅残差，从下图age的1次、log和1/2次方来看，可以看出年龄和死亡之间存在非线性关系。
5.4 利用RCS曲线确认自由度
确定RCS的自由度数量需要在灵活性和过拟合之间进行权衡。自由度数量越多，模型越灵活，可以更准确地拟合生存曲线的形状。然而，如果自由度数量过多，模型可能过度适应样本数据，导致过拟合。一种常用的方法是根据数据集的大小来选择自由度数量，通常建议在3到5之间进行试验，并选择具有最佳拟合度和预测能力的自由度数量。
#1.构建模型。【节点数为3.4.5均进行了构建和比较】
fit <- cph(Surv(rfstime, status) ~ rcs(age, 3), x=TRUE, y=TRUE,data=gbsg)
fit1 <- cph(Surv(rfstime, status) ~ rcs(age, 4), x=TRUE, y=TRUE,data=gbsg)
fit2 <- cph(Surv(rfstime, status) ~ rcs(age, 5), x=TRUE, y=TRUE,data=gbsg)
fit
fit1
fit2
结果展示：
> fit
Cox Proportional Hazards Model
cph(formula = Surv(rfstime, status) ~ rcs(age, 3), data = gbsg,
x = TRUE, y = TRUE)
Model Tests Discrimination
Indexes
Obs 686 LR chi2 6.67 R2 0.010
Events 299 d.f. 2 R2(2,686)0.007
Center -1.5409 Pr(> chi2) 0.0356 R2(2,299)0.016
Score chi2 7.18 Dxy 0.047
Pr(> chi2) 0.0275
Coef S.E. Wald Z Pr(>|Z|)
age -0.0343 0.0129 -2.66 0.0077
age' 0.0374 0.0148 2.54 0.0112
> fit1
Cox Proportional Hazards Model
cph(formula = Surv(rfstime, status) ~ rcs(age, 4), data = gbsg,
x = TRUE, y = TRUE)
Model Tests Discrimination
Indexes
Obs 686 LR chi2 21.69 R2 0.031
Events 299 d.f. 3 R2(3,686)0.027
Center -3.7621 Pr(> chi2) 0.0001 R2(3,299)0.061
Score chi2 24.92 Dxy 0.135
Pr(> chi2) 0.0000
Coef S.E. Wald Z Pr(>|Z|)
age -0.0931 0.0188 -4.96 <0.0001
age' 0.2464 0.0550 4.48 <0.0001
age'' -0.7802 0.1954 -3.99 <0.0001
> fit2
Cox Proportional Hazards Model
cph(formula = Surv(rfstime, status) ~ rcs(age, 5), data = gbsg,
x = TRUE, y = TRUE)
Model Tests Discrimination
Indexes
Obs 686 LR chi2 23.01 R2 0.033
Events 299 d.f. 4 R2(4,686)0.027
Center -4.3585 Pr(> chi2) 0.0001 R2(4,299)0.062
Score chi2 27.55 Dxy 0.125
Pr(> chi2) 0.0000
Coef S.E. Wald Z Pr(>|Z|)
age -0.1103 0.0226 -4.88 <0.0001
age' 0.3864 0.1151 3.36 0.0008
age'' -1.5493 0.6874 -2.25 0.0242
age''' 1.4925 1.1025 1.35 0.1758
3节点：R² = 0.016 和 Dxy = 0.047； 4节点：R² = 0.061 和 Dxy = 0.135； 5节点：R² = 0.062 和 Dxy = 0.125； R²和Dxy越大，拟合的模型越优。因此选择4节点。
5.5 模型拟合
dd<-datadist(gbsg)
options(datadist='dd')
fit <- cph(Surv(rfstime, status) ~ rcs(age, 4), x=TRUE, y=TRUE,data=gbsg)
ggrcs(data=gbsg,fit=fit,x="age",histcol="sky blue")
六、总结
如果想了解如何使用RCS进行建模非线性关系、优化治疗方案和预测疾病进展等，请关注和私信我，我们一起讨论学习。原创不易，如果觉得写的还行的话，请留下您的赞和再看，谢谢！
参考文献：
[1] Bhaskaran K, Dos-Santos-Silva I, Leon DA, Douglas IJ, Smeeth L. Association of BMI with overall and cause-specific mortality: a population-based cohort study of 3·6 million adults in the UK. Lancet Diabetes Endocrinol. 2018;6(12):944-953. doi:10.1016/S2213-8587(18)30288-2