PSU 2022 Inv4m Phenotype Analysis: Marginal Effects Approach

Overview

This analysis uses a two-stage approach: 1. Test all model coefficients (Genotype, Treatment, Interaction) with single FDR correction 2. For traits without significant interactions: report marginal contrasts (averaged effects) For traits with significant interactions: report conditional contrasts (within-level effects)

---

1. Setup and Data Loading

Load libraries

library(tidyverse)
library(dplyr)
library(ggplot2)
library(ggtext)
library(ggfx)
library(ggpubr)
library(ggbeeswarm)
library(rstatix)
library(emmeans)
library(knitr)

Define aesthetic constants

pal <- c("gold", "purple4")

pheno_theme <- theme_classic2(base_size = 20) +
  theme(
    plot.title = element_markdown(hjust = 0.5, face = "bold"),
    axis.title.x = element_blank(),
    axis.text.x = element_text(
      face = "bold", 
      color = "black", 
      size = 20, 
      angle = 45, 
      hjust = 1
    ),
    strip.background = element_blank(),
    strip.text = element_text(size = 25),
    legend.position = "none"
  )

pheno_theme2 <- theme_classic2(base_size = 20) +
  theme(
    plot.title = element_markdown(hjust = 0.5, face = "bold"),
    axis.title.y = element_markdown(),
    axis.title.x = element_blank(),
    axis.text.x = element_text(size = 25, face = "bold", color = "black"),
    strip.background = element_blank(),
    strip.text = element_markdown(size = 20),
    legend.position = "none"
  )

pheno_theme3 <- theme_classic2(base_size = 20) +
  theme(
    plot.title = element_markdown(hjust = 0.5, face = "bold"),
    axis.text.x = element_text(color = "black", size = 20),
    strip.background = element_blank(),
    strip.text = element_text(size = 20),
    legend.position = "none"
  )

Load spatially corrected data

# Load spatially corrected data with ionome
sp_corrected_file <- file.path(paths$intermediate, "PSU2022_spatially_corrected_both_treatments_with_ionome.csv")

sp_corrected_combined <- read.csv(sp_corrected_file) %>%
  mutate(
    Treatment = factor(Treatment, levels = c("+P", "-P")),
    genotype = factor(genotype, levels = c("CTRL", "Inv4m"))
  )

cat("Spatially corrected data dimensions:", dim(sp_corrected_combined), "\n")

## Spatially corrected data dimensions: 127 79

cat("Treatment levels:", levels(sp_corrected_combined$Treatment), "\n")

## Treatment levels: +P -P

cat("Genotype levels:", levels(sp_corrected_combined$genotype), "\n")

## Genotype levels: CTRL Inv4m

Prepare plotting data

pheno <- sp_corrected_combined %>%
  filter(genotype %in% c("CTRL", "Inv4m")) %>%
  mutate(genotype = factor(genotype, levels = c("CTRL", "Inv4m"))) %>%
  droplevels() %>%
  rename(Genotype = genotype) %>%
  mutate(
    Genotype_label = case_when(
      Genotype == "CTRL" ~ "CTRL",
      Genotype == "Inv4m" ~ "<i>Inv4m</i>",
      TRUE ~ as.character(Genotype)
    ),
    Genotype_label = factor(Genotype_label, levels = c("CTRL", "<i>Inv4m</i>"))
  )

cat("Analysis data dimensions:", dim(pheno), "\n")

## Analysis data dimensions: 64 80

---

2. Statistical Analysis

Identify phenotype variables

base_phenotypes <- c("PH", "DTA", "DTS", "STW40", "STW50", "STW60", "STWHV")
ear_phenotypes <- c("EW", "EL", "ED", "KRN", "CD", "TKW", "KW50", "TKN", "HI")

available_phenos <- intersect(
  c(base_phenotypes, ear_phenotypes),
  colnames(pheno)
)

cat("Available phenotypes for analysis:\n")

## Available phenotypes for analysis:

cat(paste(available_phenos, collapse = ", "), "\n")

## PH, DTA, DTS, STW40, STW50, STW60, STWHV, EW, EL, ED, KRN, CD, TKW, KW50, TKN, HI

vars <- available_phenos

Stage 1: Fit models and extract all effects with single FDR

# Fit models and extract ALL non-intercept coefficients
all_effects <- lapply(vars, function(x) {
  n_obs <- sum(!is.na(pheno[[x]]))
  
  if (n_obs < 20) {
    message("Skipping ", x, " - insufficient data")
    return(NULL)
  }
  
  model <- tryCatch(
    lm(as.formula(paste(x, "~ Genotype * Treatment")), data = pheno),
    error = function(e) {
      message("Model failed for ", x)
      return(NULL)
    }
  )
  
  if (is.null(model)) return(NULL)
  
  # Extract coefficients, excluding intercept
  coef_summary <- coef(summary(model))
  coef_summary <- coef_summary[!grepl("Intercept", rownames(coef_summary)), , drop = FALSE]
  
  if (nrow(coef_summary) == 0) return(NULL)
  
  out <- coef_summary %>%
    as.data.frame() %>%
    tibble::rownames_to_column("predictor")
  
  colnames(out)[3] <- "Std.Error"
  colnames(out)[5] <- "p.value"
  out$response <- x
  out %>% select(response, predictor, Estimate, Std.Error, p.value)
}) %>%
  bind_rows()

# Apply single FDR correction to all effects
all_effects <- all_effects %>%
  mutate(p.adjust = p.adjust(p.value, method = "fdr")) %>%
  arrange(p.adjust)

cat("\n=== All Model Effects (non-intercept) ===\n")

## 
## === All Model Effects (non-intercept) ===

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

## Total effects tested: 48

print(all_effects %>% head(20))

##    response                 predictor     Estimate  Std.Error      p.value
## 1     STW50               Treatment-P -18.75888523 1.58289202 1.598759e-16
## 2       DTS               Treatment-P   3.06210666 0.32199366 1.383332e-13
## 3       DTA               Treatment-P   3.43289929 0.36759849 2.668769e-13
## 4     STW40               Treatment-P  -6.51700782 1.00758227 2.137127e-08
## 5     STW60               Treatment-P -33.34346153 4.97340523 1.741324e-08
## 6     STWHV               Treatment-P -22.40748866 4.70104209 1.299584e-05
## 7       DTA             GenotypeInv4m  -1.47778357 0.36759849 1.648877e-04
## 8       DTS             GenotypeInv4m  -1.28418671 0.32199366 1.833645e-04
## 9        PH             GenotypeInv4m   5.56978074 1.48063541 3.853788e-04
## 10     KW50               Treatment-P  -1.80135139 0.53628018 1.521219e-03
## 11    STW40             GenotypeInv4m  -2.93451293 1.00758227 5.056342e-03
## 12       HI             GenotypeInv4m   0.07149298 0.02478209 5.805462e-03
## 13       CD GenotypeInv4m:Treatment-P  -2.61766738 0.97676410 1.000036e-02
## 14    STW50             GenotypeInv4m  -4.11284442 1.63844890 1.515412e-02
## 15    STWHV             GenotypeInv4m -11.65483187 4.70104209 1.609593e-02
## 16       HI               Treatment-P   0.05925653 0.02441592 1.894875e-02
## 17       EL             GenotypeInv4m   0.89364617 0.38817860 2.561592e-02
## 18      TKW             GenotypeInv4m  11.87098543 5.37554617 3.193065e-02
## 19       EL               Treatment-P   0.82722445 0.38244294 3.544784e-02
## 20    STW50 GenotypeInv4m:Treatment-P   4.61986258 2.22224856 4.248080e-02
##        p.adjust
## 1  7.674042e-15
## 2  3.319998e-12
## 3  4.270031e-12
## 4  2.051642e-07
## 5  2.051642e-07
## 6  1.039667e-04
## 7  1.100187e-03
## 8  1.100187e-03
## 9  2.055354e-03
## 10 7.301852e-03
## 11 2.206404e-02
## 12 2.322185e-02
## 13 3.692440e-02
## 14 5.150698e-02
## 15 5.150698e-02
## 16 5.684624e-02
## 17 7.232732e-02
## 18 8.514841e-02
## 19 8.955245e-02
## 20 1.019539e-01

