Inv4m Per-Plant Mineral Content (lm contrasts via emmeans)

Overview

Variant of PSU2022_ionome_content.Rmd where figure brackets come from linear-model contrasts (emmeans::emmeans + pairs) instead of marginal Welch t-tests. The brackets visualise the same kind of t-statistic that underlies each lm coefficient — same data, same model — so the asterisks in the figures match what the lm effect tables report, with full power from the joint fit (Treatment absorbed as covariate, full residual df).

Three convention notes:

• For Treatment +P vs −P brackets in the main figure: the contrast is pairs(emmeans(m, ~ Treatment | Genotype)). CTRL and Inv4m facets get facet-specific p-values (handles the S_stover concentration × Genotype interaction correctly, for example).
• For Inv4m vs CTRL brackets in the supp figure: pooled-across-treatment marginal contrast pairs(emmeans(m, ~ Genotype)).
• FDR is applied within each metric (concentration / content / HI).

Load libraries

library(dplyr)
library(tidyr)
library(ggplot2)
library(tibble)
library(ggpubr)
library(ggbeeswarm)
library(ggfx)
library(rstatix)
library(ggtext)
library(ggh4x)
library(emmeans)

Load spatially corrected data

ionome <- read.csv(
  file.path(
    paths$intermediate,
    "PSU2022_spatially_corrected_both_treatments_with_ionome.csv"
  )
) %>%
  rename(Genotype = genotype) %>%
  filter(Genotype %in% c("CTRL", "Inv4m")) %>%
  droplevels()

ionome$Treatment <- factor(ionome$Treatment, levels = c("+P", "-P"))
ionome$Genotype  <- factor(ionome$Genotype,  levels = c("CTRL", "Inv4m"))

reliable <- c("Ca", "Fe", "K", "Mg", "Mn", "P", "S", "Zn")

content_vars <- c(
  paste0(reliable, "_grain_mg"),
  paste0(reliable, "_stover_mg"),
  paste0(reliable, "_HI")
)

conc_vars <- c(
  paste0(reliable, "_grain"),
  paste0(reliable, "_stover")
)

ratio_vars <- paste0(reliable, "_ratio")

content <- ionome

cat("Sample sizes:\n")

## Sample sizes:

ionome %>% count(Genotype, Treatment) %>% print()

##   Genotype Treatment  n
## 1     CTRL        +P 16
## 2     CTRL        -P 16
## 3    Inv4m        +P 16
## 4    Inv4m        -P 16

Long-format inputs (used by both contrasts and plots)

to_plot_content <- content %>%
  select(Genotype, Treatment,
         all_of(paste0(reliable, "_grain_mg")),
         all_of(paste0(reliable, "_stover_mg"))) %>%
  pivot_longer(cols = -c(Genotype, Treatment),
               names_to = "mineral_pool", values_to = "mg") %>%
  mutate(
    mineral = sub("_(grain|stover)_mg$", "", mineral_pool),
    pool    = sub("^.+_(grain|stover)_mg$", "\\1", mineral_pool)
  ) %>%
  select(-mineral_pool) %>%
  filter(!is.na(mg)) %>%
  mutate(
    pool      = factor(pool, levels = c("stover", "grain")),
    Genotype  = gsub("Inv4m", "<i>Inv4m</i>", Genotype),
    Genotype  = factor(Genotype, levels = c("CTRL", "<i>Inv4m</i>")),
    Treatment = factor(Treatment, levels = c("+P", "-P"))
  )

to_plot_conc <- content %>%
  select(Genotype, Treatment,
         all_of(paste0(reliable, "_grain")),
         all_of(paste0(reliable, "_stover"))) %>%
  pivot_longer(cols = -c(Genotype, Treatment),
               names_to = "mineral_pool", values_to = "ppm") %>%
  mutate(
    mineral = sub("_(grain|stover)$", "", mineral_pool),
    pool    = sub("^.+_(grain|stover)$", "\\1", mineral_pool)
  ) %>%
  select(-mineral_pool) %>%
  filter(!is.na(ppm)) %>%
  mutate(
    pool      = factor(pool, levels = c("stover", "grain")),
    Genotype  = gsub("Inv4m", "<i>Inv4m</i>", Genotype),
    Genotype  = factor(Genotype, levels = c("CTRL", "<i>Inv4m</i>")),
    Treatment = factor(Treatment, levels = c("+P", "-P"))
  )

to_plot_hi <- content %>%
  select(Genotype, Treatment, all_of(paste0(reliable, "_HI"))) %>%
  pivot_longer(cols = -c(Genotype, Treatment),
               names_to = "mineral", values_to = "hi") %>%
  mutate(
    mineral   = sub("_HI$", "", mineral),
    Genotype  = gsub("Inv4m", "<i>Inv4m</i>", Genotype),
    Genotype  = factor(Genotype, levels = c("CTRL", "<i>Inv4m</i>")),
    Treatment = factor(Treatment, levels = c("+P", "-P"))
  ) %>%
  filter(!is.na(hi), is.finite(hi))

to_plot_ratio <- content %>%
  select(Genotype, Treatment, all_of(paste0(reliable, "_ratio"))) %>%
  pivot_longer(cols = -c(Genotype, Treatment),
               names_to = "mineral", values_to = "ratio") %>%
  mutate(
    mineral   = sub("_ratio$", "", mineral),
    Genotype  = gsub("Inv4m", "<i>Inv4m</i>", Genotype),
    Genotype  = factor(Genotype, levels = c("CTRL", "<i>Inv4m</i>")),
    Treatment = factor(Treatment, levels = c("+P", "-P"))
  ) %>%
  filter(!is.na(ratio), is.finite(ratio))

Fit lms per response (kept for emmeans)

# Fit lm per response and keep model objects. Same model spec everywhere.
# Results are cached to RDS; delete the RDS to force a refit.
fit_per_response <- function(responses, data) {
  models <- lapply(responses, function(resp) {
    if (!resp %in% colnames(data)) return(NULL)
    if (sum(!is.na(data[[resp]])) < 20) return(NULL)
    tryCatch(
      lm(as.formula(paste(resp, "~ Genotype * Treatment")), data = data),
      error = function(e) NULL
    )
  })
  names(models) <- responses
  models[!sapply(models, is.null)]
}

models_rds <- file.path(paths$intermediate, "ionome_content_models.rds")
if (file.exists(models_rds)) {
  cat("Loading cached lm models from", models_rds, "\n")
  .ml <- readRDS(models_rds)
  models_conc    <- .ml$conc
  models_content <- .ml$content
  models_hi      <- .ml$hi
  models_ratio   <- .ml$ratio
  rm(.ml)
} else {
  models_conc    <- fit_per_response(conc_vars, content)
  models_content <- fit_per_response(paste0(reliable, c(rep("_grain_mg", 8),
                                                       rep("_stover_mg", 8))),
                                     content)
  models_hi      <- fit_per_response(paste0(reliable, "_HI"), content)
  models_ratio   <- fit_per_response(ratio_vars, content)
  saveRDS(list(conc = models_conc, content = models_content,
               hi = models_hi, ratio = models_ratio), models_rds)
  cat("Saved lm models to", models_rds, "\n")
}

## Loading cached lm models from /Users/fvrodriguez/Library/CloudStorage/GoogleDrive-frodrig4@ncsu.edu/My Drive/repos/inv4m/results/phosphorus_paper/intermediate/ionome_content_models.rds

cat("Fitted lms — conc:", length(models_conc),
    "content:", length(models_content),
    "HI:", length(models_hi),
    "ratio:", length(models_ratio), "\n")

