限制性立方樣條 RCS
限制性立方樣條 RCS
# *?1.生成數(shù)據(jù) ------------------------------------------------------------
## 加載所需要的包
library(ggplot2)
#install.packages('rms')
library(rms)??
## 生成數(shù)據(jù)
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.做回歸、擬合曲線 ------------------------------------------------------------
# 對(duì)數(shù)據(jù)進(jìn)行打包，整理
dd <- datadist(data) #為后續(xù)程序設(shè)定數(shù)據(jù)環(huán)境
options(datadist='dd') #為后續(xù)程序設(shè)定數(shù)據(jù)環(huán)境
# 擬合模型
fit<- cph(Surv(time,death) ~ rcs(age,4) + sex,data=data)?# 節(jié)點(diǎn)數(shù)設(shè)為4
# 非線性檢驗(yàn)
# P<0.05為存在非線性關(guān)系
anova(fit)
# 查看HR預(yù)測(cè)表
# 看一下預(yù)測(cè)的HR所對(duì)因的age
HR<-Predict(fit, age,fun=exp,ref.zero = TRUE)
head(HR)
# * 3.繪圖 ------------------------------------------------------------
## 繪圖
### 無(wú)分組
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對(duì)應(yīng)的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)