Identify significant effects and traits with interactions

# Significant effects at FDR < 0.05
significant_effects <- all_effects %>%
  filter(p.adjust < 0.05)

cat("\n=== Significant Effects (FDR < 0.05) ===\n")

## 
## === Significant Effects (FDR < 0.05) ===

print(significant_effects)

##    response                 predictor     Estimate  Std.Error      p.value
## 1     STW50               Treatment-P -18.75888523 1.58289202 1.598759e-16
## 2       DTS               Treatment-P   3.06210666 0.32199366 1.383332e-13
## 3       DTA               Treatment-P   3.43289929 0.36759849 2.668769e-13
## 4     STW40               Treatment-P  -6.51700782 1.00758227 2.137127e-08
## 5     STW60               Treatment-P -33.34346153 4.97340523 1.741324e-08
## 6     STWHV               Treatment-P -22.40748866 4.70104209 1.299584e-05
## 7       DTA             GenotypeInv4m  -1.47778357 0.36759849 1.648877e-04
## 8       DTS             GenotypeInv4m  -1.28418671 0.32199366 1.833645e-04
## 9        PH             GenotypeInv4m   5.56978074 1.48063541 3.853788e-04
## 10     KW50               Treatment-P  -1.80135139 0.53628018 1.521219e-03
## 11    STW40             GenotypeInv4m  -2.93451293 1.00758227 5.056342e-03
## 12       HI             GenotypeInv4m   0.07149298 0.02478209 5.805462e-03
## 13       CD GenotypeInv4m:Treatment-P  -2.61766738 0.97676410 1.000036e-02
##        p.adjust
## 1  7.674042e-15
## 2  3.319998e-12
## 3  4.270031e-12
## 4  2.051642e-07
## 5  2.051642e-07
## 6  1.039667e-04
## 7  1.100187e-03
## 8  1.100187e-03
## 9  2.055354e-03
## 10 7.301852e-03
## 11 2.206404e-02
## 12 2.322185e-02
## 13 3.692440e-02

# Identify traits with significant interactions
traits_with_interaction <- significant_effects %>%
  filter(grepl(":", predictor)) %>%
  pull(response) %>%
  unique()

cat("\n=== Traits with Significant Interactions ===\n")

## 
## === Traits with Significant Interactions ===

if (length(traits_with_interaction) > 0) {
  cat(paste(traits_with_interaction, collapse = ", "), "\n")
} else {
  cat("None\n")
}

## CD

# Get all traits with any significant effect
traits_with_effects <- significant_effects %>%
  pull(response) %>%
  unique()

cat("\n=== All Traits with Significant Effects ===\n")

## 
## === All Traits with Significant Effects ===

cat(paste(traits_with_effects, collapse = ", "), "\n")

## STW50, DTS, DTA, STW40, STW60, STWHV, PH, KW50, HI, CD

cat("Total:", length(traits_with_effects), "traits\n")

## Total: 10 traits

Stage 2: Extract marginal or conditional contrasts for reporting

# Extract contrasts for significant traits only
contrasts_for_reporting <- lapply(traits_with_effects, function(x) {
  model <- lm(as.formula(paste(x, "~ Genotype + Treatment")), data = pheno)
  
  has_interaction <- x %in% traits_with_interaction
  
  if (has_interaction) {
    # Conditional contrasts for interaction traits
    treatment_cond <- summary(pairs(emmeans(model, ~ Treatment | Genotype)))
    genotype_cond <- summary(pairs(emmeans(model, ~ Genotype | Treatment)))
    
    bind_rows(
      treatment_cond %>%
        mutate(
          response = x,
          predictor = paste0("Treatment|", Genotype),
          contrast_type = "conditional"
        ) %>%
        select(response, predictor, contrast_type, estimate, SE, p.value),
      genotype_cond %>%
        mutate(
          response = x,
          predictor = paste0("Genotype|", Treatment),
          contrast_type = "conditional"
        ) %>%
        mutate(estimate = -estimate) %>%  # invert sign for -P vs +P I couldn't invert the contrast levels
        select(response, predictor, contrast_type, estimate, SE, p.value)
    )
  } else {
    # Marginal contrasts for non-interaction traits
    treatment_marg <- summary(pairs(emmeans(model, ~ Treatment)))
    genotype_marg <- summary(pairs(emmeans(model, ~ Genotype)))
    
    bind_rows(
      data.frame(
        response = x,
        predictor = "Treatment",
        contrast_type = "marginal",
        estimate = -treatment_marg$estimate, # invert sign for -P vs +P I couldn't invert the contrast levels
        SE = treatment_marg$SE,
        p.value = treatment_marg$p.value
      ),
      data.frame(
        response = x,
        predictor = "Genotype",
        contrast_type = "marginal",
        estimate = -genotype_marg$estimate, # invert sign for -P vs +P I couldn't invert the contrast levels
        SE = genotype_marg$SE,
        p.value = genotype_marg$p.value
      )
    )
  }
}) %>%
  bind_rows()

cat("\n=== Contrasts for Reporting ===\n")

## 
## === Contrasts for Reporting ===

print(contrasts_for_reporting)

##    response       predictor contrast_type     estimate         SE      p.value
## 1     STW50       Treatment      marginal -16.41494589 1.14466632 4.662928e-20
## 2     STW50        Genotype      marginal  -1.60148084 1.14042941 1.659616e-01
## 3       DTS       Treatment      marginal   3.41646064 0.23032280 6.674410e-22
## 4       DTS        Genotype      marginal  -0.92983272 0.23032280 1.534890e-04
## 5       DTA       Treatment      marginal   3.60161149 0.25869545 1.339050e-20
## 6       DTA        Genotype      marginal  -1.30907138 0.25869545 4.113753e-06
## 7     STW40       Treatment      marginal  -5.06257438 0.72462888 2.638380e-09
## 8     STW40        Genotype      marginal  -1.48007949 0.72462888 4.549988e-02
## 9     STW60       Treatment      marginal -27.79553850 3.51749524 2.067418e-10
## 10    STW60        Genotype      marginal   2.87358466 3.51024017 4.168069e-01
## 11    STWHV       Treatment      marginal -18.56909001 3.27431036 4.511543e-07
## 12    STWHV        Genotype      marginal  -7.81643322 3.27431036 2.020047e-02
## 13       PH       Treatment      marginal   3.42625576 1.04388584 1.706870e-03
## 14       PH        Genotype      marginal   6.40830361 1.04388584 6.819607e-08
## 15     KW50       Treatment      marginal  -1.86372143 0.36979252 6.475186e-06
## 16     KW50        Genotype      marginal   0.13864020 0.36926537 7.089149e-01
## 17       HI       Treatment      marginal   0.03424433 0.01717788 5.167509e-02
## 18       HI        Genotype      marginal   0.04572491 0.01715339 1.032085e-02
## 19       CD  Treatment|CTRL   conditional   1.53892221 0.51745913 4.515382e-03
## 20       CD Treatment|Inv4m   conditional   1.53892221 0.51745913 4.515382e-03
## 21       CD     Genotype|+P   conditional  -1.23073512 0.51672148 2.107391e-02
## 22       CD     Genotype|-P   conditional  -1.23073512 0.51672148 2.107391e-02

Calculate percentage effects

# Calculate reference means for each trait
reference_means_list <- lapply(traits_with_effects, function(x) {
  group_means <- pheno %>%
    filter(!is.na(.data[[x]])) %>%
    group_by(Genotype, Treatment) %>%
    summarise(mean_val = mean(.data[[x]], na.rm = TRUE), .groups = "drop")
  
  plus_p_marginal <- group_means %>%
    filter(Treatment == "+P") %>%
    summarise(mean = mean(mean_val)) %>%
    pull(mean)
  
  ctrl_marginal <- group_means %>%
    filter(Genotype == "CTRL") %>%
    summarise(mean = mean(mean_val)) %>%
    pull(mean)
  
  ctrl_minus_p <- group_means %>%
    filter(Genotype == "CTRL", Treatment == "-P") %>%
    pull(mean_val)
  
  if (length(ctrl_minus_p) == 0) ctrl_minus_p <- NA_real_
  
  data.frame(
    response = x,
    plus_p_marginal_mean = plus_p_marginal,
    ctrl_marginal_mean = ctrl_marginal,
    ctrl_minus_p_mean = ctrl_minus_p
  )
}) %>%
  bind_rows()

