如何在 R 中向多个计划对比添加多重测试校正？

Question

我对 R 还很陌生，正在尝试使用它为一系列实验进行大量计划对比，并且正在努力将多个测试修正应用于测试。简而言之，我有一个数据框，其中包含实验的“地点”、实验完成的“时间”、应用的“治疗”以及对治疗的“反应”的信息。在文献的支持下，我主要对治疗类型之间的两种比较（参见代码）感兴趣，并希望进行有计划的对比来研究这两种反应。问题是，当我按时间阻塞时，我想为每个测试应用多重测试校正（bonferroni），但我已经被困了几天，试图找出如何做到这一点。

虽然我试图使这个示例保持简单，但我也有许多响应测量和站点，因此我的函数能够容纳所有这些非常重要。

可重现的示例：

library(stats)

Site <- c(rep("a", 20), rep("b", 20))
Time <- rep(c("w", "x", "y", "z"), times = 10)
Treatment <- c(rep("l", 8), rep("m", 8), rep("n", 8),
               rep("o", 8), rep("p", 8))
Response <- runif(40, min = 0, max = 20)
df <- as.data.frame(cbind(Site, Time, Treatment, Response))

#Setting levels for the treatments:
df$Treatment <- factor(x = df$Treatment, levels = c("l", "m", "n", "o", "p"))

# Define a function to perform a planned contrast for a given site
perform_planned_contrast <- function(data, var) {
  site <- unique(data$Site)[1]  # Get the site name
  cat("\nSite:", site, "| Variable:", var, "\n")     # Print the site name

  # Building hypothesis vectors
  c1 <- c(0, 0, 1, -1, 0) # Contrasting treatments n & o
  c2 <- c(0, 0, 0, 1, -1) # Contrasting treatments o & p
  
  # Combining these hypotheses into a matrix
  mat <- cbind(c1, c2)
  
  # Applying the hypothesis to the treatment variable
  contrasts(data$Treatment) <- mat
  
  # Build ANOVA model
  anova_result <- aov(as.formula(paste(var, "~ Treatment + Time")), data = data)
  print(summary(anova_result))
  
  # Planned Contrast
  PC <- summary(anova_result, split=list(Treatment=list("3 vs 4"=1,
                                                        "4 vs 5"= 2)))
  
  print(PC)
}

# Apply the function to each unique site
result <- lapply(split(df, df$Site), perform_planned_contrast,
                 var = "Response")

我尝试为我的“PC”对象创建一个新列，在其中对我的 p 值应用修正，但这不起作用，原因有两个。 PC$p.adjust 的输出是 0，我认为它没有考虑到我通过 lapply 函数运行的每个唯一站点。

我在打印（PC）之前将其添加到我的 Perform_planned_contrast 函数中

 PC$p.adjust <- p.adjust(PC$p.value, method = "bonferroni")

我认为实际的解决方案可能必须修改 lapply 输出，但我真的不知道如何去做。

感谢您的帮助！

Answer 1

我稍微改变了你的函数（下面的代码），通过broom::tidy()

返回结果的“整洁”版本，这更容易以编程方式使用。我已经用 tidyverse 做到了这一点，但它也可以在基本 R 中轻松完成，以及使用 

lapply()

、

do.call(rbind, ...)

、

na.omit()

 等..

library(dplyr)
result <- (split(dd, dd$Site)
    |> purrr::map_dfr(
          ~perform_planned_contrast(., var = "Response"), .id = "Site")
    |> tidyr::drop_na()  ## drop residual terms (no p-value anyway)
    |> mutate(p.adj = p.adjust(p.value, "bonferroni"))
)

---

这里我取出了打印语句并返回了一个整理后的值。

perform_planned_contrast <- function(data, var) {
  site <- unique(data$Site)[1]  # Get the site name
  ## cat("\nSite:", site, "| Variable:", var, "\n")     # Print the site name

  # Building hypothesis vectors
  c1 <- c(0, 0, 1, -1, 0) # Contrasting treatments n & o
  c2 <- c(0, 0, 0, 1, -1) # Contrasting treatments o & p
  
  # Combining these hypotheses into a matrix
  mat <- cbind(c1, c2)
  
  # Applying the hypothesis to the treatment variable
  contrasts(data$Treatment) <- mat
  
  # Build ANOVA model
  anova_result <- aov(as.formula(paste(var, "~ Treatment + Time")), data = data)
  ## print(summary(anova_result))
  
  # Planned Contrast
  PC <- summary(anova_result, split=list(Treatment=list("3 vs 4"=1,
                                                        "4 vs 5"= 2)))
  broom::tidy(PC[[1]])
}