PSU 2022 Inv4m Phenotype Analysis with Growth Curve Parameters

Overview

This analysis integrates logistic growth curve modeling with phenotype analysis:

1. Fit logistic growth curves to stover weight time-series data
2. Extract growth parameters (AUC, carrying capacity, inflection point)
3. Add parameters to phenotype dataset
4. Test effects with FDR correction across all traits
5. Calculate marginal effects for -P treatment using proper contrast directions
6. Visualize significant 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(growthcurver)
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

sp_corrected_file <- "PSU2022_spatially_corrected_both_treatments.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 39

---

2. Growth Curve Analysis

Prepare data for growthcurver

# Create time-series dataset for growth curves
dw_data <- sp_corrected_combined %>%
  filter(genotype %in% c("CTRL", "Inv4m")) %>%
  group_by(plotId) %>%
  mutate(
    # Replace harvest weight with maximum observed across timepoints
    STW121 = pmax(STW40, STW50, STW60, STWHV, na.rm = TRUE),
    # Add initial timepoint at 0
    STW01 = 0
  ) %>%
  ungroup() %>%
  select(
    plotId,
    Treatment, 
    Genotype = genotype,
    starts_with("STW")
  ) %>%
  pivot_longer(
    cols = c("STW01", "STW40", "STW50", "STW60", "STW121"),
    names_to = "DAP_string",
    values_to = "STW"
  ) %>%
  mutate(
    DAP = as.integer(gsub("STW", "", DAP_string, perl = TRUE))
  ) %>%
  filter(!is.na(STW)) %>%
  droplevels()

# Add 1 to avoid log(0) issues in growthcurver
dw_data <- dw_data %>%
  mutate(DW1 = STW + 1)

cat("Growth data prepared\n")

## Growth data prepared

cat("Total observations:", nrow(dw_data), "\n")

## Total observations: 302

cat("Unique plots:", n_distinct(dw_data$plotId), "\n")

## Unique plots: 64

Filter samples with sufficient data

# Identify samples with more than 3 timepoints
sample_counts <- dw_data %>%
  ungroup() %>%
  group_by(Genotype, Treatment, plotId) %>%
  summarise(n = sum(!is.na(DW1)), .groups = "drop")

cat("\nTimepoint distribution:\n")

## 
## Timepoint distribution:

print(table(sample_counts$n))

## 
##  3  4  5 
##  7  4 53

morethan4 <- sample_counts %>%
  filter(n > 3) %>%
  arrange(plotId) %>%
  pull(plotId)

for_curves <- paste0("plot_", morethan4)

cat("\nSamples with >3 timepoints:", length(for_curves), "\n")

## 
## Samples with >3 timepoints: 57

cat("By treatment:\n")

## By treatment:

sample_counts %>%
  filter(plotId %in% morethan4) %>%
  count(Treatment, Genotype) %>%
  print()

## # A tibble: 4 × 3
##   Treatment Genotype     n
##   <fct>     <fct>    <int>
## 1 +P        CTRL        13
## 2 +P        Inv4m       13
## 3 -P        CTRL        15
## 4 -P        Inv4m       16

Prepare wide format for growthcurver

# Create wide format required by growthcurver
without_na <- dw_data %>%
  ungroup() %>%
  mutate(plot_id = paste0("plot_", plotId)) %>%
  filter(plot_id %in% for_curves) %>%
  select(plot_id, time = DAP, DW1) %>%
  pivot_wider(names_from = "plot_id", values_from = DW1)

cat("\nWide format dimensions:", dim(without_na), "\n")

## 
## Wide format dimensions: 5 58

cat("Timepoints:", without_na$time, "\n")

## Timepoints: 1 40 50 60 121

Fit logistic growth curves

# Fit growth curves using growthcurver
growth_sum <- SummarizeGrowthByPlate(
  without_na, 
  bg_correct = "none"
)

cat("\nGrowth curves fitted\n")

## 
## Growth curves fitted

cat("Initial fits:", nrow(growth_sum), "samples\n")

## Initial fits: 57 samples

# Check for failed fits
cat("\nFit quality check:\n")

## 
## Fit quality check:

cat("Samples with k = 0:", sum(growth_sum$k == 0, na.rm = TRUE), "\n")

## Samples with k = 0: 4

cat("Samples with k = NA:", sum(is.na(growth_sum$k)), "\n")

## Samples with k = NA: 0

cat("Samples with sigma < 0:", sum(growth_sum$sigma < 0, na.rm = TRUE), "\n")

## Samples with sigma < 0: 0

Filter and validate growth parameters

# Remove failed fits
growth_sum_valid <- growth_sum %>%
  filter(
    k > 0,                    # Carrying capacity must be positive
    !is.na(k),                # No NAs
    !is.na(auc_l),            # Valid AUC
    !is.na(t_mid),            # Valid inflection point
    sigma > 0,                # Positive sigma (growth rate)
    note == ""                # No error notes from growthcurver
  )

cat("\nAfter quality filtering:\n")

## 
## After quality filtering:

cat("Valid fits:", nrow(growth_sum_valid), "samples\n")

## Valid fits: 53 samples

cat("Removed:", nrow(growth_sum) - nrow(growth_sum_valid), "samples\n")

## Removed: 4 samples

# Report by treatment
valid_by_treatment <- sp_corrected_combined %>%
  mutate(sample = paste0("plot_", plotId)) %>%
  filter(sample %in% growth_sum_valid$sample) %>%
  count(Treatment, genotype) %>%
  rename(Genotype = genotype, n_valid = n)

cat("\nValid fits by treatment:\n")

## 
## Valid fits by treatment:

print(valid_by_treatment)

##   Treatment Genotype n_valid
## 1        +P     CTRL      12
## 2        +P    Inv4m      12
## 3        -P     CTRL      13
## 4        -P    Inv4m      16

Extract and merge growth parameters

