Nhanes数据（复杂调查数据）绘制限制立方样条（rcs）函数svyggrcs1.8尝鲜版发布

news2024/9/21 0:29:17

临床上，因变量和临床的结局有时候不是线性关系，而回归模型有一个重要的假设就是自变量和因变量呈线性关联，因此非线性关系模型用回归分析来拟合受到限制。因此，一个更好的解决方法是拟合自变量与因变量之间的非线性关系，限制性立方样条(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"]]