「转」RCS与非线性交互，怎样才是一个正确的思路

不学临床的医学生实战医学统计2023-09-02 16:04:02广东

这里是不学临床的医学生。本期推送我们关注一个新的话题，非线性交互。这个话题其实很早有朋友就提过，当时也非常感兴趣，也或多或少讨论或者问起过。最近聂师兄也发了一篇有关rcs交互项的推文，看完之后启发很大，于是便想到对于一个更为广义的非线性交互而言如何处理才算得当。

本来是想着能否在RCS交互的结果解释上有更多理解来丰富本期推文，不过时间不够加之理解不深，只能作罢。便作为下期推文的内容吧。

我们梳理一下本期推文的内容框架：01 既往文献有关非线性交互 → 02 实际应用场景下的非线性交互 → 03 实际应用场景下的分段交互 → 04 一点总结 → 05 参考内容 。

本期推文可能并不能很好回答这个问题，更多是个人的理解，欢迎大家交流和指出不对的地方。

01 既往文献有关非线性交互

既往文献我们需要借用一下聂师兄推文里提到的两篇文章，可能看过的朋友都比较熟悉了。

第一篇是发表在Circulation Research上的题为Retinol and Retinol Binding Protein 4 Levels and Cardiometabolic Disease Risk的文章。

这篇文章探索视黄醇和视黄醇结合蛋白与心脏代谢疾病风险在不同亚组之间有无差异。文中出现了如下表述：

We added the according markers and an interaction term with the main effect to the reference model.

If a significant interaction was detected, the subgroup-specific associations were presented with RCS (knots: 10th, 50th, 90th percentile). Otherwise, the interaction term was removed, and the model was rerun.

也就是说基于a*b的相乘交互是否显著来判断是否绘制亚组的RCS曲线。

然后是第二篇发表在Environmental Pollution上的题为Physical activity reduces the role of blood cadmium on depression: A cross-sectional analysis with NHANES data的文章。文中表述如下：

To explore the influence of different PA types on the relationship between blood Cd and depression, logistic regression and restricted cubic spline models were used to analyze different PA types.

In addition, we also included the multiplication of blood Cd and PA in the model to test their interaction with depression.

这篇文章的思路解释些许不同，RCS曲线其实已经在主分析的范畴，作为探索关联的辅助工具，因为他就是为了研究体力活动和重金属的联合作用。之后再看了是否存在相乘交互。

如果说第二篇文章尚且可以理解的话，第一篇文章基于相乘交互来判断是否绘制RCS曲线我觉得是存疑的。

a*b的相乘交互是基于线性尺度的，一旦在后面引入RCS，如果发现RCS里面存在非线性关联，那么就可以推翻前面线性交互的结果。

不过比较幸运的是，上面两篇文章反应出来的亚组趋势都是线性。虽然第二篇文章没有汇报详细的P值，不过看图形应该是线性。按照我的理解，RCS在这两篇文章中起到的作用仅仅是亚组剂量反应关系的可视化，并且这种可视化最好体现为线性关系，并没有涉及到非线性交互。

↑

左为第一篇文献，右为第二篇文献

a*b的相乘交互与我们传统的连续变量或者分类变量的亚组分析放在一起，这样二者才是对应的。

↑

个人认为这才是对应的

如果说我们真的要探讨非线性交互，有没有合适的方式。其实是有的，二次项其实便是一个很好的窥探这个问题的方法，在既往文献中也有报道。

那接下来分享我们第三篇发表在American journal of clinical nutrition上的题为Long-term coffee consumption, caffeine metabolism genetics, and risk of cardiovascular disease: a prospective analysis of up to 347,077 individuals and 8368 cases的文献。

关于非线性交互的表述如下。

↑

有关非线性的表述

我觉得这里面有两个非常重要的信息：1.基于既往文献背景判断咖啡摄入与CVD是非线性关系；2.首先构建高次项交互，交互项不显著再退回一次项交互。

如果说是因为既往研究提示非线性关系存在，这样非线性交互就有了先验的保障，整个研究才会更合理。

当然我们做研究很多时候是后验的，知道了结果再去调整我的假设。但是一个好的故事肯定是先验，从假设到分析到证实或者推翻假设。

02 实际应用场景下的非线性交互

这里我们用之前经常用到的analysis_mort数据集。后台回复RCS测试数据集可以获取本文使用的数据集Analysis_mort。

在我们不知道剂量反应关系的情况下，可以画RCS有个基础认知。复杂抽样我们在这就不考虑了，免得让问题更复杂。

library('rms')library('ggplot2')##导入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)
######熟悉的RCS##fit <- cph(Surv(PERyear_EXM,mort_all)~rcs(LBXTC,c(141,200,275))+RIDAGEYR+RIAGENDR+RIDRETH1,data=analysis_mort)fitanova(fit) #为后续程序设定数据环境（rms包函数很多都要做这两步）dd<- datadist(analysis_mort)options(datadist = 'dd')dd$limits$LBXTC[2]<-200fit <-update(fit)HR<- Predict(fit,LBXTC,fun=exp,ref.zero = TRUE)p1<-ggplot(HR)p1