# Extract growth parameters and merge with metadata
growth_params <- sp_corrected_combined %>%
  select(
    plotId,
    Treatment, 
    Genotype = genotype
  ) %>%
  mutate(sample = paste0("plot_", plotId)) %>%
  inner_join(
    growth_sum_valid %>% 
      select(sample, auc_e, auc_l, k, t_mid, sigma, r),
    by = "sample"
  ) %>%
  rename(
    AUCE = auc_e,
    AUCL = auc_l,
    STWMAX = k,
    TMID = t_mid,
    SIGMA = sigma,
    R = r
  ) %>%
  mutate(
    AUCE = AUCE / 1000,  # Convert to kg*day
    AUCL = AUCL / 1000   # Convert to kg*day
  ) %>%
  select(plotId, AUCE, AUCL, STWMAX, TMID, SIGMA, R)

cat("\nGrowth parameters summary:\n")

## 
## Growth parameters summary:

summary(growth_params %>% select(-plotId))

##       AUCE            AUCL           STWMAX            TMID      
##  Min.   :3.931   Min.   :4.266   Min.   : 63.84   Min.   :51.06  
##  1st Qu.:4.654   1st Qu.:5.146   1st Qu.: 82.88   1st Qu.:54.18  
##  Median :5.039   Median :5.567   Median : 88.01   Median :56.44  
##  Mean   :5.492   Mean   :5.984   Mean   : 92.60   Mean   :56.44  
##  3rd Qu.:6.524   3rd Qu.:6.960   3rd Qu.:103.46   3rd Qu.:58.61  
##  Max.   :8.381   Max.   :9.497   Max.   :151.50   Max.   :66.96  
##      SIGMA               R          
##  Min.   : 0.6371   Min.   :0.06469  
##  1st Qu.: 1.0588   1st Qu.:0.12531  
##  Median : 1.9594   Median :0.15373  
##  Mean   : 2.6417   Mean   :0.17273  
##  3rd Qu.: 3.4550   3rd Qu.:0.19955  
##  Max.   :10.3253   Max.   :0.51982

cat("\nChecking for remaining zeros:\n")

## 
## Checking for remaining zeros:

cat("AUCE zeros:", sum(growth_params$AUCE == 0, na.rm = TRUE), "\n")

## AUCE zeros: 0

cat("AUCL zeros:", sum(growth_params$AUCL == 0, na.rm = TRUE), "\n")

## AUCL zeros: 0

cat("STWMAX zeros:", sum(growth_params$STWMAX == 0, na.rm = TRUE), "\n")

## STWMAX zeros: 0

cat("TMID zeros:", sum(growth_params$TMID == 0, na.rm = TRUE), "\n")

## TMID zeros: 0

---

3. Merge with Phenotype Data

Create complete dataset

# Merge growth parameters with phenotype data
pheno <- sp_corrected_combined %>%
  filter(genotype %in% c("CTRL", "Inv4m")) %>%
  mutate(genotype = factor(genotype, levels = c("CTRL", "Inv4m"))) %>%
  droplevels() %>%
  rename(Genotype = genotype) %>%
  left_join(growth_params, by = "plotId") %>%
  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("Complete dataset dimensions:", dim(pheno), "\n")

## Complete dataset dimensions: 64 46

cat("Growth parameters added successfully\n")

## Growth parameters added successfully

---

4. 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"
)
growth_phenotypes <- c("AUCE", "AUCL", "STWMAX", "TMID")

