1行代码实现:多因素回归的亚组分析,并绘制森林图

阿越就是我 医学和生信笔记 2024-04-17 09:18

图片


图片

图片

图片

图片

图片

图片

图片



多因素回归的亚组分析和森林图绘制

多因素回归的亚组分析,或者叫协变量调整的亚组分析,在平时是非常常见的操作。比如这篇文章:

下面这篇文章(https://doi.org/10.1186/s12967-022-03839-0)中就是这样的(感谢网友Qingpeng Li提供的文献素材)。

下图展示了life-essential-8这个自变量对NAFLD这个因变量的影响,并且是经过多个协变量调整的(adjusted):

图片

我在之前的协变量调整的亚组分析和森林图绘制一文中使用了非常笨的方法,也就是手动拆分数据的做法,因为没有认真研究jstable这个神包的细节,现在我发现它也是可以实现的。

下面给大家演示如何用1行代码实现多因素回归的亚组分析(仅限逻辑回归和COX回归)。

使用survival包中的colon数据集用于演示,这是一份关于结肠癌患者的生存数据,共有1858行,16列,共分为3个组,1个观察组+2个治疗组,观察他们发生终点事件的差异。

各变量的解释如下:

  • id:患者id
  • study:没啥用,所有患者都是1
  • rx:治疗方法,共3种,Obs(观察组), Lev(左旋咪唑), Lev+5FU(左旋咪唑+5-FU)
  • sex:性别,1是男性 age:年龄
  • obstruct:肠梗阻,1是有 perfor:肠穿孔,1是有
  • adhere:和附近器官粘连,1是有
  • nodes:转移的淋巴结数量
  • status:生存状态,0代表删失,1代表发生终点事件
  • differ:肿瘤分化程度,1-well,2-moderate,3-poor
  • extent:局部扩散情况,1-submucosa,2-muscle,3-serosa,4-contiguous-structures
  • surg:手术后多久了,1-long,2-short
  • node4:是否有超过4个阳性淋巴结,1代表是
  • time:生存时间
  • etype:终点事件类型,1-复发,2-死亡

分类变量需要变为因子型,这样在进行回归时会自动进行哑变量设置。

为了演示,我们只选择Obs组和Lev+5FU组的患者,所有的分类变量都变为factor,把年龄也变为分类变量并变成factor。

rm(list = ls())
library(survival)
suppressMessages(library(tidyverse))

df <- colon %>% 
  mutate(rx=as.numeric(rx)) %>% 
  filter(etype == 1, !rx == 2) %>%  #rx %in% c("Obs","Lev+5FU"), 
  select(time, status,rx, sex, age,obstruct,perfor,adhere,differ,extent,surg,node4) %>% 
  mutate(sex=factor(sex, levels=c(0,1),labels=c("female","male")),
         age=ifelse(age >65,">65","<=65"),
         age=factor(age, levels=c(">65","<=65")),
         obstruct=factor(obstruct, levels=c(0,1),labels=c("No","Yes")),
         perfor=factor(perfor, levels=c(0,1),labels=c("No","Yes")),
         adhere=factor(adhere, levels=c(0,1),labels=c("No","Yes")),
         differ=factor(differ, levels=c(1,2,3),labels=c("well","moderate","poor")),
         extent=factor(extent, levels=c(1,2,3,4),
                       labels=c("submucosa","muscle","serosa","contiguous")),
         surg=factor(surg, levels=c(0,1),labels=c("short","long")),
         node4=factor(node4, levels=c(0,1),labels=c("No","Yes")),
         rx=ifelse(rx==3,0,1),
         rx=factor(rx,levels=c(0,1))
         )

str(df)
'data.frame':   619 obs. of  12 variables:
$ time : num 968 3087 542 245 523 ...
$ status : num 1 0 1 1 1 1 0 0 0 1 ...
$ rx : Factor w/ 2 levels "0","1": 1 1 2 1 2 1 2 1 1 2 ...
$ sex : Factor w/ 2 levels "female","male": 2 2 1 1 2 1 2 1 2 2 ...
$ age : Factor w/ 2 levels ">65","<=65": 2 2 1 1 1 2 2 1 2 2 ...
$ obstruct: Factor w/ 2 levels "No","Yes": 1 1 1 2 1 1 1 1 1 1 ...
$ perfor : Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...
$ adhere : Factor w/ 2 levels "No","Yes": 1 1 2 1 1 1 1 1 1 1 ...
$ differ : Factor w/ 3 levels "well","moderate",..: 2 2 2 2 2 2 2 2 3 2 ...
$ extent : Factor w/ 4 levels "submucosa","muscle",..: 3 3 2 3 3 3 3 3 3 3 ...
$ surg : Factor w/ 2 levels "short","long": 1 1 1 2 2 1 1 2 2 1 ...
$ node4 : Factor w/ 2 levels "No","Yes": 2 1 2 2 2 2 1 1 1 1 ...

为了演示,还是使用和协变量调整的亚组分析和森林图绘制一文中一样的变量进行演示。

还是以rxageadhere这3个变量为例。

假如rx是我主要想要研究的因素,ageadhere是协变量(或者叫混杂因素),反正后面这2个变量不是我的主要研究因素,但是又对我的主要研究因素有影响,于是就需要对这俩进行调整,大概就这意思吧。。。

jstable中,1行代码即可实现这么复杂的操作:

library(jstable)

res <- TableSubgroupMultiCox(
  
  # 指定公式
  formula = Surv(time, status) ~ rx, 
  
  # 指定哪些变量有亚组
  var_subgroups = c("sex","age","obstruct","perfor","adhere",
                    "differ","extent","surg","node4"), 
  var_cov = c("age","adhere"), # 就是在这里添加你的其他协变量!!!
  data = df #指定你的数据
)

这个结果和我们手动实现的结果是一模一样的!

那么画森林图就很简单了:

我们添加个空列用于显示可信区间,并把不想显示的NA去掉即可,还需要把P值,可信区间这些列变为数值型。

plot_df <- res
plot_df[,c(2,3,9)][is.na(plot_df[,c(2,3,9)])] <- " "
plot_df$` ` <- paste(rep(" ", nrow(plot_df)), collapse = " ")
plot_df[,4:6] <- apply(plot_df[,4:6],2,as.numeric)

画图就非常简单!

library(forestploter)
library(grid)

p <- forest(
  data = plot_df[,c(1,2,3,11,9)],
  lower = plot_df$Lower,
  upper = plot_df$Upper,
  est = plot_df$`Point Estimate`,
  ci_column = 4,
  #sizes = (plot_df$estimate+0.001)*0.3, 
  ref_line = 1
  xlim = c(0.1,4)
  )
print(p)
图片

这样就搞定了,真的是非常简单,这个就是多因素回归的亚组分析森林图了,或者叫协变量调整的亚组分析森林图。

这个forestploter在使用的时候有很多细节,我在之前的推文中强调过很多遍了,如果你的森林图总是出现:置信区间不显示,中间的方块没有,缺少某一列,想增加某列等等,赶紧去看前面几篇推文,公众号后台回复森林图即可获取所有相关合集!


森林图 · 目录
上一篇R语言亚组分析及森林图绘制手册

微信扫一扫
关注该公众号