「转」NHANES复杂抽样RCS-限制性立方样条绘制

实战医学统计2023-04-26 16:50:48广东

转载来源：月明学习小铺 R中复杂抽样RCS曲线绘制完整解决方案 http://mp.weixin.qq.com/s?__biz=MzAxNjM1ODI4OQ==&mid=2247485301&idx=1&sn=a25be3451d0229a3ce4c3efbbdd98029&scene=45#wechat_redirect

这里是不学临床的医学生。本期算是一个复杂抽样RCS绘制的完整解决方案，解决了上次提到的置信区间的计算问题，同时也解决了非线性的检验问题。在推文的开头非常感谢康宁的帮助，没有他的帮助可能我们本期推文就流产了，合作非常愉快，在交流和探讨中真的学到了很多，在思路卡壳的时候一交流就会发现原来可以这样做，整个解决方案也就是由此而来了。这应该算是真正的R语言实战了。

首先还是复习一下之前推文的预设标准：复杂抽样和不考虑复杂抽样在同样的情景下（包括自变量等）RCS曲线方向基本一致。很难去想象因为复杂抽样会直接扭曲自变量和因变量的剂量关系。

然后交代一下推文的内容架构，以及付费设置的位置。我们很早提过，会做成付费推文，付费推文收益的一半我会转给康宁作为酬谢。内容架构：01 cph协变量调整时Predict和predict的结果对比 → 02 predict在coxph和svycoxph的结果对比 → 03 复杂抽样RCS曲线方案呈现（付费）→ 04 写在最后

因为是第一篇付费推文，需要补充一下，知识付费并不是为了构建知识壁垒。而是用一个相对厚道的价格给予内容生产者一些长期生存的力量。我们设置的价格为6元，可以赠送给一位朋友。

01 cph协变量调整时Predict和predict的结果对比

这部分论证是为了解决之前遗留的问题，协变量调整时两个函数算出的结果不一致。

不一致的关键在于predict未进行协变量的统一设定。这里也是康宁在后台告诉我的。可以看下面的论证，使用的数据集仍然是我们之前构建的analysis_mort，后台回复RCS测试数据集可获取。

首先我们需要为predict设定需要用来计算预测值的数据集，analysis_mort_1：自变量LBXTC使用原数据集的LBXTC，RIDAGEYR设定为49，RIAGENDR设定为1，RIDRETH1设定为3。

library('rms')##导入NHANES测试数据集analysis_mort <- read.csv("D:/Study/公众号运营/公众号素材/RCS/Analysis_mort.csv",header=T)#为了避免其他变量的缺失影响后续数据环境的设定，我们只保留需要的那部分myvars <- c("SEQN", "RIAGENDR", "RIDAGEYR", "RIDRETH1", "WTMEC4YR", "SDMVPSU", "SDMVSTRA","LBXTC", "PERyear_EXM", "mort_all")analysis_mort <-analysis_mort[myvars]##分类变量改为因子analysis_mort$RIAGENDR <- as.factor(analysis_mort$RIAGENDR)analysis_mort$RIDRETH1 <- as.factor(analysis_mort$RIDRETH1)

##有关Predict和predict的简单论证：加入协变量analysis_mort_1<-data.frame(  LBXTC=analysis_mort$LBXTC,  RIDAGEYR=49,  RIAGENDR=1,  RIDRETH1=3)##分类变量改为因子analysis_mort_1$RIAGENDR <- as.factor(analysis_mort_1$RIAGENDR)analysis_mort_1$RIDRETH1 <- as.factor(analysis_mort_1$RIDRETH1)

然后针对cph构建的fit，分别用Predict和predict计算结果。

#为后续程序设定数据环境（rms包函数很多都要做这两步）dd<- datadist(analysis_mort)options(datadist = 'dd') #cph默认就是中位数作为参照，论证过程我们就不更新下面的fit了
fit <- cph(Surv(PERyear_EXM,mort_all)~rcs(LBXTC,c(141,200,275))+RIDAGEYR+RIAGENDR+RIDRETH1,data=analysis_mort)HR1<- Predict(fit,LBXTC=analysis_mort$LBXTC,ref.zero = TRUE,fun=exp)p2<- as.data.frame(predict(fit,se=TRUE,newdata=analysis_mort_1,ref.zero = TRUE)) ##设定newdata为analysis_mort_1，利用analysis_mort_1的参数计算预测值HR_est<-exp(p2$linear.predictors) ##计算点估计lowlimit<-exp(p2$linear.predictors-(1.96*p2$se.fit)) ##计算置信区间下限uplimit<-exp(p2$linear.predictors+(1.96*p2$se.fit))  ##计算置信区间上限HR2<-data.frame(HR=HR_est,lowlimit=lowlimit,uplimit=uplimit)