available_phenos <- intersect(
  c(base_phenotypes, ear_phenotypes, growth_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, AUCE, AUCL, STWMAX, TMID

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: 60

print(all_effects %>% head(20))

##    response                 predictor    Estimate  Std.Error      p.value
## 1     STW50               Treatment-P -18.7588852 1.58289202 1.598759e-16
## 2       DTS               Treatment-P   3.0621067 0.32199366 1.383332e-13
## 3       DTA               Treatment-P   3.4328993 0.36759849 2.668769e-13
## 4      AUCE               Treatment-P  -2.2565608 0.26520260 3.219010e-11
## 5      AUCL               Treatment-P  -2.0473606 0.28017238 2.218196e-09
## 6     STW60               Treatment-P -33.3434615 4.97340523 1.741324e-08
## 7     STW40               Treatment-P  -6.5170078 1.00758227 2.137127e-08
## 8    STWMAX               Treatment-P -25.4966201 4.72280069 1.946125e-06
## 9     STWHV               Treatment-P -22.4074887 4.70104209 1.299584e-05
## 10      DTA             GenotypeInv4m  -1.4777836 0.36759849 1.648877e-04
## 11      DTS             GenotypeInv4m  -1.2841867 0.32199366 1.833645e-04
## 12     TMID               Treatment-P   4.5546699 1.15058065 2.434637e-04
## 13       PH             GenotypeInv4m   5.5697807 1.48063541 3.853788e-04
## 14     KW50               Treatment-P  -1.8013514 0.53628018 1.521219e-03
## 15    STW40             GenotypeInv4m  -2.9345129 1.00758227 5.056342e-03
## 16       CD GenotypeInv4m:Treatment-P  -2.6176674 0.97676410 1.000036e-02
## 17       HI             GenotypeInv4m   0.1718492 0.06560131 1.168579e-02
## 18    STW50             GenotypeInv4m  -4.1128444 1.63844890 1.515412e-02
## 19    STWHV             GenotypeInv4m -11.6548319 4.70104209 1.609593e-02
## 20       EL             GenotypeInv4m   0.8936462 0.38817860 2.561592e-02
##        p.adjust
## 1  9.592553e-15
## 2  4.149997e-12
## 3  5.337538e-12
## 4  4.828515e-10
## 5  2.661835e-08
## 6  1.741324e-07
## 7  1.831823e-07
## 8  1.459594e-05
## 9  8.663896e-05
## 10 9.893262e-04
## 11 1.000170e-03
## 12 1.217319e-03
## 13 1.778671e-03
## 14 6.519511e-03
## 15 2.022537e-02
## 16 3.750135e-02
## 17 4.124396e-02
## 18 5.051374e-02
## 19 5.082926e-02
## 20 7.684777e-02

Identify significant effects

# 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.7588852 1.58289202 1.598759e-16
## 2       DTS               Treatment-P   3.0621067 0.32199366 1.383332e-13
## 3       DTA               Treatment-P   3.4328993 0.36759849 2.668769e-13
## 4      AUCE               Treatment-P  -2.2565608 0.26520260 3.219010e-11
## 5      AUCL               Treatment-P  -2.0473606 0.28017238 2.218196e-09
## 6     STW60               Treatment-P -33.3434615 4.97340523 1.741324e-08
## 7     STW40               Treatment-P  -6.5170078 1.00758227 2.137127e-08
## 8    STWMAX               Treatment-P -25.4966201 4.72280069 1.946125e-06
## 9     STWHV               Treatment-P -22.4074887 4.70104209 1.299584e-05
## 10      DTA             GenotypeInv4m  -1.4777836 0.36759849 1.648877e-04
## 11      DTS             GenotypeInv4m  -1.2841867 0.32199366 1.833645e-04
## 12     TMID               Treatment-P   4.5546699 1.15058065 2.434637e-04
## 13       PH             GenotypeInv4m   5.5697807 1.48063541 3.853788e-04
## 14     KW50               Treatment-P  -1.8013514 0.53628018 1.521219e-03
## 15    STW40             GenotypeInv4m  -2.9345129 1.00758227 5.056342e-03
## 16       CD GenotypeInv4m:Treatment-P  -2.6176674 0.97676410 1.000036e-02
## 17       HI             GenotypeInv4m   0.1718492 0.06560131 1.168579e-02
##        p.adjust
## 1  9.592553e-15
## 2  4.149997e-12
## 3  5.337538e-12
## 4  4.828515e-10
## 5  2.661835e-08
## 6  1.741324e-07
## 7  1.831823e-07
## 8  1.459594e-05
## 9  8.663896e-05
## 10 9.893262e-04
## 11 1.000170e-03
## 12 1.217319e-03
## 13 1.778671e-03
## 14 6.519511e-03
## 15 2.022537e-02
## 16 3.750135e-02
## 17 4.124396e-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, AUCE, AUCL, STW60, STW40, STWMAX, STWHV, TMID, PH, KW50, CD, HI

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

## Total: 14 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
    # Get emmeans for each factor
    emm_treatment <- emmeans(model, ~ Treatment | Genotype)
    emm_genotype <- emmeans(model, ~ Genotype | Treatment)
    
    # Treatment: -P minus +P (using revpairwise for correct direction)
    treatment_cond <- contrast(
      emm_treatment, 
      method = "revpairwise"
    ) %>% 
      summary() %>%
      as.data.frame()
    
    # Genotype: Inv4m minus CTRL (using revpairwise for correct direction)
    genotype_cond <- contrast(
      emm_genotype, 
      method = "revpairwise"
    ) %>% 
      summary() %>%
      as.data.frame()
    
    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"
        ) %>%
        select(response, predictor, contrast_type, estimate, SE, p.value)
    )
  } else {
    # Marginal contrasts for non-interaction traits
    # Get emmeans for each factor
    emm_treatment <- emmeans(model, ~ Treatment)
    emm_genotype <- emmeans(model, ~ Genotype)
    
    # Treatment: -P minus +P (using revpairwise for correct direction)
    treatment_marg <- contrast(
      emm_treatment, 
      method = "revpairwise"
    ) %>% 
      summary() %>%
      as.data.frame()
    
    # Genotype: Inv4m minus CTRL (using revpairwise for correct direction)
    genotype_marg <- contrast(
      emm_genotype, 
      method = "revpairwise"
    ) %>% 
      summary() %>%
      as.data.frame()
    
    bind_rows(
      treatment_marg %>%
        mutate(
          response = x,
          predictor = "Treatment",
          contrast_type = "marginal"
        ) %>%
        select(response, predictor, contrast_type, estimate, SE, p.value),
      genotype_marg %>%
        mutate(
          response = x,
          predictor = "Genotype",
          contrast_type = "marginal"
        ) %>%
        select(response, predictor, contrast_type, estimate, SE, p.value)
    )
  }
}) %>%
  bind_rows()

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

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

cat("Treatment contrasts: -P minus +P (negative = reduction in -P)\n")

## Treatment contrasts: -P minus +P (negative = reduction in -P)

cat("Genotype contrasts: Inv4m minus CTRL (positive = increase in Inv4m)\n")

## Genotype contrasts: Inv4m minus CTRL (positive = increase in Inv4m)

cat("Note: Using method='revpairwise' for correct contrast direction\n")

## Note: Using method='revpairwise' for correct contrast direction

print(contrasts_for_reporting)