## Fitted lms — conc: 16 content: 16 HI: 8 ratio: 8

Effect table (lm coefficients, FDR across all)

extract_coefs <- function(models, family) {
  lapply(names(models), function(resp) {
    m <- models[[resp]]
    out <- coef(summary(m)) %>%
      as.data.frame() %>%
      rownames_to_column("predictor")
    colnames(out)[3] <- "Std.Error"
    colnames(out)[5] <- "p.value"
    out$response <- resp
    out$family   <- family
    out %>% select(family, response, everything())
  }) %>% bind_rows()
}

effects <- bind_rows(
  extract_coefs(models_conc,    "conc"),
  extract_coefs(models_content, "content"),
  extract_coefs(models_hi,      "hi"),
  extract_coefs(models_ratio,   "ratio")
) %>%
  mutate(p.adjust = p.adjust(p.value, method = "fdr")) %>%
  arrange(p.adjust)

cat("Total models fitted:", length(unique(effects$response)), "\n")

## Total models fitted: 48

cat("Total effects tested:", nrow(effects), "\n")

## Total effects tested: 192

cat("\n=== Significant effects (FDR < 0.05, no Intercept) ===\n")

## 
## === Significant effects (FDR < 0.05, no Intercept) ===

effects %>%
  filter(p.adjust < 0.05, !grepl("Intercept", predictor)) %>%
  select(family, response, predictor, Estimate, p.value, p.adjust) %>%
  as_tibble() %>%
  print(n = Inf)

## # A tibble: 31 × 6
##    family  response     predictor                    Estimate  p.value p.adjust
##    <chr>   <chr>        <chr>                           <dbl>    <dbl>    <dbl>
##  1 conc    P_stover     Treatment-P               -1592.      2.24e-26 1.65e-25
##  2 ratio   P_ratio      Treatment-P                   2.09    3.45e-17 1.62e-16
##  3 content P_stover_mg  Treatment-P                -186.      4.40e-17 2.01e-16
##  4 hi      P_HI         Treatment-P                   0.280   4.10e-12 1.68e-11
##  5 conc    P_grain      Treatment-P                -672.      4.02e- 9 1.61e- 8
##  6 conc    Zn_stover    Treatment-P                   6.85    8.78e- 8 3.31e- 7
##  7 ratio   Zn_ratio     Treatment-P                  -0.306   1.33e- 7 4.91e- 7
##  8 conc    Ca_grain     Treatment-P                  19.0     1.32e- 6 4.68e- 6
##  9 conc    S_stover     Treatment-P                 113.      1.40e- 6 4.89e- 6
## 10 content Ca_stover_mg Treatment-P                 -97.4     1.49e- 5 4.94e- 5
## 11 hi      Ca_HI        Treatment-P                   0.00373 1.48e- 5 4.94e- 5
## 12 ratio   Ca_ratio     Treatment-P                   0.00416 1.83e- 5 5.95e- 5
## 13 content Ca_stover_mg GenotypeInv4m               -93.6     2.83e- 5 9.05e- 5
## 14 content Ca_grain_mg  Treatment-P                   1.26    6.12e- 5 1.90e- 4
## 15 content Mg_stover_mg Treatment-P                 -34.8     1.05e- 4 3.21e- 4
## 16 content K_stover_mg  Treatment-P                -299.      7.55e- 4 2.27e- 3
## 17 conc    S_stover     GenotypeInv4m:Treatment-P   -93.8     2.17e- 3 6.40e- 3
## 18 conc    Mg_grain     Treatment-P                 -97.5     2.46e- 3 7.15e- 3
## 19 hi      Fe_HI        GenotypeInv4m                 0.0403  3.33e- 3 9.56e- 3
## 20 conc    Mg_grain     GenotypeInv4m               -95.6     3.85e- 3 1.09e- 2
## 21 conc    Ca_stover    GenotypeInv4m              -411.      5.03e- 3 1.40e- 2
## 22 ratio   Mn_ratio     Treatment-P                  -0.0149  5.44e- 3 1.49e- 2
## 23 content P_grain_mg   Treatment-P                 -41.1     6.91e- 3 1.86e- 2
## 24 content Zn_stover_mg GenotypeInv4m                -0.440   7.06e- 3 1.86e- 2
## 25 ratio   Fe_ratio     GenotypeInv4m                 0.0531  7.02e- 3 1.86e- 2
## 26 content Fe_grain_mg  GenotypeInv4m                 0.206   7.85e- 3 2.04e- 2
## 27 ratio   Mg_ratio     Treatment-P                  -0.0673  8.95e- 3 2.29e- 2
## 28 content Mn_stover_mg Treatment-P                  -0.662   1.01e- 2 2.55e- 2
## 29 conc    S_grain      Treatment-P                  79.1     1.11e- 2 2.76e- 2
## 30 hi      Zn_HI        GenotypeInv4m                 0.0641  1.31e- 2 3.22e- 2
## 31 content K_stover_mg  GenotypeInv4m              -213.      1.39e- 2 3.39e- 2

emmeans contrast helpers

# Treatment-within-Genotype contrast (CTRL facet bracket = -P vs +P at CTRL;
# Inv4m facet bracket = -P vs +P at Inv4m). Returns one row per (response,
# Genotype).
contrasts_treatment_per_geno <- function(models) {
  lapply(names(models), function(resp) {
    m <- models[[resp]]
    emm <- emmeans(m, ~ Treatment | Genotype)
    pairs(emm, reverse = TRUE) %>%   # (-P) - (+P)
      as_tibble() %>%
      mutate(response = resp)
  }) %>% bind_rows()
}

# Pooled Genotype contrast (Inv4m vs CTRL averaged over Treatment). Returns
# one row per response.
contrasts_genotype_pooled <- function(models) {
  lapply(names(models), function(resp) {
    m <- models[[resp]]
    emm <- emmeans(m, ~ Genotype)
    pairs(emm, reverse = TRUE) %>%   # Inv4m - CTRL
      as_tibble() %>%
      mutate(response = resp)
  }) %>% bind_rows()
}

# Marginal Treatment contrast (pooled across Genotype). Returns one row per
# response. Use only when Genotype x Treatment interaction is ns.
contrasts_treatment_marginal <- function(models) {
  lapply(names(models), function(resp) {
    m <- models[[resp]]
    emm <- emmeans(m, ~ Treatment)
    pairs(emm, reverse = TRUE) %>%   # (-P) - (+P)
      as_tibble() %>%
      mutate(response = resp)
  }) %>% bind_rows()
}

# Format an emmeans contrast tibble to look like rstatix t_test output for
# stat_pvalue_manual. The function strips the contrast string and assigns
# group1/group2 explicitly. Adds FDR within the tibble (one metric).
format_treatment_contrast <- function(ct_tbl, response_to_mineral_pool, value_var) {
  ct_tbl %>%
    response_to_mineral_pool() %>%
    rename(p = p.value, statistic = t.ratio) %>%
    mutate(
      .y. = value_var,
      group1 = "+P",
      group2 = "-P",
      Genotype = gsub("Inv4m", "<i>Inv4m</i>", as.character(Genotype)),
      Genotype = factor(Genotype, levels = c("CTRL", "<i>Inv4m</i>"))
    ) %>%
    rstatix::adjust_pvalue(method = "fdr") %>%
    rstatix::add_significance("p.adj")
}