可以看到总胆固醇和死亡是非线性关联。应该是有二次项的显著趋势，那我们就可以按照American journal of clinical nutrition文章的方式进行实操了。

↑

RCS结果

与预期一致，显示存在二次项的显著趋势。

######二次项##fitx2 <- coxph(Surv(PERyear_EXM,mort_all)~LBXTC+I(LBXTC^2)+RIDAGEYR+RIAGENDR+RIDRETH1,data=analysis_mort)fitx2

↑

二次项显著

我们接下来从高次项开始拟合交互。

######高次项存在时，应该从高次项开始拟合交互##fitx2_int2<- coxph(Surv(PERyear_EXM,mort_all)~LBXTC+I(LBXTC^2)+RIDAGEYR+RIAGENDR+RIDRETH1+LBXTC:RIAGENDR+LBXTC:RIAGENDR+I(LBXTC^2):RIAGENDR,data=analysis_mort)fitx2_int2

结果显示二次项的交互不显著。

↑

二次项交互不显著

我们退回到一次项的交互。

######二次项交互不存在回到一次##fitx2_int1<- coxph(Surv(PERyear_EXM,mort_all)~LBXTC+I(LBXTC^2)+RIDAGEYR+RIAGENDR+RIDRETH1+LBXTC:RIAGENDR+LBXTC:RIAGENDR,data=analysis_mort)fitx2_int1

结果显示一次项的交互也不显著。

↑

一次项交互不显著

做到这里，我们好像可以得出总胆固醇与性别在线性和非线性尺度上都不存在显著交互作用。

03 实际应用场景下的分段交互

本来做到上面已经结束了，只是想着要不要画亚组的RCS看一下两组的趋势都是啥样。

画完之后却发现了一些有意思的地方。

######男性亚组的RCS##fit1 <- cph(Surv(PERyear_EXM,mort_all)~rcs(LBXTC,c(141,200,275))+RIDAGEYR+RIDRETH1,data=analysis_mort_male)fit1anova(fit1) dd<- datadist(analysis_mort)options(datadist = 'dd')dd$limits$LBXTC[2]<-199fit1 <-update(fit1)HR1<- Predict(fit1,LBXTC,fun=exp,ref.zero = TRUE)p2<-ggplot(HR1)p2######女性亚组的RCS##fit2 <- cph(Surv(PERyear_EXM,mort_all)~rcs(LBXTC,c(141,200,275))+RIDAGEYR+RIDRETH1,data=analysis_mort_female)fit2anova(fit2) dd<- datadist(analysis_mort)options(datadist = 'dd')dd$limits$LBXTC[2]<-201fit2 <-update(fit2)HR2<- Predict(fit2,LBXTC,fun=exp,ref.zero = TRUE)p3<-ggplot(HR2)p3
#组合两组的RCS图p4<-ggplot()+theme_classic()+  geom_line(data=HR1, aes(LBXTC,yhat),linetype="solid",size=1,alpha = 1)+  geom_line(data=HR2, aes(LBXTC,yhat),linetype="dotdash",size=1,alpha = 1)+  geom_ribbon(data=HR1, aes(LBXTC,ymin = lower, ymax = upper),alpha = 0.2,fill="red")+  geom_ribbon(data=HR2, aes(LBXTC,ymin = lower, ymax = upper),alpha = 0.2,fill="green")+  scale_x_continuous(limits=c(100,400),breaks=seq(100,400,100))+   #设置x轴的格式  scale_y_continuous(limits=c(0,3.5),breaks=seq(0,3,0.5))+  #设置Y轴的格式  geom_hline(yintercept=1,linetype=2,size=0.25)p4

男性中总关联和非线性关联的P值都<0.001。女性中总关联和非线性关联的P值都大于0.05。

↑

不同性别存在差异

然后我们看一下比较的绘图结果。可以看到效应在不同性别组之间还是有差别的，特别是在LBXTC取值比较大的时候。

↑

绘图结果

在我们对RCS交互还不清楚的情况下，有没有其他办法来刻画这种可能的交互呢。其实是有的，2018年在bmj发表的题为Predicted lean body mass, fat mass, and all cause and cause specific mortality in men: prospective US cohort study就有提到分组的思想。

基于中位数分成两组，可以看到两边存在不同的趋势。

↑

bmj文章结果截图

我们也分组看一下中位数两边不同性别的结果。