##    response       predictor contrast_type     estimate         SE      p.value
## 1     STW50       Treatment      marginal -16.44895394 1.11112428 1.792984e-20
## 2     STW50        Genotype      marginal  -1.80291313 1.11112428 1.106085e-01
## 3       DTS       Treatment      marginal   3.41646064 0.22768390 5.655325e-22
## 4       DTS        Genotype      marginal  -0.92983272 0.22768390 1.331523e-04
## 5       DTA       Treatment      marginal   3.60161149 0.25993139 2.349827e-20
## 6       DTA        Genotype      marginal  -1.30907138 0.25993139 4.631183e-06
## 7      AUCE       Treatment      marginal  -1.95727725 0.18325865 2.163741e-14
## 8      AUCE        Genotype      marginal  -0.07283517 0.18325865 6.927657e-01
## 9      AUCL       Treatment      marginal  -1.80238777 0.19360297 2.047074e-12
## 10     AUCL        Genotype      marginal  -0.14116949 0.19360297 4.693688e-01
## 11    STW60       Treatment      marginal -27.94001566 3.47063576 1.390203e-10
## 12    STW60        Genotype      marginal   2.49794403 3.47063576 4.750381e-01
## 13    STW40       Treatment      marginal  -5.08603298 0.70669919 1.255415e-09
## 14    STW40        Genotype      marginal  -1.50353809 0.70669919 3.756391e-02
## 15   STWMAX       Treatment      marginal -23.00466068 3.26352021 5.569614e-09
## 16   STWMAX        Genotype      marginal  -3.88627456 3.26352021 2.394600e-01
## 17    STWHV       Treatment      marginal -18.69703663 3.26826688 3.915191e-07
## 18    STWHV        Genotype      marginal  -7.94437985 3.26826688 1.817984e-02
## 19     TMID       Treatment      marginal   3.49297279 0.79506705 5.968589e-05
## 20     TMID        Genotype      marginal  -1.30600849 0.79506705 1.068605e-01
## 21       PH       Treatment      marginal   3.42625576 1.04696734 1.770395e-03
## 22       PH        Genotype      marginal   6.40830361 1.04696734 7.712944e-08
## 23     KW50       Treatment      marginal  -1.86189594 0.37361666 8.188092e-06
## 24     KW50        Genotype      marginal   0.14235051 0.37361666 7.048452e-01
## 25       CD  Treatment|CTRL   conditional  -0.19062548 0.70101159 7.868190e-01
## 26       CD Treatment|Inv4m   conditional  -2.80829286 0.68018443 1.413928e-04
## 27       CD     Genotype|+P   conditional   0.15830684 0.71152496 8.248576e-01
## 28       CD     Genotype|-P   conditional  -2.45936054 0.66917885 5.888667e-04
## 29       HI       Treatment      marginal   0.08142121 0.04502794 7.670857e-02
## 30       HI        Genotype      marginal   0.12456063 0.04502794 7.974106e-03

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.44895394  -48.4821874 1.11112428
## 2     STW50        Genotype      marginal  -1.80291313   -6.7766437 1.11112428
## 3       DTS       Treatment      marginal   3.41646064    4.6792972 0.22768390
## 4       DTS        Genotype      marginal  -0.92983272   -1.2367197 0.22768390
## 5       DTA       Treatment      marginal   3.60161149    5.0456112 0.25993139
## 6       DTA        Genotype      marginal  -1.30907138   -1.7729346 0.25993139
## 7      AUCE       Treatment      marginal  -1.95727725  -29.8513487 0.18325865
## 8      AUCE        Genotype      marginal  -0.07283517   -1.2972632 0.18325865
## 9      AUCL       Treatment      marginal  -1.80238777  -25.8705328 0.19360297
## 10     AUCL        Genotype      marginal  -0.14116949   -2.3005474 0.19360297
## 11    STW60       Treatment      marginal -27.94001566  -36.9820542 3.47063576
## 12    STW60        Genotype      marginal   2.49794403    4.1403838 3.47063576
## 13    STW40       Treatment      marginal  -5.08603298  -47.9522487 0.70669919
## 14    STW40        Genotype      marginal  -1.50353809  -17.0561876 0.70669919
## 15   STWMAX       Treatment      marginal -23.00466068  -21.8614092 3.26352021
## 16   STWMAX        Genotype      marginal  -3.88627456   -4.0621516 3.26352021
## 17    STWHV       Treatment      marginal -18.69703663  -18.4823746 3.26826688
## 18    STWHV        Genotype      marginal  -7.94437985   -8.2939609 3.26826688
## 19     TMID       Treatment      marginal   3.49297279    6.3976309 0.79506705
## 20     TMID        Genotype      marginal  -1.30600849   -2.2913475 0.79506705
## 21       PH       Treatment      marginal   3.42625576    2.0835569 1.04696734
## 22       PH        Genotype      marginal   6.40830361    3.9326421 1.04696734
## 23     KW50       Treatment      marginal  -1.86189594  -17.6717047 0.37361666
## 24     KW50        Genotype      marginal   0.14235051    1.4930975 0.37361666
## 25       CD  Treatment|CTRL   conditional  -0.19062548   -0.7264286 0.70101159
## 26       CD Treatment|Inv4m   conditional  -2.80829286  -10.7017398 0.68018443
## 27       CD     Genotype|+P   conditional   0.15830684    0.6073076 0.71152496
## 28       CD     Genotype|-P   conditional  -2.45936054   -9.4693927 0.66917885
## 29       HI       Treatment      marginal   0.08142121   10.8866395 0.04502794
## 30       HI        Genotype      marginal   0.12456063   17.1493038 0.04502794
##         p.value
## 1  1.792984e-20
## 2  1.106085e-01
## 3  5.655325e-22
## 4  1.331523e-04
## 5  2.349827e-20
## 6  4.631183e-06
## 7  2.163741e-14
## 8  6.927657e-01
## 9  2.047074e-12
## 10 4.693688e-01
## 11 1.390203e-10
## 12 4.750381e-01
## 13 1.255415e-09
## 14 3.756391e-02
## 15 5.569614e-09
## 16 2.394600e-01
## 17 3.915191e-07
## 18 1.817984e-02
## 19 5.968589e-05
## 20 1.068605e-01
## 21 1.770395e-03
## 22 7.712944e-08
## 23 8.188092e-06
## 24 7.048452e-01
## 25 7.868190e-01
## 26 1.413928e-04
## 27 8.248576e-01
## 28 5.888667e-04
## 29 7.670857e-02
## 30 7.974106e-03

Create summary table for -P treatment effects

minus_p_effects <- contrasts_with_pct %>%
  filter(
    predictor == "Treatment" | grepl("Treatment\\|", predictor)
  ) %>%
  mutate(
    ci_lower = estimate - 1.96 * SE,
    ci_upper = estimate + 1.96 * SE,
    pct_ci_lower = (ci_lower / reference_mean) * 100,
    pct_ci_upper = (ci_upper / reference_mean) * 100
  ) %>%
  arrange(p.value) %>%
  select(
    Trait = response,
    Predictor = predictor,
    Type = contrast_type,
    Effect = estimate,
    `Effect %` = estimate_pct,
    SE,
    `95% CI Lower` = ci_lower,
    `95% CI Upper` = ci_upper,
    `P-value` = p.value
  )

