Jstable 包另一个举例

library(jstable) #加载jstable包，用于亚组分析
library(survival) #加载survival包，调用pbc数据集
library(forestploter) #加载forestploter包，用于绘制森林图
library(grid) #加载grid包，用于基础绘图
str(pbc) #查看数据集


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)



res <- TableSubgroupMultiCox(formula = Surv(time, status) ~ sex,
                             var_subgroups = c("age", "ascites", "hepato", "edema", "spiders", "stage", "trt"), 
                             data = df
)
res


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  #查看数据


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)  #输出森林图