# Add percentages to contrasts
contrasts_with_pct <- contrasts_for_reporting %>%
  left_join(reference_means_list, by = "response") %>%
  mutate(
    reference_mean = case_when(
      predictor == "Treatment" ~ plus_p_marginal_mean,
      predictor == "Genotype" ~ ctrl_marginal_mean,
      grepl("Treatment\\|CTRL", predictor) ~ plus_p_marginal_mean,
      grepl("Treatment\\|Inv4m", predictor) ~ plus_p_marginal_mean,
      grepl("Genotype\\|\\+P", predictor) ~ ctrl_marginal_mean,
      grepl("Genotype\\|-P", predictor) ~ ctrl_minus_p_mean,
      TRUE ~ NA_real_
    ),
    estimate_pct = (estimate / reference_mean) * 100
  ) %>% select(-starts_with("ctrl_"), -starts_with("plus_"))

cat("\n=== Effects with Percentages ===\n")

## 
## === Effects with Percentages ===

print(contrasts_with_pct %>%
  select(response, predictor, contrast_type, estimate, estimate_pct, SE, p.value))

##    response       predictor contrast_type     estimate estimate_pct         SE
## 1     STW50       Treatment      marginal -16.41494589   -48.381951 1.14466632
## 2     STW50        Genotype      marginal  -1.60148084    -6.019516 1.14042941
## 3       DTS       Treatment      marginal   3.41646064     4.679297 0.23032280
## 4       DTS        Genotype      marginal  -0.92983272    -1.236720 0.23032280
## 5       DTA       Treatment      marginal   3.60161149     5.045611 0.25869545
## 6       DTA        Genotype      marginal  -1.30907138    -1.772935 0.25869545
## 7     STW40       Treatment      marginal  -5.06257438   -47.731076 0.72462888
## 8     STW40        Genotype      marginal  -1.48007949   -16.790072 0.72462888
## 9     STW60       Treatment      marginal -27.79553850   -36.790821 3.51749524
## 10    STW60        Genotype      marginal   2.87358466     4.763014 3.51024017
## 11    STWHV       Treatment      marginal -18.56909001   -18.355897 3.27431036
## 12    STWHV        Genotype      marginal  -7.81643322    -8.160384 3.27431036
## 13       PH       Treatment      marginal   3.42625576     2.083557 1.04388584
## 14       PH        Genotype      marginal   6.40830361     3.932642 1.04388584
## 15     KW50       Treatment      marginal  -1.86372143   -17.689031 0.36979252
## 16     KW50        Genotype      marginal   0.13864020     1.454181 0.36926537
## 17       HI       Treatment      marginal   0.03424433     8.277387 0.01717788
## 18       HI        Genotype      marginal   0.04572491    11.218326 0.01715339
## 19       CD  Treatment|CTRL   conditional   1.53892221     5.864469 0.51745913
## 20       CD Treatment|Inv4m   conditional   1.53892221     5.864469 0.51745913
## 21       CD     Genotype|+P   conditional  -1.23073512    -4.721431 0.51672148
## 22       CD     Genotype|-P   conditional  -1.23073512    -4.738758 0.51672148
##         p.value
## 1  4.662928e-20
## 2  1.659616e-01
## 3  6.674410e-22
## 4  1.534890e-04
## 5  1.339050e-20
## 6  4.113753e-06
## 7  2.638380e-09
## 8  4.549988e-02
## 9  2.067418e-10
## 10 4.168069e-01
## 11 4.511543e-07
## 12 2.020047e-02
## 13 1.706870e-03
## 14 6.819607e-08
## 15 6.475186e-06
## 16 7.089149e-01
## 17 5.167509e-02
## 18 1.032085e-02
## 19 4.515382e-03
## 20 4.515382e-03
## 21 2.107391e-02
## 22 2.107391e-02

Identify traits for visualization

# Traits with treatment effects (marginal or conditional)
treatment_traits <- contrasts_for_reporting %>%
  filter(grepl("Treatment", predictor)) %>%
  pull(response) %>%
  unique()

cat("\n=== Traits for Visualization ===\n")

## 
## === Traits for Visualization ===

cat(paste(treatment_traits, collapse = ", "), "\n")

## STW50, DTS, DTA, STW40, STW60, STWHV, PH, KW50, HI, CD

cat("Total:", length(treatment_traits), "traits\n")

## Total: 10 traits

Stage 3: 144-test family for figure brackets and downstream reporting

For each trait, all 9 contrasts spanning the 2x2 x phosphorus design are estimated from lm(y ~ Genotype * Treatment) via emmeans custom contrasts:

1. Treatment_at_CTRL — CTRL: -P vs +P
2. Treatment_at_Inv4m — Inv4m: -P vs +P
3. Genotype_at_plusP — +P: Inv4m vs CTRL
4. Genotype_at_minusP — -P: Inv4m vs CTRL
5. Diag_Inv4m_minusP_vs_CTRL_plusP — diagonal
6. Diag_Inv4m_plusP_vs_CTRL_minusP — diagonal
7. Treatment_marginal — -P vs +P pooled over Genotype
8. Genotype_marginal — Inv4m vs CTRL pooled over Treatment
9. Interaction — Genotype x Treatment interaction

All 144 contrasts (16 traits x 9) are FDR-adjusted in a single pool and saved to PSU2022_fig1_bracket_contrasts.csv. Backward-compatible columns contrast_level (“per_Genotype” or “marginal” for the original 3 contrast types; new values for the additional 6) and Genotype (CTRL/Inv4m for per_Genotype rows; NA otherwise) preserve filters used by the figure scripts.

# emmeans cell ordering for ~ Genotype * Treatment with levels
# Genotype = c("CTRL","Inv4m"), Treatment = c("+P","-P") is:
#   1: CTRL,+P   2: Inv4m,+P   3: CTRL,-P   4: Inv4m,-P
nine_contrasts <- list(
  "Treatment_at_CTRL"               = c(-1,  0,  1,  0),
  "Treatment_at_Inv4m"              = c( 0, -1,  0,  1),
  "Genotype_at_plusP"               = c(-1,  1,  0,  0),
  "Genotype_at_minusP"              = c( 0,  0, -1,  1),
  "Diag_Inv4m_minusP_vs_CTRL_plusP" = c(-1,  0,  0,  1),
  "Diag_Inv4m_plusP_vs_CTRL_minusP" = c( 0,  1, -1,  0),
  "Treatment_marginal"              = c(-0.5, -0.5,  0.5,  0.5),
  "Genotype_marginal"               = c(-0.5,  0.5, -0.5,  0.5),
  "Interaction"                     = c( 1,  -1,  -1,   1)
)

# Map each contrast name to (contrast_level, Genotype, group1, group2) for
# backward compatibility with existing figure filters.
contrast_metadata <- tibble::tribble(
  ~contrast_name,                    ~contrast_level,      ~Genotype, ~group1, ~group2,
  "Treatment_at_CTRL",               "per_Genotype",       "CTRL",    "+P",    "-P",
  "Treatment_at_Inv4m",              "per_Genotype",       "Inv4m",   "+P",    "-P",
  "Genotype_at_plusP",               "per_Treatment",       NA,        "CTRL",  "Inv4m",
  "Genotype_at_minusP",              "per_Treatment",       NA,        "CTRL",  "Inv4m",
  "Diag_Inv4m_minusP_vs_CTRL_plusP", "diagonal",            NA,        "CTRL+P","Inv4m-P",
  "Diag_Inv4m_plusP_vs_CTRL_minusP", "diagonal",            NA,        "CTRL-P","Inv4m+P",
  "Treatment_marginal",              "marginal",            NA,        "+P",    "-P",
  "Genotype_marginal",               "marginal_Genotype",   NA,        "CTRL",  "Inv4m",
  "Interaction",                     "interaction",         NA,        NA,      NA
)