比较HR1和HR2里面的计算结果，可以看到结果基本一致，只是在置信区间4位小数后有略微差异。

↑

左为HR1（Predict计算结果），右为HR2（predict计算结果）

Predict里面也封装了predict函数，在某种程度上我们用predict来进行函数预测值是可行的，当然由于predict本身也是泛型函数，predict.cph和predict.svycoxph是不是完全一致也算是解决方案一个小的缺陷点。rms包里面函数经常被重写，非常复杂。

02 predict在coxph和svycoxph的结果对比

这部分是为了解决有朋友提到的predict计算预测值时未考虑到复杂抽样，我们用predict进行计算然后比较一下就可以知道。predict是泛型函数，是存在predict.svycoxph的。

同样先设定预测值数据集。

##有关predict的简单论证：coxph和svycoxph。证明还是不一样，还是有复杂抽样的因素在里面，不管predict最后传递给谁了，复杂抽样应该还是考虑了analysis_mort_1<-data.frame(  LBXTC=analysis_mort$LBXTC,  RIDAGEYR=49,  RIAGENDR=1,  RIDRETH1=3)##分类变量改为因子analysis_mort_1$RIAGENDR <- as.factor(analysis_mort_1$RIAGENDR)analysis_mort_1$RIDRETH1 <- as.factor(analysis_mort_1$RIDRETH1)

计算predict在coxph的计算结果。

##coxph计算结果fitcoxph <- coxph(Surv(PERyear_EXM,mort_all)~rcs(LBXTC,c(141,200,275))+RIDAGEYR+RIAGENDR+RIDRETH1,data=analysis_mort)lnr1_est<- as.data.frame(predict(fitcoxph,se.fit= TRUE,newdata=analysis_mort_1,ref.zero = TRUE,type=c("lp"))) ##不能用大写的Predict,Predict是rms包的函数HR_est1<-exp(lnr1_est$fit)lowlimit1<-exp(lnr1_est$fit-(1.96*lnr1_est$se.fit))uplimit1<-exp(lnr1_est$fit+(1.96*lnr1_est$se.fit))HR1<-data.frame(HR=HR_est1,lowlimit=lowlimit1,uplimit=uplimit1)View(HR1)

计算predict在svycoxph的计算结果。

library('survey')##构建反映复杂抽样的变量NHANES_survey <- svydesign(data=analysis_mort, id=~SDMVPSU, strata=~SDMVSTRA, weights=~WTMEC4YR, nest=TRUE)##svycoxph计算结果fitsvycoxph <- svycoxph(Surv(PERyear_EXM,mort_all)~rcs(LBXTC,c(141,200,275))+RIDAGEYR+RIAGENDR+RIDRETH1,design=NHANES_survey,data=analysis_mort)lnr2_est<- as.data.frame(predict(fitsvycoxph,se= TRUE,newdata=analysis_mort_1,ref.zero = TRUE)) ##不能用大写的Predict,Predict是rms包的函数HR_est2<-exp(lnr2_est$fit)lowlimit2<-exp(lnr2_est$fit-(1.96*lnr2_est$se.fit))uplimit2<-exp(lnr2_est$fit+(1.96*lnr2_est$se.fit))HR2<-data.frame(HR=HR_est2,lowlimit=lowlimit2,uplimit=uplimit2)View(HR2)

比较两者结果，并不一致。可以看到右边的置信区间更宽，这和我们的预期是一致的，复杂抽样会有更宽的置信区间。

↑

左为HR1（coxph计算结果），右为HR2（svycoxph计算结果）

03 复杂抽样RCS曲线方案呈现（付费）

整个方案真的很不容易，也多亏了康宁的帮助，思路扩展，错误指出以及代码优化等等。为什么我们一直执着于要自己做，就是希望除了授人以鱼，更希望授人以渔。学会了解决方法和思路，大家也可以更好的去迁移和面对未来的问题。

而且需要说明的是整个方案其实并不难操作。没有用到自编函数，都是调用已经存在的包和函数。