format_genotype_contrast <- function(ct_tbl, response_to_mineral_pool, value_var) {
  ct_tbl %>%
    response_to_mineral_pool() %>%
    rename(p = p.value, statistic = t.ratio) %>%
    mutate(
      .y. = value_var,
      group1 = "CTRL",
      group2 = "Inv4m"
    ) %>%
    rstatix::adjust_pvalue(method = "fdr") %>%
    rstatix::add_significance("p.adj")
}

Treatment contrasts → t-test-style tibbles for the main figure

# Helpers to split response name into mineral + pool
split_content <- function(df) {
  df %>% mutate(
    mineral = sub("_(grain|stover)_mg$", "", response),
    pool    = sub("^.+_(grain|stover)_mg$", "\\1", response),
    pool    = factor(pool, levels = c("stover", "grain"))
  )
}
split_conc <- function(df) {
  df %>% mutate(
    mineral = sub("_(grain|stover)$", "", response),
    pool    = sub("^.+_(grain|stover)$", "\\1", response),
    pool    = factor(pool, levels = c("stover", "grain"))
  )
}
split_hi <- function(df) {
  df %>% mutate(
    mineral = sub("_HI$", "", response)
  )
}
split_ratio <- function(df) {
  df %>% mutate(
    mineral = sub("_ratio$", "", response)
  )
}

# Cache raw emmeans outputs (before formatting) to avoid refitting on every
# render. Delete the RDS to force a refresh when models change.
emmeans_rds <- file.path(paths$intermediate, "ionome_content_emmeans.rds")
if (file.exists(emmeans_rds)) {
  cat("Loading cached emmeans contrasts from", emmeans_rds, "\n")
  .em <- readRDS(emmeans_rds)
  ct_conc    <- .em$ct_conc
  ct_content <- .em$ct_content
  ct_hi      <- .em$ct_hi
  ct_ratio   <- .em$ct_ratio
  rm(.em)
} else {
  ct_conc    <- contrasts_treatment_per_geno(models_conc)
  ct_content <- contrasts_treatment_per_geno(models_content)
  ct_hi      <- contrasts_treatment_per_geno(models_hi)
  ct_ratio   <- contrasts_treatment_per_geno(models_ratio)
  saveRDS(list(ct_conc = ct_conc, ct_content = ct_content,
               ct_hi = ct_hi, ct_ratio = ct_ratio), emmeans_rds)
  cat("Saved emmeans contrasts to", emmeans_rds, "\n")
}

## Loading cached emmeans contrasts from /Users/fvrodriguez/Library/CloudStorage/GoogleDrive-frodrig4@ncsu.edu/My Drive/repos/inv4m/results/phosphorus_paper/intermediate/ionome_content_emmeans.rds

# Unified FDR pool across every bracket drawn on the 4-panel figure:
#   - Panel A conc (per-Genotype)
#   - Panel B ratio (per-Genotype)
#   - Panel C content (per-Genotype)
#   - Panel D HI (per-Genotype)
raw_conc <- ct_conc %>%
  split_conc() %>%
  rename(p = p.value, statistic = t.ratio) %>%
  mutate(
    .y. = "ppm", group1 = "+P", group2 = "-P", metric = "conc",
    Genotype = gsub("Inv4m", "<i>Inv4m</i>", as.character(Genotype)),
    Genotype = factor(Genotype, levels = c("CTRL", "<i>Inv4m</i>"))
  )

raw_ratio <- ct_ratio %>%
  split_ratio() %>%
  rename(p = p.value, statistic = t.ratio) %>%
  mutate(
    .y. = "ratio", group1 = "+P", group2 = "-P", metric = "ratio",
    Genotype = gsub("Inv4m", "<i>Inv4m</i>", as.character(Genotype)),
    Genotype = factor(Genotype, levels = c("CTRL", "<i>Inv4m</i>"))
  )

raw_content <- ct_content %>%
  split_content() %>%
  rename(p = p.value, statistic = t.ratio) %>%
  mutate(
    .y. = "mg", group1 = "+P", group2 = "-P", metric = "content",
    Genotype = gsub("Inv4m", "<i>Inv4m</i>", as.character(Genotype)),
    Genotype = factor(Genotype, levels = c("CTRL", "<i>Inv4m</i>"))
  )

raw_hi <- ct_hi %>%
  split_hi() %>%
  rename(p = p.value, statistic = t.ratio) %>%
  mutate(
    .y. = "hi", group1 = "+P", group2 = "-P", metric = "hi",
    Genotype = gsub("Inv4m", "<i>Inv4m</i>", as.character(Genotype)),
    Genotype = factor(Genotype, levels = c("CTRL", "<i>Inv4m</i>"))
  )

unified_bracket_pool <- bind_rows(raw_conc, raw_ratio, raw_content,
                                  raw_hi) %>%
  mutate(p.adj = p.adjust(p, method = "fdr")) %>%
  rstatix::add_significance("p.adj")

cat("\n=== Unified FDR pool (n =", nrow(unified_bracket_pool),
    " | conc =", sum(unified_bracket_pool$metric == "conc"),
    "| ratio =", sum(unified_bracket_pool$metric == "ratio"),
    "| content =", sum(unified_bracket_pool$metric == "content"),
    "| hi =", sum(unified_bracket_pool$metric == "hi"),
    ") ===\n")

## 
## === Unified FDR pool (n = 96  | conc = 32 | ratio = 16 | content = 32 | hi = 16 ) ===

treatment.conc.ttest <- unified_bracket_pool %>%
  filter(metric == "conc") %>%
  mutate(Genotype = factor(Genotype, levels = c("CTRL", "<i>Inv4m</i>")),
         pool = factor(pool, levels = c("stover", "grain"))) %>%
  group_by(pool, mineral, Genotype) %>%
  group_modify(~ {
    sub <- to_plot_conc %>%
      filter(pool == .y$pool, mineral == .y$mineral, Genotype == .y$Genotype)
    .x %>%
      rstatix::add_y_position(
        data = sub, formula = ppm ~ Treatment,
        step.increase = 0.12
      )
  }) %>%
  ungroup() %>%
  filter(!is.na(y.position), is.finite(y.position))

treatment.ratio.ttest <- unified_bracket_pool %>%
  filter(metric == "ratio") %>%
  mutate(Genotype = factor(Genotype, levels = c("CTRL", "<i>Inv4m</i>"))) %>%
  group_by(mineral, Genotype) %>%
  group_modify(~ {
    sub <- to_plot_ratio %>%
      filter(mineral == .y$mineral, Genotype == .y$Genotype)
    .x %>%
      rstatix::add_y_position(
        data = sub, formula = ratio ~ Treatment,
        step.increase = 0.12
      )
  }) %>%
  ungroup() %>%
  filter(!is.na(y.position), is.finite(y.position))

treatment.t.test <- unified_bracket_pool %>%
  filter(metric == "content") %>%
  mutate(Genotype = factor(Genotype, levels = c("CTRL", "<i>Inv4m</i>")),
         pool = factor(pool, levels = c("stover", "grain"))) %>%
  group_by(pool, mineral, Genotype) %>%
  group_modify(~ {
    sub <- to_plot_content %>%
      filter(pool == .y$pool, mineral == .y$mineral, Genotype == .y$Genotype)
    .x %>%
      rstatix::add_y_position(
        data = sub, formula = mg ~ Treatment,
        step.increase = 0.20
      )
  }) %>%
  ungroup() %>%
  filter(!is.na(y.position), is.finite(y.position)) %>%
  mutate(y.position = case_when(
    mineral == "Ca" & pool == "grain" ~ y.position + 35,
    mineral == "Mn"                   ~ y.position + 0.2,
    mineral == "K"                    ~ y.position + 100,
    TRUE                              ~ y.position
  ))

