前几天演示了如何使用purrr进行回归模型的亚组分析及森林图绘制：
R语言亚组分析及森林图绘制
本来找了好久没找到可以实现这个功能的R包，都打算自己写个包了，没想到这几天找到了！
完美实现COX回归和logistic回归的亚组分析，除此之外，还支持svyglm、svycoxph的结果，并且数据结果可直接用于绘制森林图，连NA和各种空行都给你准备好了！
包的名字叫jstable，直达网址：https://jinseob2kim.github.io/jstable/
这个包除了亚组分析外，还有其他很多函数，大家自己探索即可，我这里就演示下如何进行亚组分析！
安装
代码语言：javascript复制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)
画图就非常简单！
代码语言：javascript复制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)
这样就搞定了，真的是非常简单了，省去了大量的步骤。