#构建数据集analysis_mort_male_low<- analysis_mort[analysis_mort$RIAGENDR== 1 & analysis_mort$LBXTC<=200,]analysis_mort_male_high<- analysis_mort[analysis_mort$RIAGENDR== 1 & analysis_mort$LBXTC>200,]analysis_mort_female_low<- analysis_mort[analysis_mort$RIAGENDR== 2 & analysis_mort$LBXTC<=200,]analysis_mort_female_high<- analysis_mort[analysis_mort$RIAGENDR== 2 & analysis_mort$LBXTC>200,]#转换为pre SDanalysis_mort_male_low$LBXTC<-analysis_mort_male_low$LBXTC/sd(analysis_mort_male_low$LBXTC)analysis_mort_male_high$LBXTC<-analysis_mort_male_high$LBXTC/sd(analysis_mort_male_high$LBXTC)analysis_mort_female_low$LBXTC<-analysis_mort_female_low$LBXTC/sd(analysis_mort_female_low$LBXTC)analysis_mort_female_high$LBXTC<-analysis_mort_female_high$LBXTC/sd(analysis_mort_female_high$LBXTC)
#男性不同组结果fit_male_low <- coxph(Surv(PERyear_EXM,mort_all)~LBXTC+RIDAGEYR+RIDRETH1,data=analysis_mort_male_low)fit_male_lowfit_male_high <- coxph(Surv(PERyear_EXM,mort_all)~LBXTC+RIDAGEYR+RIDRETH1,data=analysis_mort_male_high)fit_male_high#女性不同组结果fit_female_low <- coxph(Surv(PERyear_EXM,mort_all)~LBXTC+RIDAGEYR+RIDRETH1,data=analysis_mort_female_low)fit_female_lowfit_female_high <- coxph(Surv(PERyear_EXM,mort_all)~LBXTC+RIDAGEYR+RIDRETH1,data=analysis_mort_female_high)fit_female_high

显示<中位数以及≥中位数男性女性效应均存在差距。

↑

HR值

那我们其实可以借鉴发表在European Journal of Preventive Cardiology题为Differences in HDL-related mortality risk between individuals with and without hypertension: a prospective cohort study in UK Biobank的思路进行分组交互看一下结果。

analysis_mort_low<- analysis_mort[analysis_mort$LBXTC<200,]analysis_mort_high<- analysis_mort[analysis_mort$LBXTC>=200,]#分段交互计算fit_low_int<- coxph(Surv(PERyear_EXM,mort_all)~LBXTC+RIDAGEYR+RIAGENDR+RIDRETH1+LBXTC:RIAGENDR,data=analysis_mort_low)fit_low_intfit_high_int<- coxph(Surv(PERyear_EXM,mort_all)~LBXTC+RIDAGEYR+RIAGENDR+RIDRETH1+LBXTC:RIAGENDR,data=analysis_mort_high)fit_high_int

显示在<200这一段存在显著交互，但是≥200这段并不显著。

↑

200分组交互结果

这和我们在RCS里面反映出来的结果好像有些不同，图里应该是≥200效应差距更明显。

我们仔细观察一下，其实可以发现男性女性效应相等的点并不在200。

我们查看预测值对应会发现，在185两者的取值是接近一致的。

↑

185对应yhat结果相近

我们可以修改截断值为185。

analysis_mort_low<- analysis_mort[analysis_mort$LBXTC<185,]analysis_mort_high<- analysis_mort[analysis_mort$LBXTC>=185,]#分段交互计算fit_low_int<- coxph(Surv(PERyear_EXM,mort_all)~LBXTC+RIDAGEYR+RIAGENDR+RIDRETH1+LBXTC:RIAGENDR,data=analysis_mort_low)fit_low_intfit_high_int<- coxph(Surv(PERyear_EXM,mort_all)~LBXTC+RIDAGEYR+RIAGENDR+RIDRETH1+LBXTC:RIAGENDR,data=analysis_mort_high)fit_high_int

结果显示在<185和≥185均存在交互。不过≥185显示出来的交互效应仍没有<185效应强。推测还是分段线性拟合和RCS拟合存在差距。

↑

185分组交互结果

04 一点总结

正如开头所说，并没有解决问题，只是提供了一些思路和见解。我之前做趋势的时候就做过非线性交互，如果说内容不是单独关注亚组趋势，很有可能会出现画蛇添足的问题，有时候结果很难解释。

做分析本就不是越多越好，结果解释也不存在一个千篇一律的法则，说服别人更是不易。如何在结果解释中做到突出重点但又不过度才是真正的艺术。

总觉得做分析一定是有一个思路逻辑在里面，环环相扣也好，铺开也罢，总是能承接的。

就总结这么多吧。谨函布臆，顺颂时绥~ 记得点赞、转发、在看哦~

05 参考内容

https://mp.weixin.qq.com/s/t3-rqNWhFlxwkinINpQwLA

Retinol and Retinol Binding Protein 4 Levels and Cardiometabolic Disease Risk

Physical activity reduces the role of blood cadmium on depression: A cross-sectional analysis with NHANES data

Long-term coffee consumption, caffeine metabolism genetics, and risk of cardiovascular disease: a prospective analysis of up to 347,077 individuals and 8368 cases

Predicted lean body mass, fat mass, and all cause and cause specific mortality in men: prospective US cohort study

Differences in HDL-related mortality risk between individuals with and without hypertension: a prospective cohort study in UK Biobank