I am new-ish to R and am trying to use it to conduct a large amount of planned contrasts for a bunch of experiments and am struggling to apply multiple testing corrections to the test. To briefly summarize, I have a data frame containing in information for the “Site” of my experiment, the “Time” that experiment was done, the “Treatment” applied, and the “Response” to the treatment. With support from the literature, I am primarily interested in two comparisons (see code) between treatment types and want to conduct a planned contrast to investigate those two responses. The problem is, when I block by time, I would like to apply a multiple testing correction (bonferroni) for each test, but I have been stuck for a few days trying to figure it out how to do it. It seems like the lsmeans package might be able to do what I want, but I am lost at how to apply it to my code.

While I tried to keep this example simple, I also have many response measurements and sites, so it is important that my function can accommodate them all.

Reproducible Example:

library(stats)

Site <- rep("a", 20)
Time <- rep(c("w", "x", "y", "z"), times = 5)
Treatment <- c(rep("l", 4), rep("m", 4), rep("n", 4),
               rep("o", 4), rep("p", 4))
Response <- runif(20, 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) {

  # Building hypotheses 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 weights 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 site for one var
result <- lapply(split(df, df$Site), perform_planned_contrast,
                 var = "Response")

I have tried creating a new column to my “PC” object where I am applying the correction to my p values but this does not work for two reasons. The output of the PC$p.adjust is 0, and I do not think it is accounting for each unique Site I am running through my lapply function.

I added this to my perform_planned_contrast function just before print(PC)

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

I think the actual solution might have to modify the lapply output, but I am truly lost at how I would go about doing that.

Thanks for your help!