医学数据分析技能点（01）亚组分析+交互作用（逻辑回归）
医学数据分析技能点（01）亚组分析+交互作用（逻辑回归）
诺维医学科研官网：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\\55-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）关注熊大学习社。您的一键三连是我最大的动力。