treatment.hi.ttest <- unified_bracket_pool %>%
  filter(metric == "hi") %>%
  mutate(Genotype = factor(Genotype, levels = c("CTRL", "<i>Inv4m</i>"))) %>%
  group_by(mineral, Genotype) %>%
  group_modify(~ {
    sub <- to_plot_hi %>%
      filter(mineral == .y$mineral, Genotype == .y$Genotype)
    .x %>%
      rstatix::add_y_position(
        data = sub, formula = hi ~ Treatment,
        step.increase = 0.12
      )
  }) %>%
  ungroup() %>%
  filter(!is.na(y.position), is.finite(y.position))

cat("\n=== Treatment concentration contrasts (unified FDR, p.adj < 0.05) ===\n")

## 
## === Treatment concentration contrasts (unified FDR, p.adj < 0.05) ===

treatment.conc.ttest %>% filter(p.adj < 0.05) %>%
  select(mineral, pool, Genotype, statistic, p, p.adj, p.adj.signif) %>%
  arrange(p.adj) %>% print(n = Inf)

## # A tibble: 11 × 7
##    mineral pool   Genotype     statistic        p    p.adj p.adj.signif
##    <chr>   <fct>  <fct>            <dbl>    <dbl>    <dbl> <chr>       
##  1 P       stover <i>Inv4m</i>    -20.5  3.24e-28 3.11e-26 ****        
##  2 P       stover CTRL            -18.8  2.24e-26 1.07e-24 ****        
##  3 P       grain  <i>Inv4m</i>     -7.39 1.89e- 9 2.02e- 8 ****        
##  4 P       grain  CTRL             -7.17 4.02e- 9 3.85e- 8 ****        
##  5 Zn      stover CTRL              6.12 8.78e- 8 7.66e- 7 ****        
##  6 S       stover CTRL              5.38 1.40e- 6 9.61e- 6 ****        
##  7 Ca      grain  CTRL              5.53 1.32e- 6 9.61e- 6 ****        
##  8 Ca      grain  <i>Inv4m</i>      4.92 1.07e- 5 6.42e- 5 ****        
##  9 Zn      stover <i>Inv4m</i>      4.63 2.08e- 5 9.97e- 5 ****        
## 10 Mg      grain  CTRL             -3.20 2.46e- 3 7.86e- 3 **          
## 11 S       grain  CTRL              2.64 1.11e- 2 2.80e- 2 *

cat("\n=== Treatment ratio contrasts (unified FDR, p.adj < 0.05) ===\n")

## 
## === Treatment ratio contrasts (unified FDR, p.adj < 0.05) ===

treatment.ratio.ttest %>% filter(p.adj < 0.05) %>%
  select(mineral, Genotype, statistic, p, p.adj, p.adj.signif) %>%
  arrange(p.adj) %>% print(n = Inf)

## # A tibble: 9 × 6
##   mineral Genotype     statistic        p    p.adj p.adj.signif
##   <chr>   <fct>            <dbl>    <dbl>    <dbl> <chr>       
## 1 P       CTRL             12.9  3.45e-17 1.06e-15 ****        
## 2 P       <i>Inv4m</i>     12.5  1.13e-16 2.17e-15 ****        
## 3 Zn      CTRL             -6.18 1.33e- 7 1.06e- 6 ****        
## 4 Ca      CTRL              4.76 1.83e- 5 9.23e- 5 ****        
## 5 Ca      <i>Inv4m</i>      4.24 9.98e- 5 4.17e- 4 ***         
## 6 Zn      <i>Inv4m</i>     -3.63 6.80e- 4 2.51e- 3 **          
## 7 Mn      CTRL             -2.91 5.44e- 3 1.68e- 2 *           
## 8 Fe      <i>Inv4m</i>     -2.81 7.16e- 3 2.02e- 2 *           
## 9 Mg      CTRL             -2.72 8.95e- 3 2.46e- 2 *

cat("\n=== Treatment content contrasts (unified FDR, p.adj < 0.05) ===\n")

## 
## === Treatment content contrasts (unified FDR, p.adj < 0.05) ===

treatment.t.test %>% filter(p.adj < 0.05) %>%
  select(mineral, pool, Genotype, statistic, p, p.adj, p.adj.signif) %>%
  arrange(p.adj) %>% print(n = Inf)

## # A tibble: 14 × 7
##    mineral pool   Genotype     statistic        p    p.adj p.adj.signif
##    <chr>   <fct>  <fct>            <dbl>    <dbl>    <dbl> <chr>       
##  1 P       stover CTRL            -11.8  4.40e-17 1.06e-15 ****        
##  2 P       stover <i>Inv4m</i>    -11.1  4.74e-16 7.59e-15 ****        
##  3 Ca      stover CTRL             -4.73 1.49e- 5 7.96e- 5 ****        
##  4 P       grain  <i>Inv4m</i>     -4.69 2.33e- 5 1.06e- 4 ***         
##  5 Ca      grain  CTRL              4.39 6.12e- 5 2.67e- 4 ***         
##  6 Mg      stover CTRL             -4.16 1.05e- 4 4.21e- 4 ***         
##  7 Ca      grain  <i>Inv4m</i>      4.18 1.24e- 4 4.77e- 4 ***         
##  8 K       stover CTRL             -3.56 7.55e- 4 2.68e- 3 **          
##  9 Mg      stover <i>Inv4m</i>     -3.40 1.21e- 3 4.16e- 3 **          
## 10 Ca      stover <i>Inv4m</i>     -3.37 1.34e- 3 4.43e- 3 **          
## 11 Mn      stover <i>Inv4m</i>     -2.82 6.59e- 3 1.98e- 2 *           
## 12 P       grain  CTRL             -2.82 6.91e- 3 2.01e- 2 *           
## 13 Mn      stover CTRL             -2.66 1.01e- 2 2.62e- 2 *           
## 14 Mn      grain  <i>Inv4m</i>     -2.69 9.90e- 3 2.62e- 2 *

cat("\n=== Treatment HI contrasts (unified FDR, p.adj < 0.05) ===\n")

## 
## === Treatment HI contrasts (unified FDR, p.adj < 0.05) ===

treatment.hi.ttest %>% filter(p.adj < 0.05) %>%
  select(mineral, Genotype, statistic, p, p.adj, p.adj.signif) %>%
  arrange(p.adj) %>% print(n = Inf)

## # A tibble: 4 × 6
##   mineral Genotype     statistic        p    p.adj p.adj.signif
##   <chr>   <fct>            <dbl>    <dbl>    <dbl> <chr>       
## 1 P       CTRL              9.16 4.10e-12 5.63e-11 ****        
## 2 P       <i>Inv4m</i>      7.89 3.24e-10 3.89e- 9 ****        
## 3 Ca      <i>Inv4m</i>      5.37 2.27e- 6 1.45e- 5 ****        
## 4 Ca      CTRL              4.82 1.48e- 5 7.96e- 5 ****

Genotype contrasts (pooled) → unified-FDR pool across conc / content / HI / ratio

# Cache raw pooled-Genotype emmeans outputs (including ratio).
geno_emmeans_rds <- file.path(paths$intermediate,
                              "ionome_content_geno_emmeans.rds")
