在真实世界研究研究或者临床试验中，我们经常需要进行亚组分析，并计算交互作用（P for interaction）（如图1所示）。亚组（subgroup）是指临床研究中所有研究对象的一个子集。亚组分析一般是根据研究对象的基线特征进行分组，然后在每个亚组内进行统计分析，例如不同种族、不同年龄、不同性别、是否抽烟、是否合并某种疾病等。亚组分析的目的在于观察干预措施或暴露因素对结局的影响在不同特征的亚人群中是否不同。
在亚组分析中有一个重要的指标，就是交互作用（P for interaction），其意义在于判断干预措施或暴露因素和亚组因素有无交互作用。如果没有交互作用（P for interaction > 0.05），就无需再去看每个亚组里的效应值；如果存在交互作用（P for interaction < 0.05），才有理由去观察每个亚组的效应值，判断效应值在每个亚组里是否不同！
image.png
今天我们以R内置的pbc数据集为示例，介绍如何使用R语言的jstable包实现Cox回归亚组分析和交互作用效应。pbc数据集是survival包中的一个数据集，包含原发性胆汁性肝硬化（Primary Biliary Cirrhosis，PBC）患者的生存数据。本示例所用到的具体#####变量如下：
age：年龄，连续变量
sex: 性别，m=男性，f=女性
ascites：有无腹水，1=有，0=无
hepato: 有无肝肿大，1=有，0=无
edema: 有无水肿，1=有，0=无
spiders: 有无皮肤血管畸形，1=有，0=无
stage: 组织学阶段，1=Stage I，2=Stage II，3=Stage III，4=Stage IV
trt: 药物治疗1=D-青霉素，2=安慰剂
status: 终点状态，0=存活，1=移植（存活），2=死亡
time: 生存时间
1、使用library命令加载包并查看数据
library(jstable) #加载jstable包，用于亚组分析
library(survival) #加载survival包，调用pbc数据集
library(forestploter) #加载forestploter包，用于绘制森林图
library(grid) #加载grid包，用于基础绘图
str(pbc) #查看数据集
image.png
2、选择变量并将其转换为因子变量
library(tidyverse) #加载tidyverse包，用于处理数据
df <- pbc %>%
select(time, status, sex, age, ascites, hepato, edema, spiders, stage, trt) %>% #从pbc数据集中选择需要的变量
mutate(status = as.integer(status == 2),
age = ifelse(age > 40,"> 40", "<= 40"),
age = factor(age, levels = c("> 40", "<= 40")),
ascites = factor(ascites, levels = c(0, 1), labels = c("No", "Yes")),
hepato = factor(hepato, levels = c(0, 1), labels = c("No", "Yes")),
edema = factor(edema, levels = c(0, 1), labels = c("No", "Yes")),
spiders = factor(spiders, levels = c(0, 1), labels = c("No", "Yes")),
stage = factor(stage, levels = c(1, 2, 3, 4), labels = c("I", "II", "III", "IV")),
trt = factor(trt, levels = c(1, 2), labels = c("D-penicillamine", "Placebo"))
) #将分类变量转化为因子变量
str(df)
image.png
3、进行亚组分析和交互作用效应（P for interaction）
res <- TableSubgroupMultiCox(formula = Surv(time, status) ~ sex,
var_subgroups = c("age", "ascites", "hepato", "edema", "spiders", "stage", "trt"),
data = df
)
res
代码解释：这句代码使用“jstable”包中的TableSubgroupMultiCox函数进行多变量Cox生存分析，formula = Surv(time, status) ~ sex指定生存分析的模型，以sex为分组变量，var_subgroups指定将哪些变量设置为亚组。运行直接得出包括亚组分析和交互作用效应的结果。
image.png
4、绘制森林图
（1）处理数据
plot_df <- res #对res数据重新命名
plot_df[, c(2, 3, 9, 10)][is.na(plot_df[, c(2, 3, 9, 10)])] <- " " #选取第2、3、9、10列的数据在森林图中显示，并替换缺失值NA为一个空格字符
plot_df$` ` <- paste(rep(" ", nrow(plot_df)), collapse = " ") #添加空白列，用于存放森林图的图形部分
plot_df[, 4:6] <- apply(plot_df[, 4:6], 2, as.numeric) #将第4到第6列的数据转换为数值型
plot_df$"HR (95% CI)" <- ifelse(is.na(plot_df$"Point Estimate"), "",
sprintf("%.2f (%.2f to %.2f)",
plot_df$"Point Estimate", plot_df$Lower, plot_df$Upper)) #计算HR (95% CI)，以便显示在图形中
plot_df #查看数据
image.png
（2）正式绘图
p <- forest(
data = plot_df[, c(1, 2, 3, 11, 12, 9, 10)], #选择需要用于绘图的列
lower = plot_df$Lower, #置信区间下限
upper = plot_df$Upper, #置信区间上限
est = plot_df$`Point Estimate`, #点估计值
ci_column = 4, #点估计对应的列
ref_line = 1, #设置参考线位置
xlim = c(0, 5) # x轴的范围
)
plot(p) #输出森林图
image.png