contrasts <- bind_rows(
  lapply(vars, function(v) {
    if (sum(!is.na(pheno[[v]])) < 20) return(NULL)
    m <- tryCatch(
      lm(as.formula(paste(v, "~ Genotype * Treatment")), data = pheno),
      error = function(e) NULL)
    if (is.null(m)) return(NULL)
    em <- emmeans(m, ~ Genotype * Treatment)
    contrast(em, method = nine_contrasts) %>%
      as_tibble() %>%
      mutate(trait = v, contrast_name = as.character(contrast))
  })
) %>%
  rename(p = p.value, statistic = t.ratio) %>%
  rstatix::adjust_pvalue(method = "fdr") %>%
  left_join(contrast_metadata, by = "contrast_name")

cat("Phenotype contrasts — 144-test family, BH FDR across all 144:\n")

## Phenotype contrasts — 144-test family, BH FDR across all 144:

contrasts %>%
  select(trait, contrast_name, contrast_level, statistic, p, p.adj) %>%
  arrange(p.adj) %>% as_tibble() %>% print(n = 30)

## # A tibble: 144 × 6
##    trait contrast_name                contrast_level statistic        p    p.adj
##    <chr> <chr>                        <chr>              <dbl>    <dbl>    <dbl>
##  1 DTS   Treatment_marginal           marginal           15.0  5.66e-22 8.14e-20
##  2 DTA   Treatment_marginal           marginal           13.9  2.35e-20 1.13e-18
##  3 STW50 Treatment_marginal           marginal          -14.8  1.79e-20 1.13e-18
##  4 DTS   Diag_Inv4m_plusP_vs_CTRL_mi… diagonal          -13.5  7.76e-20 2.79e-18
##  5 DTA   Diag_Inv4m_plusP_vs_CTRL_mi… diagonal          -13.4  1.24e-19 3.57e-18
##  6 DTS   Treatment_at_Inv4m           per_Genotype       11.7  3.87e-17 9.29e-16
##  7 STW50 Treatment_at_CTRL            per_Genotype      -11.9  1.60e-16 3.29e-15
##  8 STW50 Diag_Inv4m_minusP_vs_CTRL_p… diagonal          -11.7  2.61e-16 4.69e-15
##  9 DTA   Treatment_at_Inv4m           per_Genotype       10.3  8.12e-15 1.30e-13
## 10 DTS   Treatment_at_CTRL            per_Genotype        9.51 1.38e-13 1.99e-12
## 11 DTA   Treatment_at_CTRL            per_Genotype        9.34 2.67e-13 3.49e-12
## 12 STW50 Diag_Inv4m_plusP_vs_CTRL_mi… diagonal            9.25 1.18e-12 1.42e-11
## 13 STW50 Treatment_at_Inv4m           per_Genotype       -9.06 2.32e-12 2.57e-11
## 14 DTS   Diag_Inv4m_minusP_vs_CTRL_p… diagonal            7.72 1.46e-10 1.40e- 9
## 15 STW60 Treatment_marginal           marginal           -8.05 1.39e-10 1.40e- 9
## 16 STW40 Treatment_marginal           marginal           -7.20 1.26e- 9 1.13e- 8
## 17 PH    Diag_Inv4m_minusP_vs_CTRL_p… diagonal            6.64 1.02e- 8 8.61e- 8
## 18 STW40 Diag_Inv4m_minusP_vs_CTRL_p… diagonal           -6.54 1.62e- 8 1.29e- 7
## 19 STW60 Treatment_at_CTRL            per_Genotype       -6.70 1.74e- 8 1.32e- 7
## 20 STW40 Treatment_at_CTRL            per_Genotype       -6.47 2.14e- 8 1.54e- 7
## 21 DTA   Diag_Inv4m_minusP_vs_CTRL_p… diagonal            6.24 4.93e- 8 3.38e- 7
## 22 PH    Genotype_marginal            marginal_Geno…      6.12 7.71e- 8 5.05e- 7
## 23 STW60 Diag_Inv4m_plusP_vs_CTRL_mi… diagonal            6.00 2.22e- 7 1.39e- 6
## 24 STWHV Treatment_marginal           marginal           -5.72 3.92e- 7 2.35e- 6
## 25 STWHV Diag_Inv4m_minusP_vs_CTRL_p… diagonal           -5.67 4.79e- 7 2.76e- 6
## 26 STW60 Diag_Inv4m_minusP_vs_CTRL_p… diagonal           -5.37 2.02e- 6 1.12e- 5
## 27 DTA   Genotype_marginal            marginal_Geno…     -5.04 4.63e- 6 2.47e- 5
## 28 PH    Genotype_at_minusP           per_Treatment       4.89 7.77e- 6 4.00e- 5
## 29 KW50  Treatment_marginal           marginal           -4.98 8.19e- 6 4.07e- 5
## 30 STWHV Treatment_at_CTRL            per_Genotype       -4.77 1.30e- 5 6.24e- 5
## # ℹ 114 more rows

# CTRL +P mean as common reference baseline for all percent-change calculations
ctrl_plus_p_means <- vapply(vars, function(v) {
  mean(pheno[[v]][pheno$Genotype == "CTRL" & pheno$Treatment == "+P"], na.rm = TRUE)
}, numeric(1)) |> (\(x) data.frame(trait = names(x), ctrl_plus_p_mean = x))()

bracket_contrasts <- contrasts %>%
  transmute(
    trait,
    Genotype,
    contrast_level,
    contrast_name,
    group1,
    group2,
    estimate,
    SE,
    statistic,
    p,
    p.adj,
    p.adj.signif = rstatix::add_significance(
      tibble(p.adj = p.adj), "p.adj"
    )$p.adj.signif
  ) %>%
  left_join(ctrl_plus_p_means, by = "trait") %>%
  mutate(
    estimate_pct = estimate / ctrl_plus_p_mean * 100,
    SE_pct       = SE       / ctrl_plus_p_mean * 100
  ) %>%
  select(-ctrl_plus_p_mean)

cat("Per-genotype +-P contrasts (FDR over 144 tests):\n")

## Per-genotype +-P contrasts (FDR over 144 tests):

bracket_contrasts %>%
  filter(contrast_level == "per_Genotype") %>%
  select(trait, Genotype, estimate_pct, SE_pct, p.adj, p.adj.signif) %>%
  arrange(p.adj) %>% as_tibble() %>% print(n = Inf)

## # A tibble: 32 × 6
##    trait Genotype estimate_pct SE_pct    p.adj p.adj.signif
##    <chr> <chr>           <dbl>  <dbl>    <dbl> <chr>       
##  1 DTS   Inv4m           5.12   0.437 9.29e-16 ****        
##  2 STW50 CTRL          -52.1    4.40  3.29e-15 ****        
##  3 DTA   Inv4m           5.23   0.510 1.30e-13 ****        
##  4 DTS   CTRL            4.16   0.437 1.99e-12 ****        
##  5 DTA   CTRL            4.76   0.510 3.49e-12 ****        
##  6 STW50 Inv4m         -39.3    4.33  2.57e-11 ****        
##  7 STW60 CTRL          -43.3    6.46  1.32e- 7 ****        
##  8 STW40 CTRL          -54.0    8.35  1.54e- 7 ****        
##  9 STWHV CTRL          -20.9    4.39  6.24e- 5 ****        
## 10 STW60 Inv4m         -29.3    6.29  1.13e- 4 ***         
## 11 CD    Inv4m         -10.7    2.60  6.17e- 4 ***         
## 12 STW40 Inv4m         -30.3    8.21  1.88e- 3 **          
## 13 KW50  Inv4m         -18.4    4.99  2.05e- 3 **          
## 14 KW50  CTRL          -17.3    5.14  5.09e- 3 **          
## 15 STWHV Inv4m         -14.0    4.24  5.43e- 3 **          
## 16 PH    Inv4m           2.64   0.916 1.55e- 2 *           
## 17 HI    CTRL           15.7    6.46  4.43e- 2 *           
## 18 EL    CTRL            5.61   2.59  7.51e- 2 ns          
## 19 TKN   CTRL           20.4   10.3   9.64e- 2 ns          
## 20 ED    Inv4m          -5.26   2.78  1.13e- 1 ns          
## 21 PH    CTRL            1.60   0.916 1.45e- 1 ns          
## 22 TKW   Inv4m         -12.6    7.70  1.76e- 1 ns          
## 23 EW    Inv4m         -11.3    7.07  1.87e- 1 ns          
## 24 TKN   Inv4m          14.2    9.95  2.42e- 1 ns          
## 25 KRN   Inv4m          -2.80   2.19  3.01e- 1 ns          
## 26 EL    Inv4m          -1.58   2.51  6.40e- 1 ns          
## 27 HI    Inv4m           2.83   6.27  7.41e- 1 ns          
## 28 EW    CTRL           -2.88   7.29  7.70e- 1 ns          
## 29 CD    CTRL           -0.729  2.68  8.46e- 1 ns          
## 30 KRN   CTRL            0.550  2.25  8.62e- 1 ns          
## 31 ED    CTRL            0.200  2.86  9.44e- 1 ns          
## 32 TKW   CTRL           -0.562  7.94  9.44e- 1 ns

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

