Got it, let's tackle this step by step—you're dealing with multiple comparisons for two-way and three-way interactions in linear mixed models (LMMs) across 5 dependent variables, and need to do this in R. Here's a practical, actionable approach:

First, you have two levels of multiple testing to address:

• Within each dependent variable: Pairwise comparisons for your targeted two-way/three-way interactions
• Across dependent variables: Correcting for the fact you're testing 5 separate DVs (plus their interaction sub-tests)

Prioritize focusing only on the starred interactions you care about—this reduces the number of tests and boosts statistical power.

We'll use these tools, which are standard for LMMs and multiple comparisons:

• lme4: Fit linear mixed models
• emmeans: Extract marginal means and run pairwise comparisons for interactions
• dplyr/purrr: Streamline repetitive tasks (optional but makes code cleaner)

Let’s assume your dataset is named dat, with fixed effects A, B, C (for three-way interactions), random effect (1|subject) (for repeated measures on the same participants), and 5 DVs: dv1, dv2, dv3, dv4, dv5.

Step 1: Fit LMMs for Each Dependent Variable #

First, build a model for each DV. We'll use a list to organize them:

library(lme4)

# Define your 5 DVs
dv_list <- c("dv1", "dv2", "dv3", "dv4", "dv5")

# Fit models for each DV
models <- lapply(dv_list, function(dv) {
  lmer(formula(paste(dv, "~ A*B*C + (1|subject)")), data = dat)
})
names(models) <- dv_list

Step 2: Run Pairwise Comparisons for Targeted Interactions #

Use emmeans to generate pairwise comparisons for your starred interactions, and apply within-interaction correction first (we’ll use Holm’s method here—it’s more powerful than Bonferroni while still controlling Type I error).

library(emmeans)
library(purrr)

# Function to extract comparisons for your target interactions
get_target_comparisons <- function(model) {
  results <- list()
  
  # Example: Three-way interaction A*B*C (adjust if your star is a different interaction)
  emm_3way <- emmeans(model, ~ A*B*C, df = "kenward-roger")  # Use Kenward-Roger df for better accuracy
  results$three_way_AxBxC <- pairs(emm_3way, adjust = "holm") %>% as.data.frame()
  
  # Example: Two-way interaction A*B (add/remove based on your starred interactions)
  emm_2way_AxB <- emmeans(model, ~ A*B | C, df = "kenward-roger")
  results$two_way_AxB <- pairs(emm_2way_AxB, adjust = "holm") %>% as.data.frame()
  
  return(results)
}

# Apply to all models and flatten results into a single data frame
all_interaction_comps <- map(models, get_target_comparisons) %>%
  flatten() %>%
  bind_rows(.id = "test_label")

Step 3: Cross-Dependent Variable Correction #

Now, take all the p-values from your targeted tests and apply a global correction. FDR (Benjamini-Hochberg) is a good choice for exploratory analyses (controls false discovery rate), while Holm’s method is stricter for confirmatory work.

# Add globally corrected p-values
all_interaction_comps <- all_interaction_comps %>%
  mutate(
    p_adj_global = p.adjust(p.value, method = "fdr")  # Swap "fdr" with "holm" if needed
  )

• Be transparent: When writing up, specify exactly which corrections you used (e.g., "Holm’s correction for pairwise comparisons within each interaction, followed by FDR correction across all tests across 5 dependent variables").
• Only test meaningful interactions: If an interaction term isn’t significant in the initial LMM, you might skip pairwise comparisons—this aligns with confirmatory research practices.
• Random effect handling: emmeans automatically accounts for random effect variance in LMMs, so you don’t need to adjust for that manually.

If you prefer cleaner code, use purrr and dplyr to streamline the process:

library(tidyverse)

# Fit models
models <- dv_list %>%
  set_names() %>%
  map(~ lmer(formula(paste(.x, "~ A*B*C + (1|subject)")), data = dat))

# Extract and combine comparisons
all_comps <- models %>%
  map(get_target_comparisons) %>%
  list_rbind(names_to = "test_label") %>%
  mutate(p_adj_global = p.adjust(p.value, method = "fdr"))

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