cat("\n=== -P Treatment Effects Summary ===\n")

## 
## === -P Treatment Effects Summary ===

kable(
  minus_p_effects, 
  digits = 3,
  caption = "Marginal and Conditional Effects of -P Treatment"
)

Marginal and Conditional Effects of -P Treatment

Trait	Predictor	Type	Effect	Effect %	SE	95% CI Lower	95% CI Upper	P-value
DTS	Treatment	marginal	3.416	4.679	0.228	2.970	3.863	0.000
STW50	Treatment	marginal	-16.449	-48.482	1.111	-18.627	-14.271	0.000
DTA	Treatment	marginal	3.602	5.046	0.260	3.092	4.111	0.000
AUCE	Treatment	marginal	-1.957	-29.851	0.183	-2.316	-1.598	0.000
AUCL	Treatment	marginal	-1.802	-25.871	0.194	-2.182	-1.423	0.000
STW60	Treatment	marginal	-27.940	-36.982	3.471	-34.742	-21.138	0.000
STW40	Treatment	marginal	-5.086	-47.952	0.707	-6.471	-3.701	0.000
STWMAX	Treatment	marginal	-23.005	-21.861	3.264	-29.401	-16.608	0.000
STWHV	Treatment	marginal	-18.697	-18.482	3.268	-25.103	-12.291	0.000
KW50	Treatment	marginal	-1.862	-17.672	0.374	-2.594	-1.130	0.000
TMID	Treatment	marginal	3.493	6.398	0.795	1.935	5.051	0.000
CD	Treatment|Inv4m	conditional	-2.808	-10.702	0.680	-4.141	-1.475	0.000
PH	Treatment	marginal	3.426	2.084	1.047	1.374	5.478	0.002
HI	Treatment	marginal	0.081	10.887	0.045	-0.007	0.170	0.077
CD	Treatment|CTRL	conditional	-0.191	-0.726	0.701	-1.565	1.183	0.787

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, AUCE, AUCL, STW60, STW40, STWMAX, STWHV, TMID, PH, KW50, CD, HI

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

## Total: 14 traits

# Separate growth parameters from other traits
growth_traits_viz <- intersect(treatment_traits, growth_phenotypes)
other_traits_viz <- setdiff(treatment_traits, growth_phenotypes)

cat("\nGrowth parameters with effects:", length(growth_traits_viz), "\n")

## 
## Growth parameters with effects: 4

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

## AUCE, AUCL, STWMAX, TMID

cat("\nOther traits with effects:", length(other_traits_viz), "\n")

## 
## Other traits with effects: 10

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

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

---

5. Pairwise Comparisons for Visualization

Prepare time-series data

# Prepare stover dry weight time course 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()

cat("Time-series data prepared\n")

## Time-series data prepared

cat("Observations:", nrow(dw_raw), "\n")

## Observations: 236

Prepare t-tests for non-growth traits

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% other_traits_viz)

if (length(other_traits_viz) > 0) {
  treatment_pairwise <- to_t_test %>%
    group_by(Genotype, trait) %>%
    t_test(value ~ Treatment, ref.group = "+P") %>%
    adjust_pvalue(method = "fdr") %>%
    add_significance() %>%
    add_x_position(x = "Treatment") %>%
    add_y_position(scales = "free") %>%
    mutate(
      Genotype_label = factor(
        case_when(
          Genotype == "CTRL" ~ "CTRL",
          Genotype == "Inv4m" ~ "<i>Inv4m</i>",
          TRUE ~ as.character(Genotype)
        ),
        levels = c("CTRL", "<i>Inv4m</i>")
      )
    )
  
  cat("\n=== Pairwise t-test summary (non-growth traits) ===\n")
  print(
    treatment_pairwise %>% 
      select(trait, Genotype, statistic, p.adj, p.adj.signif), 
    n = 20
  )
} else {
  treatment_pairwise <- NULL
  cat("No non-growth traits with treatment effects\n")
}

## 
## === Pairwise t-test summary (non-growth traits) ===
## # A tibble: 20 × 5
##    trait Genotype statistic    p.adj p.adj.signif
##    <chr> <fct>        <dbl>    <dbl> <chr>       
##  1 CD    CTRL         0.187 8.55e- 1 ns          
##  2 DTA   CTRL        -9.54  7   e-10 ****        
##  3 DTS   CTRL       -10.4   2.79e-10 ****        
##  4 HI    CTRL        -1.70  1.24e- 1 ns          
##  5 KW50  CTRL         2.55  2.58e- 2 *           
##  6 PH    CTRL        -2.06  5.66e- 2 ns          
##  7 STW40 CTRL         5.93  8.67e- 6 ****        
##  8 STW50 CTRL        10.2   1.08e- 7 ****        
##  9 STW60 CTRL         7.18  2.40e- 6 ****        
## 10 STWHV CTRL         4.73  1.81e- 4 ***         
## 11 CD    Inv4m        6.92  9.86e- 7 ****        
## 12 DTA   Inv4m      -10.0   2.79e-10 ****        
## 13 DTS   Inv4m      -10.9   1.72e-10 ****        
## 14 HI    Inv4m       -0.585 5.95e- 1 ns          
## 15 KW50  Inv4m        5.56  1.75e- 5 ****        
## 16 PH    Inv4m       -2.54  2.19e- 2 *           
## 17 STW40 Inv4m        4.01  7.42e- 4 ***         
## 18 STW50 Inv4m        9.49  1.08e- 7 ****        
## 19 STW60 Inv4m        3.97  2.17e- 3 **          
## 20 STWHV Inv4m        3.25  5.27e- 3 **

Time-series t-tests for stover weight