cat("Saved", nrow(bracket_contrasts),
    "contrasts to PSU2022_fig1_bracket_contrasts.csv (expected 144)\n")

## Saved 144 contrasts to PSU2022_fig1_bracket_contrasts.csv (expected 144)

stopifnot(nrow(bracket_contrasts) == 144)

---

3. Pairwise Comparisons for Visualization

to_t_test <- pheno %>%
  select(plotId, Genotype, Treatment, all_of(vars)) %>%
  pivot_longer(
    cols = all_of(vars),
    names_to = "trait",
    values_to = "value"
  ) %>%
  filter(trait %in% treatment_traits)

brackets_d <- bracket_contrasts

if (length(treatment_traits) > 0) {
  y_pos <- to_t_test %>%
    group_by(trait, Genotype) %>%
    summarise(y_max   = max(value, na.rm = TRUE),
              y_range = diff(range(value, na.rm = TRUE)),
              .groups = "drop") %>%
    mutate(y.position = y_max + 0.08 * y_range)

  treatment_pairwise <- brackets_d %>%
    filter(contrast_level == "per_Genotype",
           trait %in% treatment_traits) %>%
    mutate(
      Genotype_label = factor(
        case_when(
          Genotype == "CTRL"  ~ "CTRL",
          Genotype == "Inv4m" ~ "<i>Inv4m</i>",
          TRUE ~ as.character(Genotype)
        ),
        levels = c("CTRL", "<i>Inv4m</i>")
      )
    ) %>%
    left_join(y_pos, by = c("trait", "Genotype")) %>%
    add_x_position(x = "Treatment")

  cat("\n=== Approach D contrast summary ===\n")
  print(as_tibble(treatment_pairwise) %>% select(trait, Genotype, statistic, p.adj, p.adj.signif), n = 20)
} else {
  treatment_pairwise <- NULL
  cat("No traits with treatment effects for visualization\n")
}

## 
## === Approach D contrast summary ===
## # A tibble: 20 × 5
##    trait Genotype statistic    p.adj p.adj.signif
##    <chr> <chr>        <dbl>    <dbl> <chr>       
##  1 PH    CTRL         1.75  1.45e- 1 ns          
##  2 PH    Inv4m        2.88  1.55e- 2 *           
##  3 DTA   CTRL         9.34  3.49e-12 ****        
##  4 DTA   Inv4m       10.3   1.30e-13 ****        
##  5 DTS   CTRL         9.51  1.99e-12 ****        
##  6 DTS   Inv4m       11.7   9.29e-16 ****        
##  7 STW40 CTRL        -6.47  1.54e- 7 ****        
##  8 STW40 Inv4m       -3.69  1.88e- 3 **          
##  9 STW50 CTRL       -11.9   3.29e-15 ****        
## 10 STW50 Inv4m       -9.06  2.57e-11 ****        
## 11 STW60 CTRL        -6.70  1.32e- 7 ****        
## 12 STW60 Inv4m       -4.65  1.13e- 4 ***         
## 13 STWHV CTRL        -4.77  6.24e- 5 ****        
## 14 STWHV Inv4m       -3.30  5.43e- 3 **          
## 15 CD    CTRL        -0.272 8.46e- 1 ns          
## 16 CD    Inv4m       -4.13  6.17e- 4 ***         
## 17 KW50  CTRL        -3.36  5.09e- 3 **          
## 18 KW50  Inv4m       -3.69  2.05e- 3 **          
## 19 HI    CTRL         2.43  4.43e- 2 *           
## 20 HI    Inv4m        0.452 7.41e- 1 ns

---

4. Plotting Function

plot_treatment_effect <- function(trait_name, 
                                   trait_title, 
                                   trait_ylab,
                                   y_limits = NULL,
                                   y_breaks = NULL,
                                   bracket_nudge = 0.05) {
  
  t_test_data <- NULL
  if (!is.null(treatment_pairwise)) {
    t_test_data <- treatment_pairwise %>%
      filter(trait == trait_name)
  }
  
  p <- pheno %>%
    ggplot(aes(x = Treatment, y = .data[[trait_name]], col = Treatment)) +
    ggtitle(trait_title) +
    ylab(trait_ylab) +
    geom_boxplot(width = 0.25, linewidth = 1, alpha = 0) %>% 
      with_shadow(colour = "black", x_offset = 0, y_offset = 0, sigma = 1) +
    geom_quasirandom(size = 2) %>% 
      with_shadow(colour = "black", x_offset = 0, y_offset = 0, sigma = 1) +
    scale_color_manual(values = pal) +
    facet_wrap(. ~ Genotype_label) +
    pheno_theme2 +
    theme(strip.text = element_markdown())
  
  if (!is.null(t_test_data) && nrow(t_test_data) > 0) {
    p <- p + stat_pvalue_manual(
      t_test_data,
      size = 10,
      label = "p.adj.signif"  ,
      bracket.size = 0.8,
      hide.ns = FALSE
    )
  }
  
  if (!is.null(y_limits)) {
    p <- p + coord_cartesian(ylim = y_limits)
  }
  
  if (!is.null(y_breaks)) {
    p <- p + scale_y_continuous(breaks = y_breaks)
  }
  
  return(p)
}

---

5. Individual Trait Plots

Days to Anthesis

if ("DTA" %in% treatment_traits) {
  dta_plot <- plot_treatment_effect(
    trait_name = "DTA",
    trait_title = "Anthesis",
    trait_ylab = "Days",
    y_limits = c(69, 80),
    y_breaks = 69:80,
    bracket_nudge = 0.08
  )
  print(dta_plot)
}

Days to Silking

if ("DTS" %in% treatment_traits) {
  dts_plot <- plot_treatment_effect(
    trait_name = "DTS",
    trait_title = "Silking",
    trait_ylab = "Days",
    y_limits = c(69, 80),
    y_breaks = 69:80,
    bracket_nudge = 0
  )
  print(dts_plot)
}

Plant Height

if ("PH" %in% treatment_traits) {
  ph_plot <- plot_treatment_effect(
    trait_name = "PH",
    trait_title = "Plant Height",
    trait_ylab = "cm",
    y_limits = c(140, 180),
    y_breaks = seq(140, 180, by = 10),
    bracket_nudge = 0.08
  )
  print(ph_plot)
}

50 Kernel Weight

if ("KW50" %in% treatment_traits) {
  kw50_plot <- plot_treatment_effect(
    trait_name = "KW50",
    trait_title = "50 Kernel Weight",
    trait_ylab = "g",
    y_limits = c(7, 17),
    y_breaks = 7:17,
    bracket_nudge = 0.08
  )
  print(kw50_plot)
}

Cob Diameter

if ("CD" %in% treatment_traits) {
  cd_plot <- plot_treatment_effect(
    trait_name = "CD",
    trait_title = "Cob Diameter",
    trait_ylab = "cm",
    y_limits = c(22, 31),
    y_breaks = 22:31,
    bracket_nudge = 0.08
  )
  print(cd_plot)
}

Stover Weight at 40 DAP

if ("STW40" %in% treatment_traits) {
  stw40_plot <- plot_treatment_effect(
    trait_name = "STW40",
    trait_title = "Stover Weight <br> at 40 DAP",
    trait_ylab = "g",
    bracket_nudge = 0.08
  )
  print(stw40_plot)
}

Stover Weight at 50 DAP

