限制性立方样条 RCS
限制性立方样条 RCS
# * 1.生成数据 ------------------------------------------------------------
## 加载所需要的包
library(ggplot2)
#install.packages('rms')
library(rms)
## 生成数据
n <- 1000
set.seed(731)
age <- 50 + 12*rnorm(n)
label(age) <- "age"
sex <- factor(sample(c("Male","Female"),n, rep=T,prob = c(6,4)))
cens <- 15*runif(n)
h <- 0.02*exp(0.04*(age-50) + 0.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)
head(data )
# * 2.做回归、拟合曲线 ------------------------------------------------------------
# 对数据进行打包，整理
dd <- datadist(data) #为后续程序设定数据环境
options(datadist='dd') #为后续程序设定数据环境
# 拟合模型
fit<- cph(Surv(time,death) ~ rcs(age,4) + sex,data=data) # 节点数设为4
# 非线性检验
# P<0.05为存在非线性关系
anova(fit)
# 查看HR预测表
# 看一下预测的HR所对因的age
HR<-Predict(fit, age,fun=exp,ref.zero = TRUE)
head(HR)
# * 3.绘图 ------------------------------------------------------------
## 绘图
### 无分组
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=48.89330,size=1,color = '#d40e8c')+ #查表HR=1对应的age
labs(title = "Stroke Risk", x="Age", y="HR (95%CI)")
### 性别分组
HR1 <- Predict(fit, age, sex=c('Male','Female'),
fun=exp,type="predictions",
ref.zero=TRUE,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 = "RCS", x="Age", y="HR (95%CI)")
# * 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 = T)
