医学数据分析技能点02 亚组分析+交互作用（Cox回归）
诺维医学科研官网：https://www.newboat.top 更新换版中！
bilibili：文章对应的讲解视频在此。熊大学习社 https://space.bilibili.com/475774512
微信公众号：熊大学习社、诺维之舟
公益网站，https://nwzz.xyz/ ，内有医学离线数据库、数据提取、科研神器等高质量资料库
诺维之舟AI：https://gpt4.nwzz.xyz 可在线使用GPT4|GPT3.5|Midjourney
课程说明：
（1）亚组分析+交互作用（Cox回归）代码实例已全部公开。关注公众号熊大学习社，回复med002，获取资料信息。
（2）一对一论文指导学员免费获取学习资料，了解咨询扫客服二维码。
（3）关注熊大学习社。您的一键三连是我最大的动力。
一、效果
二、代码实例
#####公众号:熊大学习社#####
# 手动设置工作目录为代码和数据所在文件夹
# 步骤方法：点菜单栏“session”->"Set Work Directory"->"Choose Directory"
# 选择代码和数据所在文件夹即可
# 查看工作目录
getwd()
# 检测是否安装了相关的库，没有则自动安装
if(!require("gtsummary")) install.packages("gtsummary")
if(!require("tidyverse")) install.packages("tidyverse")
if(!require("forestploter")) install.packages("forestploter")
if(!require("jstable")) install.packages("jstable")
if(!require("grid")) install.packages("grid")
if(!require("survival")) install.packages("survival")
# 加载库
library(gtsummary)
library(tidyverse)
library(forestploter) # 绘制森林图
library(jstable) # 于亚组分析
library(grid) # 可视化
library(survival) # 生存分析
# 代码目录，对应修改
code_path <- "E:\\医学AI自媒体\\01 医学数据分析技能点\\医学数据分析技能点（02）亚组分析+交互作用（Cox回归）"
# 9 亚组分析-----------
# 设置代码目录
setwd(code_path)
getwd()
# 读取数据
final_data <- read.csv('final_data.csv')
# 数据准备
ds <- final_data %>%
mutate(
gender = factor(gender, levels = c("F", "M"), labels = c("Female", "Male")),
co_diabetes = factor(co_diabetes, levels = c(0, 1), labels = c("No", "Yes")),
co_hypertension = factor(co_hypertension, levels = c(0, 1), labels = c("No", "Yes")),
co_neoplasm = factor(co_neoplasm, levels = c(0, 1), labels = c("No", "Yes")),
co_COPD = factor(co_COPD, levels = c(0, 1), labels = c("No", "Yes")),
co_CA_surgery = factor(co_CA_surgery, levels = c(0, 1), labels = c("No", "Yes")),
co_GI = factor(co_GI, levels = c(0, 1), labels = c("No", "Yes")),
co_ICH = factor(co_ICH, levels = c(0, 1), labels = c("No", "Yes")),
co_bleeding = factor(co_bleeding, levels = c(0, 1), labels = c("No", "Yes")),
co_VTE = factor(co_VTE, levels = c(0, 1), labels = c("No", "Yes")),
co_CI = factor(co_CI, levels = c(0, 1), labels = c("No", "Yes")),
group = factor(group, levels = c(0, 1), labels = c("No", "Yes")),
) %>%
rename(
"age" = "admission_age",
"Gender" = "gender",
"Race" = "race",
"Weight" = "weight",
"Height" = "height_imputed",
"INR" = "inr",
"PT" = "pt",
"Diabetes" = "co_diabetes",
"Hypertension" = "co_hypertension",
"Neoplasm" = "co_neoplasm",
"COPD" = "co_COPD",
"CA" = "co_CA_surgery",
"GI" = "co_GI",
"ICH" = "co_ICH",
"Bleeding" = "co_bleeding",
"VTE" = "co_VTE",
"CI" = "co_CI"
)
str(ds)
# Age2
ds$Age[ds$age>=18 & ds$age<60] ='18-60'
ds$Age[ds$age>=60] ='>=60'
ds$Age <- factor(ds$Age, levels = c('18-60', '>=60'))
table(ds$Age,useNA = 'ifan')
cox_sub_plot <- TableSubgroupMultiCox(formula = Surv(surv_90 , status_90) ~ sofa_score,
var_subgroups = c("Age", "Gender", "Race","Diabetes", "Hypertension",
"Neoplasm", "COPD","CA", "GI","ICH", "Bleeding","VTE",
"CI"),
data = ds)
cox_sub_plot %>% write.csv(file="cox_sub_plot.csv")
# Count/Percent/P value/P for interaction，4列的空值设为空格
cox_sub_plot[, c(2, 3, 7, 8)][is.na(cox_sub_plot[, c(2, 3, 7, 8)])] <- " "
# 添加空白列，用于存放森林图的图形部分
cox_sub_plot$` ` <- paste(rep(" ", nrow(cox_sub_plot)), collapse = " ")
# Count/Point Estimate/Lower/Upper, 3列数据转换为数值型
cox_sub_plot[, c(2,4,5,6)] <- apply(cox_sub_plot[, c(2,4,5,6)], 2, as.numeric)
# 计算HR (95% CI)，以便显示在图形中
cox_sub_plot$"HR (95% CI)" <- ifelse(is.na(cox_sub_plot$"Point Estimate"), "",
sprintf("%.2f (%.2f to %.2f)",
cox_sub_plot$"Point Estimate", cox_sub_plot$Lower, cox_sub_plot$Upper))
# 中间圆点的大小，与Count关联，数值型
cox_sub_plot$se <- as.numeric(ifelse(is.na(cox_sub_plot$Count), " ", round(cox_sub_plot$Count / 1500, 4)))
# Count列的空值设为空格
cox_sub_plot[, c(2)][is.na(cox_sub_plot[, c(2)])] <- " "
# 第一行Overall放在最后显示
cox_sub_plot <- rbind(cox_sub_plot[-1,],cox_sub_plot[1,])
# Percent列重命名，加上%
cox_sub_plot <- rename(cox_sub_plot, 'Percent(%)' = Percent)
str(cox_sub_plot)
# 森林图的格式设置
tm <- forest_theme(base_size = 10,
ci_pch = 15,
#ci_col = "#762a83",
ci_col = "black",
ci_fill = "black",
ci_alpha = 0.8,
ci_lty = 1,
ci_lwd = 1.5,
ci_Theight = 0.2,
refline_lwd = 1,
refline_lty = "dashed",
refline_col = "grey20",
vertline_lwd = -0.1,
vertline_lty = "dashed",
vertline_col = "red",
summary_fill = "blue",
summary_col = "#4575b4",
footnote_cex = 0.6,
footnote_fontface = "italic",
footnote_col = "grey20")
# 森林图绘制
plot_sub <- forest(
#选择需要用于绘图的列：Variable/Count/Percent/空白列/HR(95%CI)/P value/P for interaction
data = cox_sub_plot[, c(1, 2, 3, 9, 10, 7, 8)],
lower = cox_sub_plot$Lower, #置信区间下限
upper = cox_sub_plot$Upper, #置信区间上限
est = cox_sub_plot$`Point Estimate`, #点估计值
sizes = cox_sub_plot$se*0.15,
ci_column = 4, #点估计对应的列
ref_line = 1, #设置参考线位置
xlim = c(0.5, 1.5), # x轴的范围
ticks_at = c(0.5,1,1.5),
theme = tm
)
plot_sub
plot_sub <- plot_sub %>%
# 指定行加粗
edit_plot(row = c(1,4,7,12,15,18,21,24,27,30,33,36,39,42),col=c(1),gp = gpar(fontface = "bold")) %>%
# 最上侧画横线
add_border(part = c("header")) %>%
# 最下侧画横线
add_border(row=43) %>%
# 最后一行灰色背景
edit_plot(row = 43, which = "background", gp = gpar(fill = "grey"))
plot_sub
# 第一种方式:保存为图片
tiff('plot_sub.tiff',height = 2200,width = 2000,res= 200)
plot_sub
dev.off()
# 第二种方式
ggsave("plot_sub2.tiff", device='tiff', units = "cm", width = 30, height = 30, plot_sub)
三、小结
更多内容敬请关注熊大学习社，全网同名。
课程相关资料：
（1）亚组分析+交互作用（Cox回归）代码实例已全部公开。关注公众号熊大学习社，回复med002，获取资料信息。
（2）一对一论文指导学员免费获取学习资料，了解咨询扫客服二维码。
（3）关注熊大学习社。您的一键三连是我最大的动力。