if ("STW50" %in% treatment_traits) {
  stw50_plot <- plot_treatment_effect(
    trait_name = "STW50",
    trait_title = "Stover Weight <br> at 50 DAP",
    trait_ylab = "g",
    bracket_nudge = 0.08
  )
  print(stw50_plot)
}

Stover Weight at 60 DAP

if ("STW60" %in% treatment_traits) {
  stw60_plot <- plot_treatment_effect(
    trait_name = "STW60",
    trait_title = "Stover Weight <br> at 60 DAP",
    trait_ylab = "g",
    bracket_nudge = 0.08
  )
  print(stw60_plot)
}

Stover Weight at Harvest

if ("STWHV" %in% treatment_traits) {
  stwhv_plot <- plot_treatment_effect(
    trait_name = "STWHV",
    trait_title = "Stover Dry Weight <br> at Harvest",
    trait_ylab = "g",
    y_limits = c(50, 175),
    bracket_nudge = 0.08
  )
  print(stwhv_plot)
}

---

6. Stover Dry Weight Time Course

Prepare time-series data

dw_raw <- pheno %>%
  mutate(STW121 = STWHV) %>%
  select(
    plotId, 
    Genotype, 
    Genotype_label,
    Treatment, 
    starts_with("STW")
  ) %>%
  pivot_longer(
    cols = c("STW40", "STW50", "STW60", "STW121"),
    names_to = "DAP_string",
    values_to = "STW"
  ) %>%
  mutate(
    DAP = as.integer(gsub("STW", "", DAP_string)),
    Genotype = factor(Genotype, levels = c("CTRL", "Inv4m")),
    Genotype_label = factor(Genotype_label, levels = c("CTRL", "<i>Inv4m</i>"))
  ) %>%
  filter(!is.na(STW)) %>%
  droplevels()

Time-point specific t-tests

y_pos_ts <- dw_raw %>%
  filter(DAP > 1) %>%
  group_by(Genotype, DAP) %>%
  summarise(y_max   = max(STW, na.rm = TRUE),
            y_range = diff(range(STW, na.rm = TRUE)),
            .groups = "drop") %>%
  mutate(y.position = y_max + 0.08 * y_range,
         DAP        = factor(DAP))

treatment_timeseries_test <- brackets_d %>%
  filter(contrast_level == "per_Genotype",
         trait %in% c("STW40", "STW50", "STW60", "STWHV")) %>%
  mutate(
    DAP = factor(case_when(
      trait == "STW40" ~ 40,
      trait == "STW50" ~ 50,
      trait == "STW60" ~ 60,
      trait == "STWHV" ~ 121
    )),
    Genotype_label = factor(
      case_when(
        Genotype == "CTRL"  ~ "CTRL",
        Genotype == "Inv4m" ~ "<i>Inv4m</i>",
        TRUE ~ as.character(Genotype)
      ),
      levels = c("CTRL", "<i>Inv4m</i>")
    )
  ) %>%
  left_join(y_pos_ts, by = c("Genotype", "DAP")) %>%
  mutate(
    x      = as.numeric(DAP),
    xmin   = as.numeric(DAP) - 0.25,
    xmax   = as.numeric(DAP) + 0.25,
    group1 = as.character(DAP),
    group2 = as.character(DAP)
  )

cat("\n=== Time-series contrast summary (Approach D) ===\n")

## 
## === Time-series contrast summary (Approach D) ===

print(as_tibble(treatment_timeseries_test) %>%
      select(Genotype, DAP, statistic, p, p.adj, p.adj.signif))

## # A tibble: 8 × 6
##   Genotype DAP   statistic        p    p.adj p.adj.signif
##   <chr>    <fct>     <dbl>    <dbl>    <dbl> <chr>       
## 1 CTRL     40        -6.47 2.14e- 8 1.54e- 7 ****        
## 2 Inv4m    40        -3.69 4.95e- 4 1.88e- 3 **          
## 3 CTRL     50       -11.9  1.60e-16 3.29e-15 ****        
## 4 Inv4m    50        -9.06 2.32e-12 2.57e-11 ****        
## 5 CTRL     60        -6.70 1.74e- 8 1.32e- 7 ****        
## 6 Inv4m    60        -4.65 2.42e- 5 1.13e- 4 ***         
## 7 CTRL     121       -4.77 1.30e- 5 6.24e- 5 ****        
## 8 Inv4m    121       -3.30 1.66e- 3 5.43e- 3 **

Visualize stover weight trajectory

pd <- position_dodge(1)

p_dw <- dw_raw %>%
  ungroup() %>%
  filter(DAP > 1) %>%
  ggplot(aes(x = as.factor(DAP), y = STW, col = Treatment)) +
  ggtitle("Stover Dry Weight") +
  labs(
    x = "Days After Planting",
    y = "g",
    color= ""
  ) +
  geom_boxplot(
    width = 0.25, 
    linewidth = 1, 
    alpha = 0, 
    position = pd
  ) %>% 
    with_shadow(colour = "black", x_offset = 0, y_offset = 0, sigma = 1) +
  geom_quasirandom(
    size = 2, 
    width = 0.25, 
    dodge.width = 1
  ) %>% 
    with_shadow(colour = "black", x_offset = 0, y_offset = 0, sigma = 1) +
  facet_wrap(. ~ Genotype_label) +
  stat_pvalue_manual(
    treatment_timeseries_test,
    size = 10,
    bracket.size = 0.8,
    tip.length = 0.01,
    hide.ns = TRUE
  ) +
  scale_color_manual(values = pal) +
  scale_y_continuous(expand = expansion(mult = c(0.05, 0.1))) +
  pheno_theme3 +
  theme(
    legend.position = c(0.1, 0.9),
    legend.background = element_rect(fill = "transparent"),
    legend.text = element_text(face = "bold", size = 20),
    strip.text = element_markdown()
  )

print(p_dw)

---

7. Biomass Trajectory Analysis

Extract stover weight effects

# Get marginal contrasts for stover weights
stover_contrasts <- contrasts_for_reporting %>%
  filter(
    grepl("^STW", response),
    predictor %in% c("Genotype", "Treatment")
  )

# Join with reference means and calculate percentages
stover_effects <- stover_contrasts %>%
  left_join(reference_means_list, by = "response") %>%
  mutate(
    reference_mean = if_else(
      predictor == "Treatment",
      plus_p_marginal_mean,
      ctrl_marginal_mean
    ),
    ci_lower = estimate - 1.96 * SE,
    ci_upper = estimate + 1.96 * SE,
    estimate_pct = (estimate / reference_mean) * 100,
    pct_ci_lower = (ci_lower / reference_mean) * 100,
    pct_ci_upper = (ci_upper / reference_mean) * 100,
    DAP = case_when(
      response == "STW40" ~ "40",
      response == "STW50" ~ "50",
      response == "STW60" ~ "60",
      response == "STWHV" ~ "121",
      TRUE ~ NA_character_
    ),
    DAP = factor(DAP, levels = c("40", "50", "60", "121")),
    effect_type = factor(predictor, levels = c("Genotype", "Treatment"))
  ) %>%
  mutate(DAP_numeric = as.numeric(as.character(DAP))) %>%
  arrange(DAP, effect_type)

cat("\n=== Stover Effects Summary ===\n")

## 
## === Stover Effects Summary ===

print(stover_effects %>%
  select(response, DAP, predictor, estimate_pct, p.value))

##   response DAP predictor estimate_pct      p.value
## 1    STW40  40  Genotype   -16.790072 4.549988e-02
## 2    STW40  40 Treatment   -47.731076 2.638380e-09
## 3    STW50  50  Genotype    -6.019516 1.659616e-01
## 4    STW50  50 Treatment   -48.381951 4.662928e-20
## 5    STW60  60  Genotype     4.763014 4.168069e-01
## 6    STW60  60 Treatment   -36.790821 2.067418e-10
## 7    STWHV 121  Genotype    -8.160384 2.020047e-02
## 8    STWHV 121 Treatment   -18.355897 4.511543e-07

Fit regressions

minus_p_data <- stover_effects %>% filter(predictor == "Treatment")
inv4m_data   <- stover_effects %>%
  filter(predictor == "Genotype") %>%
  mutate(DAP_numeric = as.numeric(as.character(DAP)))