needs_recompute <- TRUE
if (file.exists(geno_emmeans_rds)) {
  .gm <- readRDS(geno_emmeans_rds)
  if (all(c("gp_conc", "gp_content", "gp_hi", "gp_ratio") %in% names(.gm))) {
    gp_conc    <- .gm$gp_conc
    gp_content <- .gm$gp_content
    gp_hi      <- .gm$gp_hi
    gp_ratio   <- .gm$gp_ratio
    needs_recompute <- FALSE
    cat("Loaded cached pooled-Genotype contrasts from", geno_emmeans_rds, "\n")
  }
  rm(.gm)
}

## Loaded cached pooled-Genotype contrasts from /Users/fvrodriguez/Library/CloudStorage/GoogleDrive-frodrig4@ncsu.edu/My Drive/repos/inv4m/results/phosphorus_paper/intermediate/ionome_content_geno_emmeans.rds

if (needs_recompute) {
  gp_conc    <- contrasts_genotype_pooled(models_conc)
  gp_content <- contrasts_genotype_pooled(models_content)
  gp_hi      <- contrasts_genotype_pooled(models_hi)
  gp_ratio   <- contrasts_genotype_pooled(models_ratio)
  saveRDS(list(gp_conc = gp_conc, gp_content = gp_content,
               gp_hi = gp_hi, gp_ratio = gp_ratio), geno_emmeans_rds)
  cat("Saved pooled-Genotype contrasts to", geno_emmeans_rds, "\n")
}

# Unified pooled-Genotype FDR across all four metrics. Mirrors the
# unified-bracket convention used in the main 4-panel figure.
format_raw_geno <- function(ct_tbl, response_split, value_var, metric_lbl) {
  ct_tbl %>%
    response_split() %>%
    rename(p = p.value, statistic = t.ratio) %>%
    mutate(.y. = value_var, group1 = "CTRL", group2 = "Inv4m",
           metric = metric_lbl)
}

unified_geno_pool <- bind_rows(
  format_raw_geno(gp_conc,    split_conc,    "ppm",   "conc"),
  format_raw_geno(gp_content, split_content, "mg",    "content"),
  format_raw_geno(gp_hi,      split_hi,      "hi",    "hi"),
  format_raw_geno(gp_ratio,   split_ratio,   "ratio", "ratio")
) %>%
  rstatix::adjust_pvalue(method = "fdr", p.col = "p", output.col = "p.adj") %>%
  rstatix::add_significance("p.adj")

cat(sprintf("Unified pooled-Genotype FDR: %d contrasts, %d sig at p.adj < 0.05\n",
            nrow(unified_geno_pool), sum(unified_geno_pool$p.adj < 0.05)))

## Unified pooled-Genotype FDR: 48 contrasts, 14 sig at p.adj < 0.05

# Clean copies of the long-format data with non-markdown Genotype labels
# (needed for rstatix::add_y_position which matches on Genotype values).
to_plot_conc_clean <- to_plot_conc %>%
  mutate(Genotype = gsub("<i>|</i>", "", as.character(Genotype)),
         Genotype = factor(Genotype, levels = c("CTRL", "Inv4m")))
to_plot_content_clean <- to_plot_content %>%
  mutate(Genotype = gsub("<i>|</i>", "", as.character(Genotype)),
         Genotype = factor(Genotype, levels = c("CTRL", "Inv4m")))
to_plot_hi_clean <- to_plot_hi %>%
  mutate(Genotype = gsub("<i>|</i>", "", as.character(Genotype)),
         Genotype = factor(Genotype, levels = c("CTRL", "Inv4m")))
to_plot_ratio_clean <- to_plot_ratio %>%
  mutate(Genotype = gsub("<i>|</i>", "", as.character(Genotype)),
         Genotype = factor(Genotype, levels = c("CTRL", "Inv4m")))

add_y_pos_pool <- function(tbl, data, formula_y) {
  tbl %>%
    group_by(pool, mineral) %>%
    group_modify(~ {
      sub <- data %>% filter(pool == .y$pool, mineral == .y$mineral)
      .x %>% rstatix::add_y_position(
        data = sub, formula = formula_y, step.increase = 0.12
      )
    }) %>%
    ungroup() %>%
    filter(!is.na(y.position), is.finite(y.position))
}

add_y_pos_nopool <- function(tbl, data, formula_y) {
  tbl %>%
    group_by(mineral) %>%
    group_modify(~ {
      sub <- data %>% filter(mineral == .y$mineral)
      .x %>% rstatix::add_y_position(
        data = sub, formula = formula_y, step.increase = 0.12
      )
    }) %>%
    ungroup() %>%
    filter(!is.na(y.position), is.finite(y.position))
}

genotype.conc.ttest <- unified_geno_pool %>%
  filter(metric == "conc") %>%
  add_y_pos_pool(to_plot_conc_clean, ppm ~ Genotype)

genotype.content.ttest <- unified_geno_pool %>%
  filter(metric == "content") %>%
  add_y_pos_pool(to_plot_content_clean, mg ~ Genotype)

genotype.hi.ttest <- unified_geno_pool %>%
  filter(metric == "hi") %>%
  add_y_pos_nopool(to_plot_hi_clean, hi ~ Genotype)

genotype.ratio.ttest <- unified_geno_pool %>%
  filter(metric == "ratio") %>%
  add_y_pos_nopool(to_plot_ratio_clean, ratio ~ Genotype)

cat("\n=== Inv4m contrasts (unified pooled-Genotype FDR, p.adj < 0.05) ===\n")

## 
## === Inv4m contrasts (unified pooled-Genotype FDR, p.adj < 0.05) ===

unified_geno_pool %>%
  filter(p.adj < 0.05) %>%
  select(metric, mineral, any_of("pool"), statistic, p, p.adj, p.adj.signif) %>%
  arrange(p.adj) %>%
  print(n = Inf)

## # A tibble: 14 × 7
##    metric  mineral pool   statistic           p     p.adj p.adj.signif
##    <chr>   <chr>   <fct>      <dbl>       <dbl>     <dbl> <chr>       
##  1 content Ca      stover     -5.48 0.000000959 0.0000460 ****        
##  2 conc    Ca      stover     -4.97 0.00000620  0.000149  ***         
##  3 conc    Mg      grain      -3.85 0.000344    0.00537   **          
##  4 content Zn      stover     -3.72 0.000448    0.00537   **          
##  5 conc    K       grain      -3.43 0.00125     0.0108    *           
##  6 hi      Fe      <NA>        3.41 0.00134     0.0108    *           
##  7 hi      Zn      <NA>        3.28 0.00195     0.0134    *           
##  8 conc    Zn      stover     -3.19 0.00229     0.0137    *           
##  9 conc    Mn      grain      -3.12 0.00305     0.0162    *           
## 10 conc    P       grain      -2.96 0.00478     0.0230    *           
## 11 ratio   Fe      <NA>        2.89 0.00575     0.0251    *           
## 12 content S       stover     -2.71 0.00889     0.0356    *           
## 13 content K       stover     -2.65 0.0104      0.0383    *           
## 14 content Fe      grain       2.63 0.0115      0.0394    *

Per-mineral plot functions (unchanged from base notebook)

pal <- c("+P" = "gold", "-P" = "purple4")

