临床上,因变量和临床的结局有时候不是线性关系,而回归模型有一个重要的假设就是自变量和因变量呈线性关联,因此非线性关系模型用回归分析来拟合受到限制。因此,一个更好的解决方法是拟合自变量与因变量之间的非线性关系,限制性立方样条(Restricted cubic spline,RCS)就是分析非线性关系的最常见的方法之一。
既往再文章《R语言绘制复杂抽样设计logistic回归限制立方样条图(RCS)》我已经介绍了nhanes数据(复杂调查抽样数据)怎么手动绘制限制立方样条图,但是手动相对麻烦,而且要清除异常离群值,绘制不好的话没有rms包绘制得漂亮,对新手有点难度。参考了网络上一些其他大神得画法,经过我自己改良后,svyggrcs1.8给出了两种绘图得方法,一个是常规方法,一个是重抽样绘图,重抽样绘图能很好的印证你的结果,增强结论的可靠性。因为survey包还要对残差进行加权,计算可信区间,因此我觉得重抽样绘制的图形更加可靠。总而言之,svyggrcs1.8还是一个新的函数包,还有很多待升级的空间,有意见或建议都可以私信我。
先说个结论,rms包绘图和我既往文章的绘图结果是一样的,间接证明了我文章《R语言绘制复杂抽样设计logistic回归限制立方样条图(RCS)》方法没有问题。
下面我来演示一下怎么使用svyggrcs1.8函数绘制nhanes数据(复杂调查数据)的限制立方样条图绘制。支持线性回归、逻辑回归和cox回归。
先要安装常规一些辅助包,如果没有安装的话我们要提前安装好
install.packages(c("stringr","ggplot2","rms","cowplot","survey"))
导入我写的函数,把函数放在指定的位置(我这里在E:/r/test/),然后用下面代码导入
source("E:/r/test/svyggrcs1.8.R")
正常导入的话应该生成5个函数,
我们先来做个cox回归的
先导入R包和数据
noNA<-read.csv("E:/r/test/noNA.csv",sep=',',header=TRUE)
这是一个转移性胃癌患者(Power、Capanu、Kelsen 和 Shah 2011)的数据(公众号回复:胃癌数据,可以获得数据),数据很多我们选取一部分建模,age_dx:年龄,group:分组变量,分为存活率小于2年的和大于两年的,inv_weight:概率权重,ssize:每个分组患者的人数,survival生存时间,surv_cens生存结局
建立调查函数表
dstr2 <- svydesign(id = ~1, strata = ~group, prob = ~inv_weight,
fpc = ~ssize, data = noNA)
假设我们研究的是年龄和生产率关系。生成cox回归方程,我们研究的目前是年龄age这里需要给它rcs(Age,4)
svy.cox.fit <- svycoxph(Surv(survival, surv_cens) ~ ECOG + liver_only +
Alb + Hb + rcs(Age,4) + Differentiation + Gt_1_m1site + lymph_only, x = TRUE,
design = dstr2)
生成限制立方样条图,非常简单,一句话代码
out<-svyggrcs(svy.cox.fit)
生成结果在out这里
我们打开看一下,生成3个内容,newdat是绘图数据,你可以通过它手动绘图,这样自由度大一点,boot.p是重抽样生成的图片,p2是rms包生成的图片
我们点开看一下,常规绘图
out[["p2"]]
再看下重抽样的
out[["boot.p"]]
我们可以修改X轴和Y轴标题,P值的位置也可以通过px和py参数修改
out<-svyggrcs(svy.cox.fit,xlab = "血糖",ylab="预测值",ggtitle="血糖与预测值关系")
再看看提取的数据
dat<-out[["newdat"]]
手动绘制灵活一点
ggplot(dat,aes(Age,yhat))+
geom_ribbon(aes(ymin=lower,ymax=upper),fill="red",alpha=.2)+
geom_line(linewidth=1,col="red",alpha=.2)+
theme_bw()+
theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())
接下来演示一下逻辑回归,差不多,直接上代码了
#逻辑回归
library(survey)
bc<-read.csv("E:/r/test/subtext.csv",sep=',',header=TRUE)
bc <- na.omit(bc)
bc$DMDMARTL<-ifelse(bc$DMDMARTL==1,1,0)
bc$RIAGENDR<-as.factor(bc$RIAGENDR)
bc$RIDRETH1<-as.factor(bc$RIDRETH1)
bc$DMDMARTL<-as.factor(bc$DMDMARTL)
#bc$factor.FVC<-as.factor(bc$factor.FVC)
bcSvy1<- svydesign(ids = ~ SDMVPSU, strata = ~ SDMVSTRA, weights = ~ WTMEC2YR,
nest=TRUE,data = bc)
svy.glm<- svyglm(factor.FVC ~ RIAGENDR + RIDAGEYR+RIDRETH1+ DMDMARTL + rcs(LBXGH,4),
design = bcSvy1)
也是一句话生成图,
out<-svyggrcs(svy.glm)
out[["p2"]]
最后展示线性回归,也是直接上代码
##线性回归
library(survey)
bc<-read.csv("E:/nhanes/nhanes.csv",sep=',',header=TRUE)
bc <- na.omit(bc)
bc$DMDMARTL<-ifelse(bc$DMDMARTL==1,1,0)
bc$RIAGENDR<-as.factor(bc$RIAGENDR)
bc$RIDRETH1<-as.factor(bc$RIDRETH1)
bc$DMDMARTL<-as.factor(bc$DMDMARTL)
bcSvy2<- svydesign(ids = ~ SDMVPSU, strata = ~ SDMVSTRA, weights = ~ WTMEC2YR,
nest=TRUE,data = bc)
svy.glm.fit <- svyglm(SPXNFVC ~ RIAGENDR + RIDAGEYR+RIDRETH1+ DMDMARTL + LBXGH,
design = bcSvy2)
svy.lm<- svyglm(SPXNFVC ~ RIAGENDR + RIDAGEYR+RIDRETH1+ DMDMARTL + rcs(LBXGH,4),
design = bcSvy2)
out<-svyggrcs(svy.lm)
out[["p2"]]
如果你还不会,可以看下面视频
Nhanes数据(复杂调查数据)绘制限制立方样条图(rcs)
需要获得svyggrcs1.8函数尝鲜版,可以看下面这篇文章