【统计学代码包计划】箱线图

一次输出多个箱线图，自动根据参数的分布选择组间比较统计学方法

# 加载必要的库
library(ggplot2)
library(gridExtra)  # 用于将多个图合并
library(grid)

# 设置工作目录
setwd("D:/")

# 读取数据
bio <- read.csv(".csv")

# 确保诊断列是因子类型
bio$diagnosis <- as.factor(bio$diagnosis)

# 选择有显著差异的参数
parameters <- c("SSI", "cTBI", "HCR", "A2L", "A1T", "A2T")

# 筛选数据
tmp <- subset(bio, diagnosis == "0" | diagnosis == "2")

# 创建一个空列表，用于存储每个参数的p值
p_values <- c()

# 创建一个空列表，用于存储每个箱线图
plot_list <- list()

# 进行每个参数的正态性检验，并选择相应的组间比较方法
for (param in parameters) {
  # 正态性检验：Shapiro-Wilk检验
  shapiro_test <- shapiro.test(tmp[[param]])
  
  # 判断正态性
  if (shapiro_test$p.value > 0.05) {
    # 如果是正态分布，使用t检验
    test_result <- t.test(tmp[[param]] ~ tmp$diagnosis)
  } else {
    # 如果是非正态分布，使用Mann-Whitney U检验
    test_result <- wilcox.test(tmp[[param]] ~ tmp$diagnosis)
  }
  
  # 将p值添加到列表中
  p_values <- c(p_values, test_result$p.value)
  
  # 设置显著性星号
  if (test_result$p.value < 0.001) {
    star <- "***"
  } else if (test_result$p.value < 0.01) {
    star <- "**"
  } else if (test_result$p.value < 0.05) {
    star <- "*"
  } else {
    star <- "ns"  # 如果p值大于0.05，标记为"ns"
  }
  
  # 计算每个参数的最大值并设置合适的y位置
  max_y <- max(tmp[[param]], na.rm = TRUE)
  offset <- 0.1 * max_y  # 偏移量确保p值不会超出图的顶部
  
  # 绘制箱线图，并添加p值和星号标注
  p <- ggplot(tmp, aes(x = diagnosis, y = get(param), fill = diagnosis)) +
    geom_boxplot() +
    labs(title = paste("Boxplot of", param, "by Diagnosis"), x = "Diagnosis", y = param) +
    theme_minimal() +
    scale_fill_manual(values = c("lightblue", "lightgreen")) +
    theme(plot.title = element_text(hjust = 0.5)) +
    annotate("text", x = 1.5, y = max_y + offset, label = paste("p =", round(test_result$p.value, 3), star), size = 5, hjust = 0.5) +
    ylim(0, max_y + offset)  # 设置y轴的范围，避免p值标签超出图形顶部
  　　print(p)
  # 将图表添加到图表列表中
  plot_list[[param]] <- p
}

# 使用grid.arrange将所有图表合并到一个图形中，这一步总是输出有误，待调试，不要用！
png("boxplots_all_parameters.png", height = 7, width = 12, units = "in", res = 300)
grid.arrange(grobs = plot_list, ncol = 3)  # 每行3个图表
dev.off()

# 打印每个参数的p值
for (i in 1:length(parameters)) {
  cat("p-value for", parameters[i], ":", p_values[i], "\n")
}

　　

 结果如图：