plot_mineral_content <- function(mineral_name, base_size = 18) {
  plot_data <- to_plot_content %>% filter(mineral == mineral_name)
  ttest_subset <- treatment.t.test %>% filter(mineral == mineral_name)

  ggplot(plot_data, aes(x = Treatment, y = mg, color = Treatment)) +
    geom_boxplot(width = 0.5, linewidth = 1, outlier.shape = NA) %>%
    with_shadow(colour = "black", x_offset = 0, y_offset = 0, sigma = 1) +
    geom_quasirandom(size = 1) %>%
    with_shadow(colour = "black", x_offset = 0, y_offset = 0, sigma = 1) +
    stat_pvalue_manual(ttest_subset, label = "p.adj.signif",
                       size = 8, bracket.size = 1, hide.ns = TRUE) +
    scale_color_manual(values = pal) +
    scale_y_continuous(expand = expansion(mult = c(0.05, 0.25))) +
    ggh4x::facet_nested(~ pool + Genotype, scales = "free_y") +
    labs(title = mineral_name, x = NULL, y = "Content [mg / plant]") +
    theme_classic(base_size = base_size) +
    theme(
      legend.position = "none",
      axis.text.x = element_text(color = "black", face = "bold", size = rel(1.0)),
      axis.title.y = element_text(face = "bold", size = rel(0.9)),
      plot.title = element_text(face = "bold", size = rel(1.55), hjust = 0.5),
      strip.text = element_markdown(face = "bold", size = rel(0.9)),
      strip.background = element_blank(),
      panel.spacing = unit(0.1, "lines"),
      ggh4x.facet.nestline = element_line(color = "black", linewidth = 0.5)
    )
}

plot_mineral_conc <- function(mineral_name, log_scale = FALSE, base_size = 18) {
  plot_data <- to_plot_conc %>% filter(mineral == mineral_name)
  ttest_subset <- treatment.conc.ttest %>% filter(mineral == mineral_name)

  p <- ggplot(plot_data, aes(x = Treatment, y = ppm, color = Treatment)) +
    geom_boxplot(width = 0.5, linewidth = 1, outlier.shape = NA) %>%
    with_shadow(colour = "black", x_offset = 0, y_offset = 0, sigma = 1) +
    geom_quasirandom(size = 1) %>%
    with_shadow(colour = "black", x_offset = 0, y_offset = 0, sigma = 1) +
    scale_color_manual(values = pal) +
    labs(title = mineral_name, x = NULL, y = "Concentration [mg / kg]")

  if (log_scale) {
    ttest_log <- ttest_subset
    ttest_log$y.position <- log2(ttest_log$y.position)
    p <- p +
      scale_y_continuous(trans = "log2",
        breaks = c(5, 10, 25, 50, 100, 200, 400, 800, 1600, 3200, 6400),
        expand = expansion(mult = c(0.05, 0.15))) +
      stat_pvalue_manual(ttest_log, label = "p.adj.signif",
                         size = 8, bracket.size = 1, hide.ns = TRUE)
  } else {
    p <- p + scale_y_continuous(expand = expansion(mult = c(0.05, 0.15))) +
      stat_pvalue_manual(ttest_subset, label = "p.adj.signif",
                         size = 8, bracket.size = 1, hide.ns = TRUE)
  }

  p + ggh4x::facet_nested(~ pool + Genotype, scales = "free_y") +
    theme_classic(base_size = base_size) +
    theme(
      legend.position = "none",
      axis.text.x = element_text(color = "black", face = "bold", size = rel(1.0)),
      axis.title.y = element_text(face = "bold", size = rel(0.9)),
      plot.title = element_text(face = "bold", size = rel(1.55), hjust = 0.5),
      strip.text = element_markdown(face = "bold", size = rel(0.9)),
      strip.background = element_blank(),
      panel.spacing = unit(0.1, "lines"),
      ggh4x.facet.nestline = element_line(color = "black", linewidth = 0.5)
    )
}

plot_mineral_ratio <- function(mineral_name, base_size = 18) {
  plot_data <- to_plot_ratio %>% filter(mineral == mineral_name)
  ttest_subset <- treatment.ratio.ttest %>% filter(mineral == mineral_name)

  ggplot(plot_data, aes(x = Treatment, y = ratio, color = Treatment)) +
    geom_boxplot(width = 0.5, linewidth = 1, outlier.shape = NA) %>%
    with_shadow(colour = "black", x_offset = 0, y_offset = 0, sigma = 1) +
    geom_quasirandom(size = 1) %>%
    with_shadow(colour = "black", x_offset = 0, y_offset = 0, sigma = 1) +
    stat_pvalue_manual(ttest_subset, label = "p.adj.signif",
                       size = 8, bracket.size = 1, hide.ns = TRUE) +
    scale_color_manual(values = pal) +
    scale_y_continuous(expand = expansion(mult = c(0.05, 0.15))) +
    facet_wrap(~ Genotype, nrow = 1) +
    labs(title = mineral_name, x = NULL, y = "Grain / stover ratio") +
    theme_classic(base_size = base_size) +
    theme(
      legend.position = "none",
      axis.text.x = element_text(color = "black", face = "bold", size = rel(1.0)),
      axis.title.y = element_text(face = "bold", size = rel(0.9)),
      plot.title = element_text(face = "bold", size = rel(1.55), hjust = 0.5),
      strip.text = element_markdown(face = "bold", size = rel(0.9)),
      strip.background = element_blank(),
      panel.spacing = unit(0.1, "lines")
    )
}

plot_mineral_HI_marginal <- function(mineral_name, base_size = 18) {
  plot_data <- to_plot_hi %>% filter(mineral == mineral_name)
  ttest_subset <- treatment.hi.marginal.ttest %>% filter(mineral == mineral_name)

  p <- ggplot(plot_data, aes(x = Treatment, y = hi, color = Treatment)) +
    geom_boxplot(width = 0.5, linewidth = 1, outlier.shape = NA) %>%
    with_shadow(colour = "black", x_offset = 0, y_offset = 0, sigma = 1) +
    geom_quasirandom(size = 1) %>%
    with_shadow(colour = "black", x_offset = 0, y_offset = 0, sigma = 1) +
    scale_color_manual(values = pal) +
    scale_y_continuous(expand = expansion(mult = c(0.05, 0.15))) +
    labs(title = mineral_name, x = NULL, y = "Mineral HI [grain / total]") +
    theme_classic(base_size = base_size) +
    theme(
      legend.position = "none",
      axis.text.x = element_text(color = "black", face = "bold", size = rel(1.0)),
      axis.title.y = element_text(face = "bold", size = rel(0.9)),
      plot.title = element_text(face = "bold", size = rel(1.55), hjust = 0.5),
      strip.text = element_markdown(face = "bold", size = rel(0.9)),
      strip.background = element_blank()
    )

  if (nrow(ttest_subset) > 0) {
    p <- p + stat_pvalue_manual(ttest_subset, label = "p.adj.signif",
                                size = 8, bracket.size = 1, hide.ns = TRUE)
  }
  p
}

plot_mineral_HI <- function(mineral_name, base_size = 18) {
  plot_data <- to_plot_hi %>% filter(mineral == mineral_name)
  ttest_subset <- treatment.hi.ttest %>% filter(mineral == mineral_name)

  p <- ggplot(plot_data, aes(x = Treatment, y = hi, color = Treatment)) +
    geom_boxplot(width = 0.5, linewidth = 1, outlier.shape = NA) %>%
    with_shadow(colour = "black", x_offset = 0, y_offset = 0, sigma = 1) +
    geom_quasirandom(size = 1) %>%
    with_shadow(colour = "black", x_offset = 0, y_offset = 0, sigma = 1) +
    scale_color_manual(values = pal) +
    scale_y_continuous(expand = expansion(mult = c(0.05, 0.15))) +
    labs(title = mineral_name, x = NULL, y = "Mineral HI [grain / total]") +
    facet_wrap(~ Genotype) +
    theme_classic(base_size = base_size) +
    theme(
      legend.position = "none",
      axis.text.x = element_text(color = "black", face = "bold", size = rel(1.0)),
      axis.title.y = element_text(face = "bold", size = rel(0.9)),
      plot.title = element_text(face = "bold", size = rel(1.55), hjust = 0.5),
      strip.text = element_markdown(face = "bold", size = rel(0.9)),
      strip.background = element_blank()
    )

  if (nrow(ttest_subset) > 0) {
    p <- p + stat_pvalue_manual(ttest_subset, label = "p.adj.signif",
                                size = 8, bracket.size = 1, hide.ns = TRUE)
  }
  p
}