treatment_timeseries_test <- dw_raw %>%
  ungroup() %>%
  filter(DAP > 1) %>%
  mutate(DAP = as.factor(DAP)) %>%
  group_by(Genotype, DAP) %>%
  t_test(STW ~ Treatment, ref.group = "+P") %>%
  adjust_pvalue(method = "fdr") %>%
  add_significance() %>%
  add_y_position(scales = "free_y") %>%
  add_x_position(x = "DAP", group = NULL, dodge = 1) %>%
  mutate(
    Genotype_label = factor(
      case_when(
        Genotype == "CTRL" ~ "CTRL",
        Genotype == "Inv4m" ~ "<i>Inv4m</i>",
        TRUE ~ as.character(Genotype)
      ),
      levels = c("CTRL", "<i>Inv4m</i>")
    )
  )

data_range <- range(dw_raw$STW[dw_raw$DAP > 1], na.rm = TRUE)
y_range <- diff(data_range)

treatment_timeseries_test <- treatment_timeseries_test %>%
  mutate(y.position = y.position + (0.08 * y_range))

cat("\n=== Time-series t-test summary ===\n")

## 
## === Time-series t-test summary ===

print(
  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
##   <fct>    <fct>     <dbl>        <dbl>       <dbl> <chr>       
## 1 CTRL     40         5.93 0.0000039    0.0000078   ****        
## 2 CTRL     50        10.2  0.0000000324 0.000000130 ****        
## 3 CTRL     60         7.18 0.000000961  0.00000256  ****        
## 4 CTRL     121        4.73 0.0000995    0.000159    ***         
## 5 Inv4m    40         4.01 0.000445     0.000593    ***         
## 6 Inv4m    50         9.49 0.0000000271 0.000000130 ****        
## 7 Inv4m    60         3.97 0.00141      0.00161     **          
## 8 Inv4m    121        3.25 0.00369      0.00369     **

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") +
  ylab("g") +
  xlab("Days After Planting") +
  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.12, 0.8),
    strip.text = element_markdown()
  )

print(p_dw)

---

6. 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)
    
    if (!is.null(t_test_data) && nrow(t_test_data) > 0) {
      if (!is.null(y_limits)) {
        y_range <- diff(y_limits)
        y_max <- y_limits[2]
        t_test_data <- t_test_data %>%
          mutate(y.position = y_max - (bracket_nudge * y_range))
      } else {
        data_range <- pheno %>%
          filter(!is.na(.data[[trait_name]])) %>%
          pull(.data[[trait_name]]) %>%
          range()
        y_range <- diff(data_range)
        y_max <- max(data_range)
        t_test_data <- t_test_data %>%
          mutate(y.position = y_max + (bracket_nudge * y_range))
      }
    }
  }
  
  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)
}

---

7. 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(68, 80),
    y_breaks = 68: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(68, 80),
    y_breaks = 68: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(5, 17),
    y_breaks = 5: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(21, 31),
    y_breaks = 21: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)
}

---

8. Growth Curve Parameter Visualization

Identify which growth parameters to plot

# Check if growth_traits_viz exists from section 4
# If not, create it now
if (!exists("growth_traits_viz")) {
  cat("Creating growth_traits_viz (not found from previous section)\n")
  
  # First ensure treatment_traits exists
  if (!exists("treatment_traits")) {
    treatment_traits <- contrasts_for_reporting %>%
      filter(grepl("Treatment", predictor)) %>%
      pull(response) %>%
      unique()
  }
  
  # Now create growth_traits_viz
  growth_traits_viz <- intersect(treatment_traits, growth_phenotypes)
  other_traits_viz <- setdiff(treatment_traits, growth_phenotypes)
}

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

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

cat("Available growth phenotypes:", paste(growth_phenotypes, collapse = ", "), "\n")

## Available growth phenotypes: AUCE, AUCL, STWMAX, TMID

cat("Traits with treatment effects:", paste(treatment_traits, collapse = ", "), "\n")

## Traits with treatment effects: STW50, DTS, DTA, AUCE, AUCL, STW60, STW40, STWMAX, STWHV, TMID, PH, KW50, CD, HI

cat("\ngrowth_traits_viz contains:", length(growth_traits_viz), "traits\n")

## 
## growth_traits_viz contains: 4 traits

if (length(growth_traits_viz) > 0) {
  cat("Trait names:", paste(growth_traits_viz, collapse = ", "), "\n")
} else {
  cat("No growth parameters have significant treatment effects\n")
}

## Trait names: AUCE, AUCL, STWMAX, TMID

Prepare t-tests for growth parameters

# Only run if there are growth parameters with significant effects
if (length(growth_traits_viz) > 0) {
  
  # Prepare long format data for the growth parameters we want to visualize
  growth_data_long <- pheno %>%
    select(plotId, Genotype, Treatment, all_of(growth_traits_viz)) %>%
    pivot_longer(
      cols = all_of(growth_traits_viz),
      names_to = "trait",
      values_to = "value"
    ) %>%
    filter(!is.na(value))  # Remove NA values
  
  cat("\nGrowth data prepared for t-tests:\n")
  cat("Total observations:", nrow(growth_data_long), "\n")
  growth_data_long %>%
    count(trait, Treatment, Genotype) %>%
    print()
  
  # Perform pairwise t-tests for each growth parameter
  growth_pairwise <- growth_data_long %>%
    group_by(Genotype, trait) %>%
    t_test(value ~ Treatment, ref.group = "+P") %>%
    adjust_pvalue(method = "fdr") %>%
    add_significance() %>%
    add_x_position(x = "Treatment") %>%
    add_y_position(scales = "free") %>%
    mutate(
      Genotype_label = factor(
        case_when(
          Genotype == "CTRL" ~ "CTRL",
          Genotype == "Inv4m" ~ "<i>Inv4m</i>",
          TRUE ~ as.character(Genotype)
        ),
        levels = c("CTRL", "<i>Inv4m</i>")
      ),
      y.position = y.position * 1.05
    )
  
  cat("\n=== Growth Parameter t-tests ===\n")
  print(
    growth_pairwise %>% 
      select(trait, Genotype, n1, n2, statistic, p, p.adj, p.adj.signif)
  )
  
} else {
  growth_pairwise <- NULL
  cat("\nNo growth parameters with significant treatment effects to test\n")
}

