在GGally的ggpairs相关性绘图中实现多重检验显著性水平校正

Great question—this is a common gotcha with ggpairs since its default behavior doesn’t handle multiple testing adjustments automatically. Let’s tackle your two questions one by one:

Short answer: Yes, that p.adjust argument in aes() is being ignored. #

ggplot2 aesthetics (aes()) are designed for mapping data to visual properties (like color, size, or position), not statistical test parameters. ggpairs doesn’t recognize p.adjust as a valid aesthetic, so it simply ignores the argument entirely.

How to verify this & view raw/corrected p-values: #

You can extract the raw p-values from the ggpairs output and compare them to manually corrected values. Then, replace the default significance stars with actual p-values using a custom plotting function:

library(GGally)
library(ggplot2)
data(iris)

# First, run your original ggpairs call
gg <- ggpairs(iris, columns = 1:4, aes(p.adjust = "Bonferroni"))

# Extract raw p-values from the plot object (example for upper triangle plots)
raw_p_values <- lapply(gg$plots, function(plot) {
  if (!is.null(plot$data$p.value)) plot$data$p.value else NA
})

# Manually calculate Bonferroni-corrected p-values for comparison
corrected_p <- p.adjust(unlist(raw_p_values[!is.na(raw_p_values)]), method = "bonferroni")

# To display actual p-values (instead of stars), use a custom correlation function
my_cor_with_p <- function(data, mapping, method = "pearson", adjust = "bonferroni") {
  x <- eval_data_col(data, mapping$x)
  y <- eval_data_col(data, mapping$y)
  
  # Run correlation test and adjust p-value
  test_result <- cor.test(x, y, method = method)
  adj_p <- p.adjust(test_result$p.value, method = adjust)
  
  # Create text to display
  display_text <- paste0(
    "r = ", round(test_result$estimate, 2), "\n",
    "p (adj) = ", round(adj_p, 4)
  )
  
  ggplot(data, mapping) +
    geom_text(aes(label = display_text), size = 3) +
    theme_void()
}

# Use the custom function in ggpairs
ggpairs(iris, columns = 1:4, upper = list(continuous = my_cor_with_p))

This will show you both the correlation coefficient and the corrected p-value directly in each plot cell, confirming that the original p.adjust argument had no effect.

When working with grouped data (e.g., splitting by species), you need to account for all correlation tests across all groups when applying corrections. The default ggpairs workflow doesn’t handle this, so you’ll need to precompute all p-values, apply your chosen correction, then feed those corrected values back into a custom plotting function.

Step-by-step solution: #

library(GGally)
library(ggplot2)
library(dplyr)
library(tidyr)
data(iris)

# 1. Precompute all correlation results (raw p-values + coefficients) across groups
cor_summary <- iris %>%
  pivot_longer(-Species, names_to = "var1", values_to = "val1") %>%
  pivot_longer(-c(Species, var1), names_to = "var2", values_to = "val2") %>%
  filter(var1 < var2) %>%  # Avoid duplicate variable pairs
  group_by(Species, var1, var2) %>%
  summarise(
    r = cor(val1, val2, method = "pearson"),
    raw_p = cor.test(val1, val2)$p.value,
    .groups = "drop"
  )

# 2. Apply multiple testing corrections (FDR, Bonferroni, etc.)
# Choose: correct across ALL tests, or per group
cor_summary <- cor_summary %>%
  # Correct across all tests in the dataset
  mutate(
    p_bonferroni = p.adjust(raw_p, method = "bonferroni"),
    p_fdr = p.adjust(raw_p, method = "fdr")
  ) %>%
  # Optional: correct within each species group
  group_by(Species) %>%
  mutate(
    p_bonferroni_group = p.adjust(raw_p, method = "bonferroni"),
    p_fdr_group = p.adjust(raw_p, method = "fdr")
  ) %>%
  ungroup()

# 3. Custom function to plot corrected results in ggpairs
grouped_cor_with_adj <- function(data, mapping, adjust_col = "p_fdr") {
  # Get current group and variable pair
  current_group <- unique(data$Species)
  x_var <- as_label(mapping$x)
  y_var <- as_label(mapping$y)
  
  # Fetch precomputed corrected results
  result <- cor_summary %>%
    filter(
      Species == current_group,
      (var1 == x_var & var2 == y_var) | (var1 == y_var & var2 == x_var)
    )
  
  # Add significance stars (optional)
  sig_star <- case_when(
    result[[adjust_col]] < 0.001 ~ "***",
    result[[adjust_col]] < 0.01 ~ "**",
    result[[adjust_col]] < 0.05 ~ "*",
    TRUE ~ ""
  )
  
  display_text <- paste0(
    "r = ", round(result$r, 2), "\n",
    "p (adj) = ", round(result[[adjust_col]], 4), sig_star
  )
  
  ggplot(data, mapping) +
    geom_point(alpha = 0.5) +
    geom_text(aes(label = display_text), 
              x = mean(data[[x_var]]), y = mean(data[[y_var]]),
              size = 3) +
    theme_bw()
}

# 4. Generate the final grouped ggpairs plot
ggpairs(iris, columns = 1:4, aes(color = Species),
        upper = list(continuous = grouped_cor_with_adj),
        lower = list(continuous = wrap("points", alpha = 0.5)))

This approach gives you full control over the correction method and scope (global vs. per-group), and ensures all tests are accounted for in the adjustment.

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