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

学术动态

医学数据分析技能点01 亚组分析+交互作用(逻辑回归)

时间:2024-09-05 19:59:58 阅读:93


诺维医学科研官网:https://www.newboat.top 更新换版中!

bilibili:文章对应的讲解视频在此。熊大学习社 https://space.bilibili.com/475774512

微信公众号:熊大学习社、诺维之舟

公益网站,https://nwzz.xyz/ ,内有医学离线数据库、数据提取、科研神器等高质量资料库

诺维之舟AI:https://gpt4.nwzz.xyz 可在线使用GPT4|GPT3.5|Midjourney

课程相关资料:

(1)亚组分析+交互作用(逻辑回归)代码实例已全部公开。资料包放在公益网站(https://www.nwzz.xyz)。或关注公众号熊大学习社,回复med001,获取资料信息。

(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")


# 加载库
library(gtsummary)
library(tidyverse)
library(forestploter) # 绘制森林图
library(jstable)      # 于亚组分析
library(grid)          # 可视化



# 数据准备-----------

# 代码目录,对应修改
code_path <- "D:Courses02--2 亚组分析+交互作用-逻辑回归"
# 设置代码目录
setwd(code_path)
getwd

# 读取数据
data  <-  read.csv('data.csv')

# 查看数据
head(data)

data <- data %>% mutate(
 sex = factor(sex, levels = c(0, 1), labels = c("Female", "Male")),
 diabetes = factor(diabetes, levels = c(0, 1), labels = c('No', 'Yes')),
 age_group = factor(ifelse(age > 60, ">60", "<=60"), levels = c("<=60", ">60")),
 BMI_group = factor(cut(BMI, breaks = c(-Inf, 25, 30, Inf), labels = c("<25", "25-30", ">30")))) %>%
 rename(
   Sex = sex,
   Diabetes = diabetes,
   Age = age_group,
 )
 


str(data)


sub_plot <- TableSubgroupMultiGLM(formula = Hypertension ~ NHHR,
                                     var_subgroups = c("Age", "Sex", "Diabetes", "BMI_group"),
                                     data = data)
sub_plot %>% write.csv(file="sub_plot.csv")


sub_plot



# Count/Percent/P value/P for interaction,4列的空值设为空格
sub_plot[, c(2, 3, 7, 8)][is.na(sub_plot[, c(2, 3, 7, 8)])] <- " "
# 添加空白列,用于存放森林图的图形部分
sub_plot$` ` <- paste(rep(" ", nrow(sub_plot)), collapse = " ")
# Count/Point Estimate/Lower/Upper, 3列数据转换为数值型
sub_plot[, c(2,4,5,6)] <- apply(sub_plot[, c(2,4,5,6)], 2, as.numeric)
# 计算HR (95% CI),以便显示在图形中
sub_plot$"OR (95% CI)" <- ifelse(is.na(sub_plot$"OR"), "",
                                    sprintf("%.2f (%.2f to %.2f)",
                                            sub_plot$"OR", sub_plot$Lower, sub_plot$Upper))
# 中间圆点的大小,与Count关联,数值型
sub_plot$se <- as.numeric(ifelse(is.na(sub_plot$Count), " ", round(sub_plot$Count / 1500, 4)))
# Count列的空值设为空格
sub_plot[, c(2)][is.na(sub_plot[, c(2)])] <- " "
# 第一行Overall放在最后显示
sub_plot <- rbind(sub_plot[-1,],sub_plot[1,])
# Percent列重命名,加上%
sub_plot <- rename(sub_plot, 'Percent(%)' = Percent)

str(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")
sub_plot
# 森林图绘制
plot_sub <- forest(
 #选择需要用于绘图的列:Variable/Count/Percent/空白列/HR(95%CI)/P value/P for interaction
 data = sub_plot[, c(1, 2, 3, 9, 10, 7, 8)],  
 lower = sub_plot$Lower,  #置信区间下限
 upper = sub_plot$Upper,  #置信区间上限
 est = sub_plot$OR,  #点估计值
 sizes = sub_plot$se*0.2,
 ci_column = 4,  #点估计对应的列
 ref_line = 1,   #设置参考线位置
 xlim = c(0.8, 1.2),  # x轴的范围
 ticks_at = c(0.8,1,1.2),
 theme = tm
)
plot_sub

plot_sub <- plot_sub %>%
 # 指定行加粗
 edit_plot(row = c(1,4,7,10,14),col=c(1),gp = gpar(fontface = "bold")) %>%
 # 最上侧画横线
 add_border(part = c("header")) %>%
 # 最下侧画横线
 add_border(row=14) %>%
 # 最后一行灰色背景
 edit_plot(row = 14, which = "background", gp = gpar(fill = "grey"))
plot_sub

# 第一种方式:保存为图片
tiff('plot_sub.tiff',height = 1200,width = 1500,res= 200)
plot_sub
dev.off()

# 第二种方式
ggsave("plot_sub2.tiff", device='tiff', units = "cm", width = 30, height = 30, plot_sub)

三、小结

更多内容敬请关注熊大学习社,全网同名。

课程相关资料:

(1)亚组分析+交互作用(逻辑回归)代码实例已全部公开。资料包放在公益网站(https://www.nwzz.xyz)。或关注公众号熊大学习社,回复med001,获取资料信息。

(2)一对一论文指导学员免费获取学习资料,了解咨询扫客服二维码。

(3)关注熊大学习社。您的一键三连是我最大的动力。