## 
## Growth data prepared for t-tests:
## Total observations: 212 
## # A tibble: 16 × 4
##    trait  Treatment Genotype     n
##    <chr>  <fct>     <fct>    <int>
##  1 AUCE   +P        CTRL        12
##  2 AUCE   +P        Inv4m       12
##  3 AUCE   -P        CTRL        13
##  4 AUCE   -P        Inv4m       16
##  5 AUCL   +P        CTRL        12
##  6 AUCL   +P        Inv4m       12
##  7 AUCL   -P        CTRL        13
##  8 AUCL   -P        Inv4m       16
##  9 STWMAX +P        CTRL        12
## 10 STWMAX +P        Inv4m       12
## 11 STWMAX -P        CTRL        13
## 12 STWMAX -P        Inv4m       16
## 13 TMID   +P        CTRL        12
## 14 TMID   +P        Inv4m       12
## 15 TMID   -P        CTRL        13
## 16 TMID   -P        Inv4m       16
## 
## === Growth Parameter t-tests ===
## # A tibble: 8 × 8
##   trait  Genotype    n1    n2 statistic           p      p.adj p.adj.signif
##   <chr>  <fct>    <int> <int>     <dbl>       <dbl>      <dbl> <chr>       
## 1 AUCE   CTRL        12    13      8.54 0.000000192 0.00000154 ****        
## 2 AUCL   CTRL        12    13      7.08 0.0000029   0.0000116  ****        
## 3 STWMAX CTRL        12    13      5.08 0.000139    0.000222   ***         
## 4 TMID   CTRL        12    13     -3.29 0.00352     0.00402    **          
## 5 AUCE   Inv4m       12    16      5.75 0.0000551   0.000147   ***         
## 6 AUCL   Inv4m       12    16      5.35 0.0000803   0.000161   ***         
## 7 STWMAX Inv4m       12    16      4.33 0.000525    0.0007     ***         
## 8 TMID   Inv4m       12    16     -2.72 0.0127      0.0127     *

Plot growth parameters

if (!is.null(growth_pairwise) && nrow(growth_pairwise) > 0) {
  
  # Function to create individual growth parameter plots
  plot_growth_param <- function(trait_name, 
                                 trait_title, 
                                 trait_ylab) {
    
    # Get t-test results for this trait
    t_test_data <- growth_pairwise %>%
      filter(trait == trait_name)
    
    # Create plot
    p <- pheno %>%
      filter(!is.na(.data[[trait_name]])) %>%
      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) +
      stat_pvalue_manual(
        t_test_data,
        size = 10,
        label = "p.adj.signif",
        bracket.size = 0.8,
        hide.ns = FALSE
      ) +
      scale_y_continuous(expand = expansion(mult = c(0.05, 0.1))) +
      pheno_theme2 +
      theme(strip.text = element_markdown())
    
    return(p)
  }
  
  # Create list to store plots
  growth_plots <- list()
  
  # Generate plot for each significant growth parameter
  if ("AUCE" %in% growth_traits_viz) {
    growth_plots$auce <- plot_growth_param(
      "AUCE", 
      "AUC<sub>empirical</sub>", 
      "kg × day"
    )
    cat("Created AUCE plot\n")
  }
  
  if ("AUCL" %in% growth_traits_viz) {
    growth_plots$aucl <- plot_growth_param(
      "AUCL", 
      "AUC<sub>logistic</sub>", 
      "kg × day"
    )
    cat("Created AUCL plot\n")
  }
  
  if ("STWMAX" %in% growth_traits_viz) {
    growth_plots$stwmax <- plot_growth_param(
      "STWMAX", 
      "STW<sub>max</sub>", 
      "g"
    )
    cat("Created STWMAX plot\n")
  }
  
  if ("TMID" %in% growth_traits_viz) {
    growth_plots$tmid <- plot_growth_param(
      "TMID", 
      "T<sub>1/2</sub>", 
      "Days"
    )
    cat("Created TMID plot\n")
  }
  
  # Arrange plots in a panel
  if (length(growth_plots) > 0) {
    p_gc <- ggarrange(
      plotlist = growth_plots,
      nrow = 1,
      ncol = length(growth_plots)
    )
    
    cat("\nCreated combined growth parameter panel with", 
        length(growth_plots), "plots\n")
    print(p_gc)
  }
  
} else {
  cat("\nNo growth parameter plots to generate\n")
  cat("This means none of the growth parameters (AUCE, AUCL, STWMAX, TMID)\n")
  cat("showed significant treatment effects in the linear models.\n")
}

## Created AUCE plot
## Created AUCL plot
## Created STWMAX plot
## Created TMID plot
## 
## Created combined growth parameter panel with 4 plots

---

9. Fitted Growth Curves Visualization

Prepare fitted curve data

# Generate predictions from fitted models - only for valid samples
valid_samples <- growth_sum_valid$sample

growth_fit <- lapply(valid_samples, function(x) {
  # Check if sample exists in data
  if (!x %in% colnames(without_na)) {
    return(NULL)
  }
  
  growth <- tryCatch(
    SummarizeGrowth(without_na$time, without_na[[x]]),
    error = function(e) {
      message("Failed to fit ", x)
      return(NULL)
    }
  )
  
  if (is.null(growth)) return(NULL)
  
  t <- 40:121
  dw <- predict(growth$model, newdata = data.frame(t = t))
  data.frame(sample = x, time = t, fittedDW = dw)
}) %>%
  bind_rows()

cat("Fitted curves generated for", length(unique(growth_fit$sample)), "samples\n")

## Fitted curves generated for 53 samples

