R语言亚组分析1行代码实现！

医学和生信笔记

发布于 2023-08-30 13:06:33

1.3K10

代码可运行

文章被收录于专栏：医学和生信笔记医学和生信笔记

运行总次数：0

代码可运行

前几天演示了如何使用purrr进行回归模型的亚组分析及森林图绘制：

R语言亚组分析及森林图绘制

本来找了好久没找到可以实现这个功能的R包，都打算自己写个包了，没想到这几天找到了！

完美实现COX回归和logistic回归的亚组分析，除此之外，还支持svyglm、svycoxph的结果，并且数据结果可直接用于绘制森林图，连NA和各种空行都给你准备好了！

包的名字叫jstable，直达网址：https://jinseob2kim.github.io/jstable/

这个包除了亚组分析外，还有其他很多函数，大家自己探索即可，我这里就演示下如何进行亚组分析！

安装

install.packages("jstable")

## From github: latest version
remotes::install_github('jinseob2kim/jstable')
library(jstable)

准备数据

还是使用上次演示的数据。

使用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-死亡

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

str(colon)

可以使用cox回归探索危险因素。分类变量需要变为因子型，这样在进行回归时会自动进行哑变量设置。

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

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)

亚组分析

使用jstable，只要1行代码即可！！！

library(jstable)

res <- TableSubgroupMultiCox(
  
  # 指定公式
  formula = Surv(time, status) ~ rx, 
  
  # 指定哪些变量有亚组
  var_subgroups = c("sex","age","obstruct","perfor","adhere",
                    "differ","extent","surg","node4"), 
  data = df #指定你的数据
  )
res

直接就得出了结果！除了亚组分析的各种结果，还给出了交互作用的P值！

画森林图

这个结果不需要另存为csv也能直接使用（除非你是细节控，需要修改各种大小写等信息），当然如果你需要HR(95%CI)这种信息，还是需要自己添加一下的。

我们添加个空列用于显示可信区间，并把不想显示的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)