如何将裂区分析的R代码改写为循环化可复现脚本？

问题背景 #

现有裂区分析R代码需手动修改响应变量名、图表轴标签、输出文件名，且需手动复制Anova与LSD检验结果，面对100+变量时效率极低，需实现全流程自动化。

自动化实现方案 #

1. 前期准备：加载包与数据预处理 #

# 加载所需工具包
library(agricolae)
library(dplyr)
library(ggplot2)

# 读取数据并转换因子类型
data <- read.csv("file.csv", header = TRUE)
data$Genotype <- as.factor(data$Genotype)
data$N.level <- as.factor(data$N.level)
data$Block <- as.factor(data$Block) # 确保区组为因子类型

2. 定义自动化分析函数 #

将模型拟合、LSD检验、统计量计算、绘图、结果保存全流程封装为函数，通过参数传入响应变量信息：

split_plot_analysis <- function(response_var, y_label, output_prefix) {
  # 拟合裂区模型
  model <- sp.plot(block = data$Block, 
                   pplot = data$Genotype, 
                   splot = data$N.level, 
                   Y = data[[response_var]])
  
  # 自动提取误差自由度与均方
  Edf_a <- model$gl.a
  Edf_b <- model$gl.b
  EMS_a <- model$Ea
  EMS_b <- model$Eb
  
  # 执行三次LSD检验
  out1 <- LSD.test(y = data[[response_var]], 
                   trt = data$Genotype,
                   DFerror = Edf_a, 
                   MSerror = EMS_a,
                   alpha = 0.05,
                   group = TRUE,
                   console = FALSE)
  
  out2 <- LSD.test(y = data[[response_var]], 
                   trt = data$N.level,
                   DFerror = Edf_b, 
                   MSerror = EMS_b,
                   alpha = 0.05,
                   group = TRUE,
                   console = FALSE)
  
  out3 <- LSD.test(y = data[[response_var]], 
                   trt = interaction(data$Genotype, data$N.level, sep = ":"),
                   DFerror = Edf_b, 
                   MSerror = EMS_b,
                   alpha = 0.05,
                   group = TRUE,
                   console = FALSE)
  
  # 整理交互作用分组结果，匹配均值数据
  ascend_AB <- out3$groups %>%
    rownames_to_column("interaction") %>%
    separate(interaction, into = c("Genotype", "N.level"), sep = ":") %>%
    arrange(Genotype, N.level)
  
  # 计算均值与标准差
  MeanSD_AB <- data %>%
    group_by(Genotype, N.level) %>%
    summarise(avg_AB = mean(.data[[response_var]]),
              sd = sd(.data[[response_var]]),
              .groups = "drop") %>%
    left_join(ascend_AB, by = c("Genotype", "N.level"))
  
  # 绘制交互作用柱状图
  plot <- ggplot(MeanSD_AB, aes(x = N.level,
                                y = avg_AB,
                                fill = Genotype)) +
    geom_bar(stat = "identity", color = "black", position = position_dodge(width=0.9)) +
    geom_errorbar(aes(ymax = avg_AB + sd, ymin = avg_AB - sd), 
                  position = position_dodge(width=0.9), width = 0.25) +
    labs(title = paste(response_var, "基因型×施氮水平交互作用"),
         x = "施氮水平",
         y = y_label,
         fill = "基因型") +
    geom_text(aes(label = groups), 
              position = position_dodge(width = 0.9),
              vjust = -(0.5), size = 3) +
    theme_bw()
  
  # 自动保存所有结果文件
  write.table(model$anova, file = paste0(output_prefix, "_anova.csv"), sep = ",", row.names = TRUE)
  write.table(out1$groups, file = paste0(output_prefix, "_genotype_lsd.csv"), sep = ",", row.names = TRUE)
  write.table(out2$groups, file = paste0(output_prefix, "_nlevel_lsd.csv"), sep = ",", row.names = TRUE)
  write.table(out3$groups, file = paste0(output_prefix, "_interaction_lsd.csv"), sep = ",", row.names = TRUE)
  write.table(MeanSD_AB, file = paste0(output_prefix, "_mean_sd.csv"), sep = ",", row.names = FALSE)
  ggsave(paste0(output_prefix, "_interaction_plot.jpeg"), plot = plot, 
         width = 10, height = 10, units = "cm", dpi = 300)
  
  # 返回结果列表（可选，用于后续调试或查看）
  return(list(model = model, lsd_genotype = out1, lsd_nlevel = out2, lsd_interaction = out3, plot = plot))
}

3. 批量运行分析 #

定义需分析的响应变量列表与对应轴标签，循环调用函数完成全量分析：

# 替换为你的实际响应变量与对应标签
response_vars <- c("GPC.30DAP", "GPC.40DAP", "Yield", "Protein.Content")
y_labels <- c("GPC (%) 花后30天", "GPC (%) 花后40天", "产量 (kg/亩)", "蛋白质含量 (%)")

# 批量执行分析
for (i in seq_along(response_vars)) {
  split_plot_analysis(response_var = response_vars[i], 
                      y_label = y_labels[i], 
                      output_prefix = response_vars[i])
}

核心优化点 #

• 用data[[response_var]]实现响应变量动态调用，无需手动替换变量名
• 自动提取模型误差参数，避免手动复制
• LSD检验结果与统计量自动保存为CSV文件，省去手动复制步骤
• 图表标题、轴标签随参数自动生成，无需手动修改
• 循环批量处理所有变量，一键完成100+变量的全流程分析

内容的提问来源于stack exchange，提问作者Luqman