欢迎访问诺维之舟医学科研平台,我们的课程比丁香园更香!
Language:

学术动态

​医学数据分析技能点02 亚组分析+交互作用(Cox回归)

时间:2024-09-05 21:01:57 阅读:84

诺维医学科研官网: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自媒体 医学数据分析技能点医学数据分析技能点(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(01), labels = c("No""Yes")),
    co_hypertension = factor(co_hypertension, levels = c(01), labels = c("No""Yes")),
    co_neoplasm = factor(co_neoplasm, levels = c(01), labels = c("No""Yes")),
    co_COPD = factor(co_COPD, levels = c(01), labels = c("No""Yes")),
    co_CA_surgery = factor(co_CA_surgery, levels = c(01), labels = c("No""Yes")),
    co_GI = factor(co_GI, levels = c(01), labels = c("No""Yes")),
    co_ICH = factor(co_ICH, levels = c(01), labels = c("No""Yes")),
    co_bleeding = factor(co_bleeding, levels = c(01), labels = c("No""Yes")),
    co_VTE = factor(co_VTE, levels = c(01), labels = c("No""Yes")),
    co_CI = factor(co_CI, levels = c(01), labels = c("No""Yes")),
    group = factor(group, levels = c(01), 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(2378)][is.na(cox_sub_plot[, c(2378)])] <- " "  
# 添加空白列,用于存放森林图的图形部分
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 / 15004)))
# 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(12391078)],  
  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.51.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)关注熊大学习社。您的一键三连是我最大的动力。