Main 4-panel figure

rowA_main <- ggarrange(
  plot_mineral_conc("P"),
  plot_mineral_conc("Zn") + ggpubr::rremove("ylab"),
  plot_mineral_conc("Ca", log_scale = TRUE) + ggpubr::rremove("ylab"),
  plot_mineral_conc("S")  + ggpubr::rremove("ylab"),
  ncol = 4, nrow = 1, align = "h"
)

rowB_main <- ggarrange(
  plot_mineral_ratio("P") + labs(title = NULL) +
    scale_y_continuous(breaks = 1:5,
                       expand = expansion(mult = c(0.05, 0.15))),
  plot_mineral_ratio("Zn") + labs(title = NULL) + ggpubr::rremove("ylab"),
  plot_mineral_ratio("Ca") + labs(title = NULL) + ggpubr::rremove("ylab"),
  plot_mineral_ratio("S")  + labs(title = NULL) + ggpubr::rremove("ylab"),
  ncol = 4, nrow = 1, align = "h"
)

rowC_main <- ggarrange(
  plot_mineral_content("P")  + labs(title = NULL),
  plot_mineral_content("Zn") + labs(title = NULL) + ggpubr::rremove("ylab"),
  plot_mineral_content("Ca") + labs(title = NULL) + ggpubr::rremove("ylab"),
  plot_mineral_content("S")  + labs(title = NULL) + ggpubr::rremove("ylab"),
  ncol = 4, nrow = 1, align = "h"
)

rowD_main <- ggarrange(
  plot_mineral_HI("P")  + labs(title = NULL),
  plot_mineral_HI("Zn") + labs(title = NULL) + ggpubr::rremove("ylab"),
  plot_mineral_HI("Ca") + labs(title = NULL) + ggpubr::rremove("ylab"),
  plot_mineral_HI("S")  + labs(title = NULL) + ggpubr::rremove("ylab"),
  ncol = 4, nrow = 1, align = "h"
)

p_main_4panel <- ggarrange(
  rowA_main, rowB_main, rowC_main, rowD_main,
  ncol = 1, nrow = 4,
  heights = c(1, 0.8, 0.9, 0.8),
  labels = c("A", "B", "C", "D"),
  font.label = list(size = 30),
  align = "hv"
)
print(p_main_4panel)

Supplementary 4-panel figure

rowA_supp <- ggarrange(
  plot_mineral_conc("Mg"),
  plot_mineral_conc("Mn", log_scale = TRUE) + ggpubr::rremove("ylab"),
  plot_mineral_conc("K")  + ggpubr::rremove("ylab"),
  plot_mineral_conc("Fe") + ggpubr::rremove("ylab"),
  ncol = 4, nrow = 1, align = "h"
)

rowB_supp <- ggarrange(
  plot_mineral_ratio("Mg") + labs(title = NULL),
  plot_mineral_ratio("Mn") + labs(title = NULL) + ggpubr::rremove("ylab"),
  plot_mineral_ratio("K")  + labs(title = NULL) + ggpubr::rremove("ylab"),
  plot_mineral_ratio("Fe") + labs(title = NULL) + ggpubr::rremove("ylab"),
  ncol = 4, nrow = 1, align = "h"
)

rowC_supp <- ggarrange(
  plot_mineral_content("Mg") + labs(title = NULL),
  plot_mineral_content("Mn") + labs(title = NULL) + ggpubr::rremove("ylab"),
  plot_mineral_content("K")  + labs(title = NULL) + ggpubr::rremove("ylab"),
  plot_mineral_content("Fe") + labs(title = NULL) + ggpubr::rremove("ylab"),
  ncol = 4, nrow = 1, align = "h"
)

rowD_supp <- ggarrange(
  plot_mineral_HI("Mg") + labs(title = NULL),
  plot_mineral_HI("Mn") + labs(title = NULL) + ggpubr::rremove("ylab"),
  plot_mineral_HI("K")  + labs(title = NULL) + ggpubr::rremove("ylab"),
  plot_mineral_HI("Fe") + labs(title = NULL) + ggpubr::rremove("ylab"),
  ncol = 4, nrow = 1, align = "h"
)

p_supp_4panel <- ggarrange(
  rowA_supp, rowB_supp, rowC_supp, rowD_supp,
  ncol = 1, nrow = 4,
  heights = c(1, 0.8, 0.9, 0.8),
  labels = c("A", "B", "C", "D"),
  font.label = list(size = 30),
  align = "hv"
)
print(p_supp_4panel)

Inv4m supp figure (4 metrics x 5 ions; Figure-2-style organization)

Rows are metrics (conc, content, HI, ratio) and columns are ions (Mg, Ca, Zn, Fe, K). For conc and content, a single pool is shown per cell — the pool where the Inv4m effect is significant under the unified FDR, or stover as a fallback. Ion symbols are shown as plot titles on the top row (A, conc); pool (grain / stover) appears as a facet strip on A and B. Brackets are drawn only when the unified pooled-Genotype FDR is < 0.05.

pal_geno <- c("CTRL" = "gold", "Inv4m" = "purple4")

plot_inv4m_cell <- function(mineral_name, metric, pool_name = NULL,
                            show_title = FALSE, show_ylab = TRUE,
                            base_size = 18, log_scale = FALSE) {
  if (metric == "conc") {
    d <- to_plot_conc %>% filter(mineral == mineral_name, pool == pool_name)
    s <- genotype.conc.ttest %>% filter(mineral == mineral_name, pool == pool_name)
    y_var <- "ppm"
    ylab_text <- "Concentration [mg / kg]"
  } else if (metric == "content") {
    d <- to_plot_content %>% filter(mineral == mineral_name, pool == pool_name)
    s <- genotype.content.ttest %>% filter(mineral == mineral_name, pool == pool_name)
    y_var <- "mg"
    ylab_text <- "Content [mg / plant]"
  } else if (metric == "hi") {
    d <- to_plot_hi %>% filter(mineral == mineral_name)
    s <- genotype.hi.ttest %>% filter(mineral == mineral_name)
    y_var <- "hi"
    ylab_text <- "Mineral HI [grain / total]"
  } else if (metric == "ratio") {
    d <- to_plot_ratio %>% filter(mineral == mineral_name)
    s <- genotype.ratio.ttest %>% filter(mineral == mineral_name)
    y_var <- "ratio"
    ylab_text <- "Grain / stover ratio"
  }

  d <- d %>%
    mutate(Genotype = gsub("<i>|</i>", "", as.character(Genotype)),
           Genotype = factor(Genotype, levels = c("CTRL", "Inv4m")))

  p <- ggplot(d, aes(x = Genotype, y = .data[[y_var]], colour = Genotype)) +
    geom_boxplot(width = 0.5, linewidth = 1, outlier.shape = NA) %>%
      with_shadow(colour = "black", x_offset = 0, y_offset = 0, sigma = 1) +
    geom_quasirandom(size = 1.2) %>%
      with_shadow(colour = "black", x_offset = 0, y_offset = 0, sigma = 1) +
    scale_colour_manual(values = pal_geno) +
    scale_x_discrete(labels = c("CTRL" = "CTRL", "Inv4m" = "<i>Inv4m</i>")) +
    labs(
      title = if (show_title) mineral_name
              else if (metric %in% c("hi", "ratio")) " "
              else NULL,
      x = NULL,
      y = if (show_ylab) ylab_text else NULL
    )

  if (log_scale) {
    p <- p + scale_y_continuous(
      trans = "log10", expand = expansion(mult = c(0.05, 0.15))
    )
  } else {
    p <- p + scale_y_continuous(expand = expansion(mult = c(0.05, 0.15)))
  }

  if (metric %in% c("conc", "content")) {
    p <- p + facet_wrap(~ pool, nrow = 1, strip.position = "top")
  }

  if (nrow(s) > 0) {
    p <- p + stat_pvalue_manual(
      s, label = "p.adj.signif",
      size = 8, bracket.size = 1, hide.ns = TRUE
    )
  }

  p + theme_classic(base_size = base_size) +
    theme(
      legend.position = "none",
      axis.text.x = element_markdown(face = "bold", size = rel(1.0), color = "black"),
      axis.title.y = element_text(face = "bold", size = rel(0.9)),
      plot.title = element_text(face = "bold", size = rel(1.55), hjust = 0.5),
      strip.text = element_text(face = "bold", size = rel(0.9)),
      strip.background = element_blank()
    )
}