lm_minus_p <- lm(estimate_pct ~ DAP_numeric, data = minus_p_data)
lm_inv4m   <- lm(estimate_pct ~ DAP_numeric, data = inv4m_data)

r_squared_minus_p <- summary(lm_minus_p)$r.squared
p_value_minus_p   <- summary(lm_minus_p)$coefficients[2, 4]
r_squared_inv4m   <- summary(lm_inv4m)$r.squared
p_value_inv4m     <- summary(lm_inv4m)$coefficients[2, 4]

regression_label_minus_p <- sprintf(
  "italic(R)^2 == %.3f~~italic(p) == %.3f",
  r_squared_minus_p, p_value_minus_p
)
regression_label_inv4m <- sprintf(
  "italic(R)^2 == %.3f~~italic(p) == %.3f",
  r_squared_inv4m, p_value_inv4m
)

Plot biomass trajectory

 p_stw_lag <- stover_effects %>%
  mutate(DAP_numeric = as.numeric(as.character(DAP))) %>%
  ggplot(aes(x = DAP_numeric, y = -estimate_pct, group = effect_type)) +
  geom_hline(yintercept = 0, linetype = "solid", color = "gray50", linewidth = 0.8) +
  geom_smooth(
    data = minus_p_data,
    method = "lm",
    se = FALSE,
    linewidth = 1,
    linetype = "dashed",
    color = "black"
  ) +
  geom_smooth(
    data = inv4m_data,
    method = "lm",
    se = FALSE,
    linewidth = 1,
    linetype = "dotted",
    color = "black"
  ) +
  geom_pointrange(
    aes(ymin = -pct_ci_lower, ymax = -pct_ci_upper, shape = effect_type),
    color = "black",
    size = 1,
    linewidth = 1,
    fatten = 5,
    position = position_dodge(width = 3)
  ) +
  annotate(
    "text",
    x = 80,
    y = 40,
    label = regression_label_minus_p,
    parse = TRUE,
    size = 6,
    hjust = 0
  ) +
  annotate(
    "text",
    x = 80,
    y = 10,
    label = regression_label_inv4m,
    parse = TRUE,
    size = 6,
    hjust = 0
  ) +
  scale_x_continuous(
    breaks = c(40, 50, 60, 121),
    labels = c("40", "50", "60", "121")
  ) +
  scale_shape_manual(
    values = c(16, 18),
    labels = c("<i>Inv4m</i> vs CTRL", "-P vs +P")
  ) +
  guides(shape = guide_legend(reverse = TRUE)) +
  labs(
    title = "Stover Biomass Lag",
    x = "Days After Planting",
    y = "% reference mean",
    shape = ""
  ) +
  theme_classic2(base_size = 30) +
  theme(
    plot.title = element_markdown(hjust = 0.5, face = "bold",size=25),
    axis.text.x = element_text(color = "black"),
        legend.background = element_rect(fill = "transparent"),
    legend.position = c(0.85,0.95),
    legend.text = element_markdown()
  )

#the odd value at DAP60 is real
summary(pairs(emmeans(lm(data=pheno,  STW60 ~ Genotype * Treatment), ~ Genotype), reverse = TRUE))

##  contrast     estimate   SE df t.ratio p.value
##  Inv4m - CTRL      2.5 3.47 50   0.720  0.4750
## 
## Results are averaged over the levels of: Treatment

print( p_stw_lag)

---

8. Combined Panels

Main phenotype panel

# Load growth curve plots
plot_growth_params_file <- file.path(paths$intermediate, "plot_growth_params.rds")
plot_fitted_curves_file <- file.path(paths$intermediate, "plot_fitted_curves.rds")

if (file.exists(plot_growth_params_file)) {
  p_gc <- readRDS(plot_growth_params_file)
  cat("Loaded growth parameter plot\n")
} else {
  warning("Growth parameter plot not found. Run PSU2022_growthcurves.Rmd first.")
}

## Loaded growth parameter plot

if (file.exists(plot_fitted_curves_file)) {
  p_fitted <- readRDS(plot_fitted_curves_file)
  cat("Loaded fitted curves plot\n")
} else {
  warning("Fitted curves plot not found. Run PSU2022_growthcurves.Rmd first.")
}

## Loaded fitted curves plot

plot_list <- list()

if (exists("dta_plot")) plot_list$dta <- dta_plot
if (exists("dts_plot")) plot_list$dts <- dts_plot
# if (exists("ph_plot")) plot_list$ph <- ph_plot
if (exists("kw50_plot")) plot_list$kw50 <- kw50_plot
if (exists("cd_plot")) plot_list$cd <- cd_plot

if (length(plot_list) >= 4) {
  main_plots <- plot_list[1:min(8, length(plot_list))]
  
  combined_plot <- ggarrange(
    plotlist = main_plots,
    ncol = 4,
    nrow = 1,
    labels = LETTERS[1:length(main_plots)+2],
    font.label = list(size = 30)
  )
  
  print(combined_plot)
} else if (length(plot_list) > 0) {
  combined_plot <- ggarrange(
    plotlist = plot_list,
    ncol = 2,
    nrow = ceiling(length(plot_list) / 2),
    labels = LETTERS[1:length(plot_list)+2],
    font.label = list(size = 30)
  )
  
  print(combined_plot)
}

Combined stover weight panel

# Only create panel if p_fitted exists
if (exists("p_fitted") && exists("p_dw")) {
  p_stw <- ggarrange(
    p_dw,
    p_fitted,
    nrow = 1,
    ncol = 2,
    align="hv",
    labels = LETTERS[(length(plot_list)+3):(length(plot_list)+4)],
    font.label = list(size = 30)
  )

  print(p_stw)
} else {
  cat("Skipping stover panel - p_fitted or p_dw not available. Run PSU2022_growthcurves.Rmd first.\n")
}

### Combined figure with time series

# Create comprehensive figure including time series
if (exists("combined_plot") && exists("p_stw")) {
  full_figure <- ggarrange(
    combined_plot,
    p_stw,
    ncol = 1,
    nrow = 2,
    # heights = c(1, 1.2),
    #labels = c("", "E"),
    font.label = list(size = 30)
  )

  print(full_figure)

  ggsave(
    filename = file.path(paths$figures, "pheno.png"),
    plot = full_figure,
    height = 10,
    width = 16,
    units = "in",
    dpi = 150
  )
} else {
  cat("Skipping full figure - required plots not available.\n")
}

Supplementary Figure S2: Stover Dry Weight Growth Curves and Derived Parameters

if (exists("p_gc") && exists("p_stw_lag")) {
  growth_supp_figure <- ggarrange(
    p_gc,
    p_stw_lag,
    nrow = 2, ncol = 1,
    labels = LETTERS[1:2],
    font.label = list(size = 35),
    heights = c(0.75, 1)
  )

  print(growth_supp_figure)

  ggsave(
    filename = file.path(paths$figures, "biomass.png"),
    plot = growth_supp_figure,
    height = 16, width = 12,
    units = "in", dpi = 150
  )
} else {
  cat("Skipping Figure S2 — run PSU2022_growthcurves.Rmd first.\n")
}

Inv4m differences in Anthesis and Plant Height

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

# Inv4m - CTRL contrasts at each Treatment level for DTA and PH
geno_contrasts <- lapply(c("DTA", "PH"), function(v) {
  m <- lm(as.formula(paste(v, "~ Genotype * Treatment")), data = pheno)
  pairs(emmeans(m, ~ Genotype | Treatment), reverse = TRUE) %>%
    as_tibble() %>%
    mutate(trait = v)
}) %>%
  bind_rows() %>%
  rename(p = p.value, statistic = t.ratio) %>%
  mutate(group1 = "CTRL", group2 = "<i>Inv4m</i>") %>%
  rstatix::adjust_pvalue(method = "fdr") %>%
  rstatix::add_significance("p.adj")

cat("=== Inv4m - CTRL contrasts at each Treatment level ===\n")

## === Inv4m - CTRL contrasts at each Treatment level ===

print(as_tibble(geno_contrasts) %>%
  select(trait, Treatment, estimate, SE, statistic, p, p.adj, p.adj.signif))