# Merge with metadata
gc_fitted <- sp_corrected_combined %>%
  select(plotId, Treatment, Genotype = genotype) %>%
  mutate(
    sample = paste0("plot_", plotId),
    Genotype = factor(Genotype, levels = c("CTRL", "Inv4m")),
    Genotype_label = factor(
      case_when(
        Genotype == "CTRL" ~ "CTRL",
        Genotype == "Inv4m" ~ "<i>Inv4m</i>",
        TRUE ~ as.character(Genotype)
      ),
      levels = c("CTRL", "<i>Inv4m</i>")
    )
  ) %>%
  right_join(growth_fit, by = "sample") %>%
  filter(!is.na(Treatment)) %>%
  arrange(Treatment, sample) %>%
  mutate(
    plot_order = row_number(),
    sample = forcats::fct_reorder(
      as.factor(sample), 
      plot_order, 
      .fun = first
    )
  )

cat("Fitted curves with metadata:", nrow(gc_fitted), "observations\n")

## Fitted curves with metadata: 4346 observations

cat("By treatment:\n")

## By treatment:

gc_fitted %>%
  group_by(Treatment, Genotype) %>%
  summarise(n_samples = n_distinct(sample), .groups = "drop") %>%
  print()

## # A tibble: 4 × 3
##   Treatment Genotype n_samples
##   <fct>     <fct>        <int>
## 1 +P        CTRL            12
## 2 +P        Inv4m           12
## 3 -P        CTRL            13
## 4 -P        Inv4m           16

Plot fitted growth curves

p_fitted <- gc_fitted %>%
  ggplot(aes(x = time, y = fittedDW, color = Treatment, group = sample)) +
  ggtitle("Fitted Stover Weight") +
  ylab("g") +
  xlab("Days After Planting") +
  xlim(c(40, 130)) +
  ylim(c(0, 155)) +
  geom_line(linewidth = 0.8) %>% 
    with_shadow(colour = "black", x_offset = 0, y_offset = 0, sigma = 1) +
  scale_color_manual(values = pal) +
  labs(color = "") +
  facet_wrap(~ Genotype_label) +
  theme_classic2(base_size = 20) +
  theme(
    legend.position = c(0.8, 0.9),
    plot.title = element_text(hjust = 0.5, face = "bold", size = 25),
    strip.text = element_markdown(size = 20),
    legend.text = element_text(face = "bold", size = 20),
    strip.background = element_blank()
  )

print(p_fitted)

---

10. Combined Figures

Main phenotype panel with growth parameters

# Combine time series and growth parameters
if (exists("p_dw") && exists("growth_plots") && length(growth_plots) > 0) {
  # Create growth parameter panel
  p_gc <- ggarrange(
    plotlist = growth_plots,
    nrow = 1,
    ncol = length(growth_plots)
  )
  
  # Check if fitted curves exist
  if (exists("p_fitted")) {
    out_plot <- ggarrange(
      p_dw,
      p_fitted,
      p_gc,
      nrow = 3,
      ncol = 1,
      heights = c(1, 1, 1),
      labels = c("A", "B", "C"),
      font.label = list(size = 30)
    )
  } else {
    out_plot <- ggarrange(
      p_dw,
      p_gc,
      nrow = 2,
      ncol = 1,
      heights = c(1, 1),
      labels = c("A", "B"),
      font.label = list(size = 30)
    )
  }
  
  print(out_plot)
} else if (exists("p_dw")) {
  print(p_dw)
  cat("\nNote: Growth parameter plots not available\n")
}

---

11. Session Info

sessionInfo()

## R version 4.5.1 (2025-06-13)
## Platform: x86_64-apple-darwin20
## Running under: macOS Sequoia 15.6.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.5-x86_64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.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.50         growthcurver_0.3.1 emmeans_1.11.2     rstatix_0.7.2     
##  [5] ggbeeswarm_0.7.2   ggpubr_0.6.1       ggfx_1.0.2         ggtext_0.1.2      
##  [9] lubridate_1.9.4    forcats_1.0.0      stringr_1.5.1      dplyr_1.1.4       
## [13] purrr_1.1.0        readr_2.1.5        tidyr_1.3.1        tibble_3.3.0      
## [17] ggplot2_3.5.2      tidyverse_2.0.0   
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6       beeswarm_0.4.0     xfun_0.53          bslib_0.9.0       
##  [5] tzdb_0.5.0         vctrs_0.6.5        tools_4.5.1        generics_0.1.4    
##  [9] pkgconfig_2.0.3    RColorBrewer_1.1-3 lifecycle_1.0.4    compiler_4.5.1    
## [13] farver_2.1.2       textshaping_1.0.1  minpack.lm_1.2-4   carData_3.0-5     
## [17] litedown_0.7       vipor_0.4.7        htmltools_0.5.8.1  sass_0.4.10       
## [21] yaml_2.3.10        Formula_1.2-5      pillar_1.11.0      car_3.1-3         
## [25] jquerylib_0.1.4    cachem_1.1.0       magick_2.8.7       abind_1.4-8       
## [29] commonmark_2.0.0   tidyselect_1.2.1   digest_0.6.37      mvtnorm_1.3-3     
## [33] stringi_1.8.7      labeling_0.4.3     cowplot_1.2.0      fastmap_1.2.0     
## [37] grid_4.5.1         cli_3.6.5          magrittr_2.0.3     utf8_1.2.6        
## [41] broom_1.0.9        withr_3.0.2        scales_1.4.0       backports_1.5.0   
## [45] timechange_0.3.0   estimability_1.5.1 rmarkdown_2.29     ggsignif_0.6.4    
## [49] ragg_1.4.0         hms_1.1.3          evaluate_1.0.4     markdown_2.0      
## [53] rlang_1.1.6        gridtext_0.1.5     Rcpp_1.1.0         xtable_1.8-4      
## [57] glue_1.8.0         xml2_1.4.0         rstudioapi_0.17.1  jsonlite_2.0.0    
## [61] R6_2.6.1           systemfonts_1.2.3