ions <- c("Mg", "Ca", "Zn", "Fe", "K")

# Selected pool per cell for conc / content. Sig pool if any, else stover.
pool_conc    <- c(Mg = "grain",  Ca = "stover", Zn = "stover",
                  Fe = "stover", K  = "stover")
pool_content <- c(Mg = "stover", Ca = "stover", Zn = "stover",
                  Fe = "grain",  K  = "stover")

# Row order: A = conc (with ion titles), B = ratio, C = content, D = HI
row_A <- lapply(ions, function(m) {
  plot_inv4m_cell(m, "conc", pool_name = pool_conc[[m]],
                  show_title = TRUE, show_ylab = (m == "Mg"))
})
row_B <- lapply(ions, function(m) {
  plot_inv4m_cell(m, "ratio",
                  show_title = FALSE, show_ylab = (m == "Mg"))
})
row_C <- lapply(ions, function(m) {
  plot_inv4m_cell(m, "content", pool_name = pool_content[[m]],
                  show_title = FALSE, show_ylab = (m == "Mg"))
})
row_D <- lapply(ions, function(m) {
  plot_inv4m_cell(m, "hi",
                  show_title = FALSE, show_ylab = (m == "Mg"))
})

p_inv4m_effects <- ggarrange(
  ggarrange(plotlist = row_A, ncol = 5, align = "h"),
  ggarrange(plotlist = row_B, ncol = 5, align = "h"),
  ggarrange(plotlist = row_C, ncol = 5, align = "h"),
  ggarrange(plotlist = row_D, ncol = 5, align = "h"),
  ncol = 1, nrow = 4,
  heights = c(1, 1, 1, 1),
  labels = c("A", "B", "C", "D"),
  font.label = list(size = 30)
)

print(p_inv4m_effects)

Save outputs

write.csv(
  effects,
  file.path(paths$intermediate, "ionome_content_linear_model_effects.csv"),
  row.names = FALSE
)

write_ttest <- function(x, name) {
  if (nrow(x) == 0) return(invisible(NULL))
  write.csv(
    x %>% select(where(~ !is.list(.))),
    file.path(paths$intermediate, name),
    row.names = FALSE
  )
}

write_ttest(treatment.t.test,       "ionome_content_treatment_ttests.csv")
write_ttest(treatment.conc.ttest,   "ionome_content_conc_treatment_ttests.csv")
write_ttest(treatment.hi.ttest,     "ionome_content_hi_treatment_ttests.csv")
write_ttest(treatment.ratio.ttest,  "ionome_content_ratio_treatment_ttests.csv")
write_ttest(genotype.conc.ttest,    "ionome_content_conc_genotype_ttests.csv")
write_ttest(genotype.content.ttest, "ionome_content_genotype_ttests.csv")
write_ttest(genotype.hi.ttest,      "ionome_content_hi_genotype_ttests.csv")
write_ttest(genotype.ratio.ttest,   "ionome_content_ratio_genotype_ttests.csv")
write_ttest(unified_geno_pool %>% select(where(~ !is.list(.))),
            "ionome_content_unified_genotype_ttests.csv")

ggsave(file.path(paths$figures, "ionome_4panel_main_figure.png"),
       plot = p_main_4panel, width = 15, height = 16, units = "in", dpi = 300)
ggsave(file.path(paths$figures, "ionome_4panel_supplementary_figure.png"),
       plot = p_supp_4panel, width = 15, height = 16, units = "in", dpi = 300)
ggsave(file.path(paths$figures, "ionome_inv4m_effects_figure.png"),
       plot = p_inv4m_effects, width = 14, height = 16, units = "in", dpi = 300)

Session info

sessionInfo()

## R version 4.5.2 (2025-10-31)
## Platform: aarch64-apple-darwin20
## Running under: macOS Tahoe 26.3.1
## 
## Matrix products: default
## BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
## 
## time zone: America/New_York
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] emmeans_2.0.1    ggh4x_0.3.1      ggtext_0.1.2     rstatix_0.7.3   
##  [5] ggfx_1.0.3       ggbeeswarm_0.7.3 ggpubr_0.6.2     tibble_3.3.0    
##  [9] ggplot2_4.0.1    tidyr_1.3.2      dplyr_1.1.4     
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6       beeswarm_0.4.0     xfun_0.55          bslib_0.9.0       
##  [5] lattice_0.22-7     vctrs_0.6.5        tools_4.5.2        generics_0.1.4    
##  [9] pkgconfig_2.0.3    RColorBrewer_1.1-3 S7_0.2.1           lifecycle_1.0.4   
## [13] stringr_1.6.0      compiler_4.5.2     farver_2.1.2       textshaping_1.0.4 
## [17] carData_3.0-5      litedown_0.9       vipor_0.4.7        htmltools_0.5.9   
## [21] sass_0.4.10        yaml_2.3.12        Formula_1.2-5      pillar_1.11.1     
## [25] car_3.1-3          jquerylib_0.1.4    cachem_1.1.0       magick_2.9.0      
## [29] abind_1.4-8        commonmark_2.0.0   tidyselect_1.2.1   digest_0.6.39     
## [33] stringi_1.8.7      mvtnorm_1.3-3      purrr_1.2.0        labeling_0.4.3    
## [37] cowplot_1.2.0      rprojroot_2.1.1    fastmap_1.2.0      grid_4.5.2        
## [41] here_1.0.2         cli_3.6.5          magrittr_2.0.4     utf8_1.2.6        
## [45] dichromat_2.0-0.1  broom_1.0.11       withr_3.0.2        scales_1.4.0      
## [49] backports_1.5.0    estimability_1.5.1 rmarkdown_2.30     otel_0.2.0        
## [53] ggsignif_0.6.4     ragg_1.5.0         coda_0.19-4.1      evaluate_1.0.5    
## [57] knitr_1.51         markdown_2.0       rlang_1.1.6        gridtext_0.1.5    
## [61] Rcpp_1.1.0         xtable_1.8-4       glue_1.8.0         xml2_1.5.1        
## [65] jsonlite_2.0.0     R6_2.6.1           systemfonts_1.3.1