## # A tibble: 4 × 8
##   trait Treatment estimate    SE statistic          p     p.adj p.adj.signif
##   <chr> <fct>        <dbl> <dbl>     <dbl>      <dbl>     <dbl> <chr>       
## 1 DTA   +P           -1.48 0.368     -4.02 0.000165   0.000330  ***         
## 2 DTA   -P           -1.14 0.368     -3.10 0.00293    0.00293   **          
## 3 PH    +P            5.57 1.48       3.76 0.000385   0.000514  ***         
## 4 PH    -P            7.25 1.48       4.89 0.00000777 0.0000311 ****

supp_geno_panel <- function(trait_name, ylab_text, y_limits, y_breaks) {
  d <- pheno %>% select(Genotype, Genotype_label, Treatment,
                         value = all_of(trait_name)) %>%
    filter(!is.na(value))

  s <- geno_contrasts %>%
    filter(trait == trait_name) %>%
    group_by(Treatment) %>%
    mutate(y.position = max(d$value[d$Treatment == Treatment], na.rm = TRUE) +
                        0.08 * diff(range(d$value[d$Treatment == Treatment], na.rm = TRUE))) %>%
    ungroup()

  ggplot(d, aes(x = Genotype_label, y = value, colour = Genotype)) +
    ggtitle(trait_name) +
    ylab(ylab_text) +
    geom_boxplot(width = 0.25, linewidth = 1, alpha = 0) %>%
      with_shadow(colour = "black", x_offset = 0, y_offset = 0, sigma = 1) +
    geom_quasirandom(size = 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", "<i>Inv4m</i>" = "<i>Inv4m</i>")) +
    facet_wrap(~ Treatment) +
    stat_pvalue_manual(s, label = "p.adj.signif",
                       size = 10, bracket.size = 0.8, hide.ns = FALSE) +
    scale_y_continuous(breaks = y_breaks, limits = y_limits,
                       expand = expansion(mult = 0)) +
    pheno_theme2 +
    theme(
      axis.text.x  = element_markdown(size = 25, face = "bold", color = "black"),
      strip.text   = element_markdown(size = 25, face = "bold")
    )
}

p_dta <- supp_geno_panel("DTA", "Days", c(69, 80), 69:80)
p_ph  <- supp_geno_panel("PH",  "cm",   c(150, 185), seq(150, 185, by = 5))

supp_global_fdr_plot <- ggarrange(
  p_dta, p_ph,
  nrow = 1,
  labels = c("A", "B"),
  font.label = list(size = 30)
)

ggsave(file.path(paths$figures, "supp_global_fdr.png"), supp_global_fdr_plot,
       width = 12, height = 6, units = "in", dpi = 150)

print(supp_global_fdr_plot)

Marginal effects of phosphorus deficiency on plant height and kernel traits

Boxplots pool CTRL and Inv4m; bracket labels come from the marginal emmeans contrast FDR (48-contrast D pool). PH has a robust marginal Treatment effect and non-significant Genotype × Treatment interaction. KW50, TKN, and TKW are shown together to address Sawers’ Figure 1 comment: under \(-P\) the plant trades kernel size for kernel number (KW50 drops significantly, TKN rises significantly in the pooled marginal contrast), leaving total kernel weight (TKW) statistically unchanged.

pal_trt     <- c("+P" = "gold", "-P" = "purple4")
supp_traits <- c("PH", "KW50", "TKN", "TKW")

marg_stats <- bracket_contrasts %>%
  filter(contrast_level == "marginal", trait %in% supp_traits)

supp_panel <- function(trait_name, ylab_text, y_limits = NULL) {
  d <- pheno %>%
    select(Treatment, value = all_of(trait_name)) %>%
    filter(!is.na(value))

  s <- marg_stats %>%
    filter(trait == trait_name) %>%
    mutate(y.position = max(d$value, na.rm = TRUE) +
                        0.08 * diff(range(d$value, na.rm = TRUE)))

  p <- ggplot(d, aes(x = Treatment, y = value, colour = Treatment)) +
    ggtitle(trait_name) +
    ylab(ylab_text) +
    geom_boxplot(width = 0.25, linewidth = 1, alpha = 0) %>%
      with_shadow(colour = "black", x_offset = 0, y_offset = 0, sigma = 1) +
    geom_quasirandom(size = 2) %>%
      with_shadow(colour = "black", x_offset = 0, y_offset = 0, sigma = 1) +
    scale_colour_manual(values = pal_trt) +
    theme_classic2(base_size = 20) +
    theme(
      plot.title      = element_markdown(hjust = 0.5, face = "bold"),
      axis.title.x    = element_blank(),
      axis.text.x     = element_text(size = 25, face = "bold", color = "black"),
      legend.position = "none"
    )

  if (!is.null(y_limits)) p <- p + coord_cartesian(ylim = y_limits)

  p <- p + scale_y_continuous(expand = expansion(mult = c(0.05, 0.15)))

  if (nrow(s) > 0)
    p <- p + stat_pvalue_manual(s, label = "p.adj.signif",
                                size = 10, bracket.size = 0.8, hide.ns = FALSE)
  p
}

p_PH   <- supp_panel("PH",   "Plant Height (cm)")
p_KW50 <- supp_panel("KW50", "50-kernel weight (g)")
p_TKN  <- supp_panel("TKN",  "Total Kernel Number")
p_TKW  <- supp_panel("TKW",  "Total Kernel Weight (g)")

supp_marginal <- ggarrange(
  p_PH, p_KW50, p_TKN, p_TKW,
  nrow = 1, ncol = 4, align = "h",
  labels = c("A", "B", "C", "D"),
  font.label = list(size = 30)
)

print(supp_marginal)

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] knitr_1.51       emmeans_2.0.1    rstatix_0.7.3    ggbeeswarm_0.7.3
##  [5] ggpubr_0.6.2     ggfx_1.0.3       ggtext_0.1.2     lubridate_1.9.4 
##  [9] forcats_1.0.1    stringr_1.6.0    dplyr_1.1.4      purrr_1.2.0     
## [13] readr_2.1.6      tidyr_1.3.2      tibble_3.3.0     ggplot2_4.0.1   
## [17] tidyverse_2.0.0  here_1.0.2      
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.1   vipor_0.4.7        farver_2.1.2       S7_0.2.1          
##  [5] fastmap_1.2.0      digest_0.6.39      timechange_0.3.0   estimability_1.5.1
##  [9] lifecycle_1.0.4    magrittr_2.0.4     compiler_4.5.2     rlang_1.1.6       
## [13] sass_0.4.10        tools_4.5.2        utf8_1.2.6         yaml_2.3.12       
## [17] ggsignif_0.6.4     labeling_0.4.3     xml2_1.5.1         RColorBrewer_1.1-3
## [21] abind_1.4-8        withr_3.0.2        grid_4.5.2         xtable_1.8-4      
## [25] scales_1.4.0       dichromat_2.0-0.1  cli_3.6.5          mvtnorm_1.3-3     
## [29] rmarkdown_2.30     ragg_1.5.0         generics_0.1.4     otel_0.2.0        
## [33] tzdb_0.5.0         commonmark_2.0.0   cachem_1.1.0       splines_4.5.2     
## [37] vctrs_0.6.5        Matrix_1.7-4       jsonlite_2.0.0     carData_3.0-5     
## [41] car_3.1-3          litedown_0.9       hms_1.1.4          Formula_1.2-5     
## [45] beeswarm_0.4.0     systemfonts_1.3.1  magick_2.9.0       jquerylib_0.1.4   
## [49] glue_1.8.0         cowplot_1.2.0      stringi_1.8.7      gtable_0.3.6      
## [53] pillar_1.11.1      htmltools_0.5.9    R6_2.6.1           textshaping_1.0.4 
## [57] rprojroot_2.1.1    evaluate_1.0.5     lattice_0.22-7     markdown_2.0      
## [61] backports_1.5.0    gridtext_0.1.5     broom_1.0.11       bslib_0.9.0       
## [65] Rcpp_1.1.0         coda_0.19-4.1      nlme_3.1-168       mgcv_1.9-4        
## [69] xfun_0.55          pkgconfig_2.0.3