1 Purpose

This R Markdown workflow imports raw OD values from one or more Tecan Infinite 200Pro 96-well plate Excel files and a sample metadata CSV file, then converts raw absorbance values into final cortisol concentrations in ng/mL.

When more than one plate is analyzed, each plate is processed independently through blank subtraction, B0 calculation, %B/B0 calculation, 4-parameter standard curve fitting, and concentration back-calculation. Only the final corrected sample concentrations are combined across plates for downstream plots and statistical comparisons.

The workflow does the following:

  1. Imports the raw 96-well plate OD values from one or more Tecan Excel exports.
  2. Imports the sample metadata file containing sample IDs and duplicate well locations.
  3. Finds the blank wells from the metadata row where SampleID == "BLK".
  4. Calculates the average raw OD of the blanks separately for each plate.
  5. Subtracts the plate-specific blank average from wells on that same plate.
  6. Averages duplicate wells for standards and samples.
  7. Calculates %B/B0 values using the plate-specific standard S0/B0 as the maximum binding reference.
  8. Fits a separate 4-parameter logistic standard curve for each plate using standards S1-S7.
  9. Back-calculates sample cortisol concentrations from that sample’s plate-specific standard curve.
  10. Applies dilution and/or extraction correction factors if needed.
  11. Outputs an updated metadata CSV with calculated values.
  12. Combines final corrected sample concentrations across plates for CSV output, plots, and statistics.
  13. Outputs a PDF with jittered boxplots for grouping variables in the metadata file.

2 Required packages

# Install these once if needed:
# install.packages(c("readxl", "readr", "dplyr", "tidyr", "stringr", "purrr", "ggplot2", "minpack.lm", "scales", "multcompView"))

library(readxl)
library(readr)
library(dplyr)
library(tidyr)
library(stringr)
library(purrr)
library(ggplot2)
library(minpack.lm)
library(scales)
library(multcompView)

3 User inputs

Edit the file paths below for each new ELISA plate.

# Raw Tecan Excel export file(s) containing the 96-well plate OD values.
#
# For a single plate, use one named entry:
# plate_excel_files <- c("1" = "Plate1_Tecan_RawData.xlsx")
#
# For multiple plates, add one named entry per plate. The names must match the
# values in the metadata column titled `Plate Number`.
# Example:
# plate_excel_files <- c(
#   "1" = "Plate1_Tecan_RawData.xlsx",
#   "2" = "Plate2_Tecan_RawData.xlsx"
# )
plate_excel_files <- c("1" = "Plate1_Tecan_RawData.xlsx")


# Sample metadata CSV.
# Required columns:
#   SampleID
#   Cell Position 1
#   Cell Position 2
#   Plate Number
# Additional columns such as Treatment, SampleDay, Time, etc. are allowed.
#
# Important: Cell positions are interpreted within each plate. For example, A3
# on Plate 1 and A3 on Plate 2 are treated as different wells because they have
# different Plate Number values in the metadata.
metadata_csv_file <- "ThermalPriming_Cortisol_Metadata_9July2026-.csv"

# Output files.
processed_metadata_output <- "processed_cortisol_metadata.csv"
sample_processing_diagnostics_output <- "sample_processing_diagnostics.csv"
samples_excluded_from_statistics_output <- "samples_excluded_from_statistics_diagnostic.csv"
group_representation_diagnostics_output <- "group_representation_diagnostics.csv"
standard_curve_plot_output <- "cortisol_standard_curve.pdf"
standard_curve_with_samples_plot_output <- "cortisol_standard_curve_with_samples.pdf"
competitive_elisa_interpretation_plot_output <- "cortisol_competitive_elisa_interpretation_with_samples.pdf"
grouped_boxplots_output <- "cortisol_grouped_boxplots.pdf"
statistical_summary_output <- "cortisol_group_statistics.csv"
anova_kruskal_summary_output <- "cortisol_overall_group_tests.csv"
posthoc_tukey_output <- "cortisol_tukey_posthoc_tests.csv"
posthoc_wilcoxon_output <- "cortisol_wilcoxon_posthoc_tests.csv"

# Quick-reference files that include only statistically significant post-hoc comparisons.
significant_tukey_output <- "cortisol_significant_tukey_posthoc_tests.csv"
significant_wilcoxon_output <- "cortisol_significant_wilcoxon_posthoc_tests.csv"
excluded_well_notes_output <- "cortisol_excluded_well_notes.csv"
significance_letters_output <- "cortisol_significance_letters.csv"
plate_processing_summary_output <- "cortisol_plate_processing_summary.csv"

# Outlier-screening outputs. These do not alter the primary analysis.
outlier_screening_output <- "cortisol_extreme_outlier_screening.csv"
grouped_boxplots_extreme_outliers_removed_output <- "cortisol_grouped_boxplots_extreme_outliers_removed.pdf"
statistical_summary_extreme_outliers_removed_output <- "cortisol_group_statistics_extreme_outliers_removed.csv"
anova_kruskal_summary_extreme_outliers_removed_output <- "cortisol_overall_group_tests_extreme_outliers_removed.csv"
posthoc_tukey_extreme_outliers_removed_output <- "cortisol_tukey_posthoc_tests_extreme_outliers_removed.csv"
posthoc_wilcoxon_extreme_outliers_removed_output <- "cortisol_wilcoxon_posthoc_tests_extreme_outliers_removed.csv"
significant_tukey_extreme_outliers_removed_output <- "cortisol_significant_tukey_posthoc_tests_extreme_outliers_removed.csv"
significant_wilcoxon_extreme_outliers_removed_output <- "cortisol_significant_wilcoxon_posthoc_tests_extreme_outliers_removed.csv"

# Significance threshold used for the quick-reference post-hoc outputs and plot letters.
posthoc_significance_alpha <- 0.05

# Add compact letter displays to final jittered boxplots.
# Groups that share a letter are not significantly different.
# Groups with no letters in common are significantly different at the alpha level above.
show_significance_letters <- TRUE

# Choose which post-hoc method should be used for the plot letters.
# Options:
#   "tukey"    = ANOVA/Tukey-based letters; best when ANOVA assumptions are reasonable.
#   "wilcoxon" = Pairwise Wilcoxon/BH-adjusted letters; useful when emphasizing non-parametric results.
significance_letter_method <- "wilcoxon"

# Concentration back-calculation/extrapolation options.
#
# ELISA concentrations are most reliable when sample %B/B0 values fall within the
# measured standard curve range, and especially within the common 20-80% B/B0
# working range. For preliminary work, you may still want to calculate values
# outside the standard range using the fitted 4PL equation.
#
# TRUE  = keep finite 4PL back-calculated values even if they are outside the
#         observed standard curve range, and flag them as extrapolated/cautious.
# FALSE = set values outside the observed standard curve range to NA.
allow_curve_extrapolation <- TRUE

# Include extrapolated values in downstream plots/statistics?
# TRUE is useful for exploratory/preliminary analyses, but final/publication
# analyses should usually inspect these flags carefully.
include_extrapolated_concentrations_in_statistics <- TRUE

# Quantification-confidence labels in the output table and standard-curve-with-samples plot.
# These labels do not alter the concentration calculation; they are interpretive flags.
# Preferred assay range is commonly 20-80% B/B0. Values within the measured
# standard range but outside 20-80% are still interpolation values, but are less ideal.
preferred_percent_B_B0_lower <- 20
preferred_percent_B_B0_upper <- 80

# When TRUE, the sample-on-standard-curve plot will show the measured standard
# concentration range as a lightly shaded vertical band. Points outside this band
# are extrapolated or tail-estimated and should be interpreted cautiously.
show_standard_range_shading <- TRUE

# Add a teaching/interpretation plot showing the inverse competitive ELISA
# relationship, overlaid with actual samples. This is separate from the formal
# standard-curve plot and is meant to help interpret dilution/rerun decisions.
make_competitive_elisa_interpretation_plot <- TRUE

# Label the most extreme actual sample points on the interpretation plot.
label_extreme_samples_on_interpretation_plot <- TRUE
number_of_extreme_samples_to_label <- 10

# If a sample %B/B0 value is outside the observed standard range but still inside
# the fitted 4PL curve/asymptotes, the script will use the 4PL equation to
# extrapolate a concentration. This preserves variation among samples instead of
# assigning all values to the nearest standard concentration.
#
# If a sample is outside the fitted 4PL asymptotes, the algebraic 4PL inverse has
# no exact real-valued solution. For exploratory/preliminary work, the script can
# optionally use a local log-linear tail extrapolation based on the two nearest
# standards at that end of the curve. This is less defensible than true 4PL
# interpolation/extrapolation, but it avoids artificial identical boundary values.
use_tail_extrapolation_outside_model_asymptotes <- TRUE

# Enforce the expected competitive ELISA direction during tail extrapolation.
# For cortisol competitive ELISA: higher %B/B0 must produce lower concentration,
# even when values are extrapolated beyond the fitted curve.
enforce_competitive_monotonic_tail_extrapolation <- TRUE

# Boundary estimates are now OFF by default because they collapse all out-of-range
# samples to the same minimum/maximum standard value and can visually erase real
# variation. Turn this on only if you specifically want censored values such as
# "at or below lowest standard" or "at or above highest standard."
use_boundary_estimates_outside_model_asymptotes <- FALSE

# Boundary values used only when use_boundary_estimates_outside_model_asymptotes = TRUE.
# These are retained as an optional fallback but are not used by default.
lower_boundary_concentration_ng_mL <- 0.005
upper_boundary_concentration_ng_mL <- 50

# Optional manually excluded wells.
# Add any wells that should be excluded because of a visible pipetting issue, bubble, plate reader artifact,
# or any other concern with the raw OD value.
#
# For one plate, you can list wells directly:
# excluded_wells <- c("C7", "D7")
#
# For multiple plates, use plate-specific labels written as "PlateNumber:Well":
# excluded_wells <- c("1:C7", "2:D7")
#
# Direct well labels without a plate number are applied to ALL plates, so plate-specific
# labels are safer when analyzing more than one plate.
# Leave empty if no wells should be excluded.
excluded_wells <- c()

# Plot options for the grouped concentration plots.
# Set use_log10_y_axis <- TRUE if cortisol concentrations span a wide range.
use_log10_y_axis <- FALSE

# Metadata-driven grouping options.
# The first three metadata columns are always treated as identifiers/well positions:
#   SampleID, Cell Position 1, Cell Position 2
# Every additional metadata column is treated as a possible grouping variable unless excluded below.
#
# Leave grouping_variables_to_analyze <- NULL to analyze all eligible metadata grouping columns.
# Or provide a character vector to analyze only selected metadata columns, for example:
# grouping_variables_to_analyze <- c("Treatment", "SampleDay")
grouping_variables_to_analyze <- NULL

# Column used to count unique research samples in statistical summaries.
# Default = SampleID. If your metadata has a true biological individual ID
# column, such as OysterID, AnimalID, IndividualID, etc., change this to that
# column name so n reflects unique biological individuals rather than rows.
#
# Important: statistical tests still use one concentration value per metadata row.
# The unique-n column is mainly for accurate reporting and diagnostics.
statistics_unique_id_column <- "SampleID"

# Optional column identifying the biological individual/organism.
# For your oyster dataset this can be "Oyster". If the column is present,
# diagnostics will report both tissue/sample-level n and oyster-level n.
# If absent, the script will simply skip oyster-level diagnostic counts.
biological_individual_id_column <- "Oyster"

# Optional expected design checks. Set to NA to skip a check.
# These are only diagnostic warnings; they do not change the calculations.
expected_individuals_per_simple_group <- 8
expected_min_tissues_per_individual <- 4
expected_max_tissues_per_individual <- 5

# Exclude standards from plots/statistics. Standards are needed for the curve,
# but they are not research samples and should not contribute to group n values.
exclude_standards_from_plots_and_statistics <- TRUE

# Optional: manually exclude additional SampleID values from plots/statistics.
# This is useful for QC decisions after processing. Leave empty if not needed.
exclude_sample_ids_from_plots_and_statistics <- c()

# Optional extreme outlier screening.
# This does not change the primary figures/statistics. If enabled, the report will
# print and save a list of samples that appear to be extreme outliers, then create
# a second exploratory set of figures/statistics with those rows excluded.
#
# Method: within each grouping variable, extreme outliers are flagged using the
# conservative boxplot rule: values below Q1 - 3*IQR or above Q3 + 3*IQR.
# This is stricter than the usual 1.5*IQR rule and is intended only as a QC screen.
screen_extreme_outliers <- TRUE
extreme_outlier_iqr_multiplier <- 3
produce_secondary_outputs_without_extreme_outliers <- TRUE

# Minimum number of observations required within a group before outliers are screened.
# Groups smaller than this are skipped because quartile/IQR outlier calls are unstable.
minimum_group_n_for_outlier_screening <- 6

# If TRUE, a row flagged as an extreme outlier for ANY grouping variable will be
# excluded from the secondary exploratory outputs. The primary outputs still include it.
exclude_if_extreme_outlier_in_any_grouping <- TRUE

# Optional combined grouping variables.
# This allows comparisons such as TreatmentDay, DoseTime, SiteTreatment, etc., without hard-coding
# any specific metadata structure. Each named item becomes a new grouping variable by pasting together
# the listed metadata columns.
#
# Examples:
# combined_grouping_variables <- list(TreatmentDay = c("Treatment", "SampleDay"))
# combined_grouping_variables <- list(DoseTime = c("Dose", "Time"), SiteTreatment = c("Site", "Treatment"))
# Leave as list() if no combined grouping variable is needed.
combined_grouping_variables <- list(Ploidy = c("Ploidy", "Tissue"))

# Optional focused plotting layout.
# This creates an additional plot with one metadata variable on the x-axis and another shown by color.
# Leave either value as NULL to skip this extra plot.
# This is useful for designs like x = SampleDay and color = Treatment, but can be changed
# to match any metadata columns in a future CSV.
focused_plot_x_variable <- "Treatment"
focused_plot_color_variable <- "Tissue"

# Set run_group_statistics <- TRUE to output statistical comparisons for each grouping variable.
# For grouping variables with 2 groups, this uses a Welch two-sample t-test.
# For grouping variables with 3 or more groups, this uses a one-way ANOVA.
run_group_statistics <- TRUE

# Optional correction factors.
# If every sample has the same dilution factor, set it here.
# Example: if samples were diluted 1:200 before loading, use 200.
default_sample_dilution_factor <- 50

# If every sample has the same extraction dilution/concentration correction factor, set it here.
# Leave as 1 if no extraction correction is needed.
default_extraction_correction_factor <- 1

# If your metadata file contains columns named DilutionFactor and/or ExtractionCorrectionFactor,
# those sample-specific values will override the defaults above.

4 Worked example: dilution and extraction corrections

The script supports two separate correction factors so the output remains transparent:

  1. DilutionFactor corrects for any dilution made before loading the ELISA plate.
  2. ExtractionCorrectionFactor corrects for any dilution, concentration, or scaling that occurred during sample extraction/preparation before the ELISA dilution step.

The final corrected concentration is calculated as:

FinalConcentration_ng_mL <- Concentration_ng_mL * DilutionFactor * ExtractionCorrectionFactor

4.1 Example

Suppose a tissue extract is diluted 1:4 before being loaded onto the ELISA plate. The ELISA standard curve estimates the diluted sample concentration as 0.75 ng/mL.

If no additional extraction correction is needed:

measured_concentration <- 0.75
DilutionFactor <- 4
ExtractionCorrectionFactor <- 1

final_concentration <- measured_concentration * DilutionFactor * ExtractionCorrectionFactor
final_concentration
[1] 3

So the corrected concentration would be 3.0 ng/mL.

If the extraction process also requires a 2x correction, then:

measured_concentration <- 0.75
DilutionFactor <- 4
ExtractionCorrectionFactor <- 2

final_concentration <- measured_concentration * DilutionFactor * ExtractionCorrectionFactor
final_concentration
[1] 6

So the corrected concentration would be 6.0 ng/mL.

4.2 Per-sample correction factors

Yes: this R Markdown is designed to use per-sample correction factors when they are present in the metadata file.

If your metadata CSV contains a column named DilutionFactor, those sample-specific values are used. If that column is absent, the script uses default_sample_dilution_factor for every sample.

If your metadata CSV contains a column named ExtractionCorrectionFactor, those sample-specific values are used. If that column is absent, the script uses default_extraction_correction_factor for every sample.

For example, your metadata CSV can include:

SampleID Cell Position 1 Cell Position 2 Treatment DilutionFactor ExtractionCorrectionFactor
Sample_001 A3 A4 Control 4 1
Sample_002 B3 B4 HeatStress 10 2
Sample_003 C3 C4 HeatStress 1 1

This allows each sample to be corrected individually when needed.

5 Helper functions

# Convert a well location such as "A3" or "H12" into row and column components.
parse_well_position <- function(well_position) {
  well_position <- str_trim(as.character(well_position))
  well_position <- str_to_upper(well_position)
  
  tibble(
    Well = well_position,
    PlateRow = str_extract(well_position, "^[A-H]"),
    PlateCol = as.integer(str_extract(well_position, "[0-9]+$"))
  )
}

# Normalize plate numbers so metadata values and plate_excel_files names match reliably.
normalize_plate_number <- function(plate_number) {
  str_trim(as.character(plate_number))
}

# Read one Tecan Excel export and return the 96-well OD table in long format.
read_tecan_plate <- function(tecan_excel_file, plate_number) {
  tecan_raw <- read_excel(tecan_excel_file, col_names = FALSE)

  plate_header_row <- which(tecan_raw[[1]] == "<>")

  if (length(plate_header_row) != 1) {
    stop("Could not uniquely identify the plate header row marked with '<>' in file: ", tecan_excel_file)
  }

  plate_matrix <- read_excel(
    tecan_excel_file,
    skip = plate_header_row - 1,
    n_max = 9,
    col_names = TRUE
  )

  colnames(plate_matrix)[1] <- "PlateRow"

  plate_matrix %>%
    pivot_longer(
      cols = -PlateRow,
      names_to = "PlateCol",
      values_to = "RawOD"
    ) %>%
    mutate(
      `Plate Number` = normalize_plate_number(plate_number),
      PlateRow = str_to_upper(as.character(PlateRow)),
      PlateCol = as.integer(PlateCol),
      Well = paste0(PlateRow, PlateCol),
      PlateWellID = paste(`Plate Number`, Well, sep = ":"),
      RawOD = as.numeric(RawOD)
    ) %>%
    select(`Plate Number`, PlateWellID, Well, PlateRow, PlateCol, RawOD)
}

# Extract an OD value from the long-format plate table using a plate number and well location like "A3".
get_od_from_well <- function(well_position, plate_number, plate_long, value_column = "RawOD") {
  parsed <- parse_well_position(well_position)
  plate_number <- normalize_plate_number(plate_number)

  value <- plate_long %>%
    filter(`Plate Number` == plate_number, PlateRow == parsed$PlateRow, PlateCol == parsed$PlateCol) %>%
    pull({{ value_column }})

  if (length(value) != 1) {
    stop("Could not uniquely identify well position ", well_position, " on plate ", plate_number, ".")
  }

  value
}

# Convert manually entered excluded wells into a plate-aware table.
# Entries like "1:C7" exclude C7 only on Plate 1.
# Entries like "C7" exclude C7 on every plate.
make_excluded_well_table <- function(excluded_wells, plate_numbers) {
  excluded_wells <- str_to_upper(str_trim(as.character(excluded_wells)))
  excluded_wells <- excluded_wells[!is.na(excluded_wells) & excluded_wells != ""]
  plate_numbers <- normalize_plate_number(plate_numbers)

  if (length(excluded_wells) == 0) {
    return(tibble(`Plate Number` = character(), Well = character(), PlateWellID = character()))
  }

  map_dfr(excluded_wells, function(entry) {
    if (str_detect(entry, ":")) {
      parts <- str_split_fixed(entry, ":", 2)
      tibble(
        `Plate Number` = normalize_plate_number(parts[, 1]),
        Well = str_to_upper(str_trim(parts[, 2]))
      )
    } else {
      tibble(
        `Plate Number` = plate_numbers,
        Well = entry
      )
    }
  }) %>%
    distinct() %>%
    mutate(PlateWellID = paste(`Plate Number`, Well, sep = ":"))
}

# Summarize a duplicate pair while respecting manually excluded wells.
# If both wells are usable, the two values are averaged.
# If one well is excluded, the usable single value is carried forward and a note is created.
# If both wells are excluded, the result is set to NA and a warning note is created.
summarize_duplicate_wells <- function(well1, well2, plate_number, plate_long, excluded_well_table = tibble()) {
  plate_number <- normalize_plate_number(plate_number)
  wells <- str_to_upper(str_trim(c(as.character(well1), as.character(well2))))
  plate_well_ids <- paste(plate_number, wells, sep = ":")

  raw_values <- map_dbl(wells, get_od_from_well, plate_number = plate_number, plate_long = plate_long, value_column = "RawOD")
  corrected_values <- map_dbl(wells, get_od_from_well, plate_number = plate_number, plate_long = plate_long, value_column = "BlankCorrectedOD")
  excluded_flags <- plate_well_ids %in% excluded_well_table$PlateWellID
  usable_flags <- !excluded_flags

  n_excluded <- sum(excluded_flags)
  n_used <- sum(usable_flags)
  excluded_list <- paste(plate_well_ids[excluded_flags], collapse = ";")
  used_list <- paste(plate_well_ids[usable_flags], collapse = ";")

  note <- dplyr::case_when(
    n_excluded == 0 ~ NA_character_,
    n_excluded == 1 ~ paste0("One duplicate well excluded (", excluded_list, "); value calculated from single usable well: ", used_list, "."),
    n_excluded == 2 ~ paste0("Both duplicate wells excluded (", excluded_list, "); average and concentration set to NA."),
    TRUE ~ NA_character_
  )

  tibble(
    RawOD_1 = raw_values[1],
    RawOD_2 = raw_values[2],
    BlankCorrectedOD_1 = corrected_values[1],
    BlankCorrectedOD_2 = corrected_values[2],
    Well1_Excluded = excluded_flags[1],
    Well2_Excluded = excluded_flags[2],
    ExcludedWellsForSample = if_else(excluded_list == "", NA_character_, excluded_list),
    WellsUsedForAverage = if_else(used_list == "", NA_character_, used_list),
    NumberOfWellsUsedForAverage = n_used,
    AverageRawOD = if_else(n_used > 0, mean(raw_values[usable_flags], na.rm = TRUE), NA_real_),
    AverageBlankCorrectedOD = if_else(n_used > 0, mean(corrected_values[usable_flags], na.rm = TRUE), NA_real_),
    ExcludedWellNote = note
  )
}

# Four-parameter logistic model.
# Bottom = lower asymptote
# Top    = upper asymptote
# IC50   = concentration where response is halfway between Top and Bottom
# Hill   = slope/steepness parameter
four_parameter_logistic <- function(x, Bottom, Top, IC50, Hill) {
  Bottom + ((Top - Bottom) / (1 + (x / IC50)^Hill))
}

# Inverse of the 4-parameter logistic model.
# This converts %B/B0 values back into concentrations.
invert_four_parameter_logistic <- function(y, Bottom, Top, IC50, Hill) {
  # Algebraic inverse of:
  # y = Bottom + ((Top - Bottom) / (1 + (x / IC50)^Hill))
  #
  # Important: if y is outside the fitted model asymptotes, there is no exact
  # real-valued solution. In that case this function returns NA rather than
  # a misleading complex/NaN value.
  ratio <- ((Top - Bottom) / (y - Bottom)) - 1
  x <- IC50 * (ratio^(1 / Hill))
  x[!is.finite(x) | !is.numeric(x)] <- NA_real_
  as.numeric(x)
}

# Local tail extrapolation used only when a sample response is outside the fitted
# 4PL asymptotes and therefore cannot be inverted exactly. This extends the
# relationship between %B/B0 and log10(concentration) using the two nearest
# nonzero standards at the relevant end of the curve. Values are flagged clearly
# and should be considered exploratory only.
tail_extrapolate_log_linear <- function(y, slope, intercept) {
  x <- 10^(intercept + slope * y)
  x[!is.finite(x) | x <= 0] <- NA_real_
  as.numeric(x)
}

6 Multi-plate processing logic

For multi-plate analyses, the script uses Plate Number as a fixed metadata column. This means the same well address can safely appear on multiple plates. For example, A3 on Plate 1 and A3 on Plate 2 are treated as different observations.

The plate-specific calculations are intentionally kept separate until each sample has a final concentration:

7 Import and tidy the raw Tecan plate data

The Tecan export contains instrument settings before the actual plate matrix. This code automatically finds the plate data by locating the row where the first cell is <>, which is followed by columns 1-12 and rows A-H.

# Each Excel file is read separately and tagged with its Plate Number before being combined.
# This prevents identical well names on different plates, such as A3 on Plate 1 and A3 on Plate 2,
# from being confused.

if (is.null(names(plate_excel_files)) || any(names(plate_excel_files) == "")) {
  stop("plate_excel_files must be a named vector, for example c('1' = 'Plate1.xlsx', '2' = 'Plate2.xlsx').")
}

plate_file_map <- tibble(
  `Plate Number` = normalize_plate_number(names(plate_excel_files)),
  TecanExcelFile = as.character(plate_excel_files)
)

plate_long <- pmap_dfr(
  plate_file_map,
  function(`Plate Number`, TecanExcelFile) {
    read_tecan_plate(TecanExcelFile, plate_number = `Plate Number`)
  }
)

plate_long
# A tibble: 96 × 6
   `Plate Number` PlateWellID Well  PlateRow PlateCol RawOD
   <chr>          <chr>       <chr> <chr>       <int> <dbl>
 1 1              1:A1        A1    A               1 0.374
 2 1              1:A2        A2    A               2 0.517
 3 1              1:A3        A3    A               3 0.498
 4 1              1:A4        A4    A               4 0.572
 5 1              1:A5        A5    A               5 0.520
 6 1              1:A6        A6    A               6 0.507
 7 1              1:A7        A7    A               7 0.507
 8 1              1:A8        A8    A               8 0.586
 9 1              1:A9        A9    A               9 0.551
10 1              1:A10       A10   A              10 0.548
# ℹ 86 more rows

8 Import metadata

The metadata file must contain at least these columns:

The row with SampleID == "BLK" is used to identify the two blank wells.

metadata <- read_csv(metadata_csv_file, show_col_types = FALSE)

required_columns <- c("SampleID", "Cell Position 1", "Cell Position 2", "Plate Number")
missing_columns <- setdiff(required_columns, colnames(metadata))

if (length(missing_columns) > 0) {
  stop("The metadata file is missing required column(s): ", paste(missing_columns, collapse = ", "))
}

metadata <- metadata %>%
  mutate(
    SampleID = as.character(SampleID),
    `Plate Number` = normalize_plate_number(`Plate Number`),
    `Cell Position 1` = str_to_upper(str_trim(as.character(`Cell Position 1`))),
    `Cell Position 2` = str_to_upper(str_trim(as.character(`Cell Position 2`)))
  )

metadata_plates <- sort(unique(metadata$`Plate Number`))
plate_files_available <- sort(unique(plate_file_map$`Plate Number`))
missing_plate_files <- setdiff(metadata_plates, plate_files_available)

if (length(missing_plate_files) > 0) {
  stop(
    "The metadata contains Plate Number value(s) with no matching Excel file in plate_excel_files: ",
    paste(missing_plate_files, collapse = ", ")
  )
}

metadata
# A tibble: 40 × 11
   SampleID `Plate Number` `Cell Position 1` `Cell Position 2` `Dilution Factor`
   <chr>    <chr>          <chr>             <chr>                         <dbl>
 1 D-19_06… 1              A3                A4                              0  
 2 D-19_06… 1              B3                B4                             35  
 3 D-19_06… 1              C3                C4                             30  
 4 D-48_06… 1              D3                D4                              0  
 5 D-48_06… 1              E3                E4                             26.5
 6 D-48_06… 1              F3                F4                              0  
 7 D-CTL-4… 1              G3                G4                             50  
 8 D-43_06… 1              H3                H4                              0  
 9 D-CTL-1… 1              A5                A6                             50  
10 D-CTL-1… 1              B5                B6                             45.5
# ℹ 30 more rows
# ℹ 6 more variables:
#   SampleTissue_MassorVolume_ForCortisolExtraction_mg_or_ul <dbl>,
#   ExtractionCorrectionFactor <dbl>, Tissue <chr>, Ploidy <chr>,
#   SampleTime <chr>, Treatment <chr>

9 Register manually excluded wells

Any wells listed in excluded_wells are removed from the averaging step. If only one well from a duplicate pair is excluded, the code uses the remaining single well and records a note in the outputs.

excluded_well_table <- make_excluded_well_table(
  excluded_wells = excluded_wells,
  plate_numbers = plate_file_map$`Plate Number`
)

if (nrow(excluded_well_table) == 0) {
  cat("No wells were manually excluded.\n")
} else {
  cat("Manually excluded plate/well combinations:\n")
  print(excluded_well_table)
}
No wells were manually excluded.

10 Identify blanks and calculate average raw OD of blanks

blank_rows <- metadata %>%
  filter(str_to_upper(SampleID) == "BLK")

blank_counts <- blank_rows %>%
  count(`Plate Number`, name = "NumberOfBlankRows")

plates_missing_blank <- setdiff(metadata_plates, blank_counts$`Plate Number`)
plates_with_multiple_blanks <- blank_counts %>%
  filter(NumberOfBlankRows != 1) %>%
  pull(`Plate Number`)

if (length(plates_missing_blank) > 0) {
  stop("Each plate must have exactly one metadata row where SampleID is 'BLK'. Missing for plate(s): ", paste(plates_missing_blank, collapse = ", "))
}

if (length(plates_with_multiple_blanks) > 0) {
  stop("Each plate must have exactly one metadata row where SampleID is 'BLK'. Multiple blank rows found for plate(s): ", paste(plates_with_multiple_blanks, collapse = ", "))
}

blank_summary <- blank_rows %>%
  mutate(
    BlankWells = map2(`Cell Position 1`, `Cell Position 2`, ~ c(.x, .y)),
    BlankRawValues = map2(
      BlankWells,
      `Plate Number`,
      ~ map_dbl(.x, get_od_from_well, plate_number = .y, plate_long = plate_long, value_column = "RawOD")
    ),
    BlankPlateWellIDs = map2(`Plate Number`, BlankWells, ~ paste(.x, .y, sep = ":")),
    BlankExcludedFlags = map(BlankPlateWellIDs, ~ .x %in% excluded_well_table$PlateWellID)
  ) %>%
  rowwise() %>%
  mutate(
    AverageBlankRawOD = {
      flags <- unlist(BlankExcludedFlags)
      vals <- unlist(BlankRawValues)
      if (all(flags)) NA_real_ else mean(vals[!flags], na.rm = TRUE)
    },
    ExcludedBlankWells = paste(unlist(BlankPlateWellIDs)[unlist(BlankExcludedFlags)], collapse = ";"),
    UsedBlankWells = paste(unlist(BlankPlateWellIDs)[!unlist(BlankExcludedFlags)], collapse = ";")
  ) %>%
  ungroup() %>%
  select(`Plate Number`, AverageBlankRawOD, ExcludedBlankWells, UsedBlankWells)

if (any(is.na(blank_summary$AverageBlankRawOD))) {
  bad_plates <- blank_summary %>% filter(is.na(AverageBlankRawOD)) %>% pull(`Plate Number`)
  stop("Both blank wells were manually excluded for plate(s): ", paste(bad_plates, collapse = ", "), ". At least one blank well per plate is required.")
}

cat("Average Raw OD of Blanks by plate:\n")
Average Raw OD of Blanks by plate:
print(blank_summary %>% select(`Plate Number`, AverageBlankRawOD, UsedBlankWells, ExcludedBlankWells))
# A tibble: 1 × 4
  `Plate Number` AverageBlankRawOD UsedBlankWells ExcludedBlankWells
  <chr>                      <dbl> <chr>          <chr>             
1 1                         0.0541 1:H11;1:H12    ""

11 Blank-correct all wells

plate_long <- plate_long %>%
  left_join(blank_summary %>% select(`Plate Number`, AverageBlankRawOD), by = "Plate Number") %>%
  mutate(
    BlankCorrectedOD = RawOD - AverageBlankRawOD,
    ExcludedFromAnalysis = PlateWellID %in% excluded_well_table$PlateWellID
  )

plate_long
# A tibble: 96 × 9
   `Plate Number` PlateWellID Well  PlateRow PlateCol RawOD AverageBlankRawOD
   <chr>          <chr>       <chr> <chr>       <int> <dbl>             <dbl>
 1 1              1:A1        A1    A               1 0.374            0.0541
 2 1              1:A2        A2    A               2 0.517            0.0541
 3 1              1:A3        A3    A               3 0.498            0.0541
 4 1              1:A4        A4    A               4 0.572            0.0541
 5 1              1:A5        A5    A               5 0.520            0.0541
 6 1              1:A6        A6    A               6 0.507            0.0541
 7 1              1:A7        A7    A               7 0.507            0.0541
 8 1              1:A8        A8    A               8 0.586            0.0541
 9 1              1:A9        A9    A               9 0.551            0.0541
10 1              1:A10       A10   A              10 0.548            0.0541
# ℹ 86 more rows
# ℹ 2 more variables: BlankCorrectedOD <dbl>, ExcludedFromAnalysis <lgl>

12 Process standards

The kit standards are assumed to be loaded in duplicate in columns 1 and 2:

Known cortisol concentrations are:

Standard Concentration ng/mL
S0 0.000
S1 0.005
S2 0.020
S3 0.100
S4 0.500
S5 2.000
S6 10.000
S7 50.000 
standard_map <- expand_grid(
  `Plate Number` = plate_file_map$`Plate Number`,
  StandardIndex = 0:7
) %>%
  mutate(
    Standard = paste0("S", StandardIndex),
    BindingName = paste0("B", StandardIndex),
    PlateRow = LETTERS[StandardIndex + 1],
    Well1 = paste0(PlateRow, "1"),
    Well2 = paste0(PlateRow, "2"),
    KnownConcentration_ng_mL = c(0.000, 0.005, 0.020, 0.100, 0.500, 2.000, 10.000, 50.000)[StandardIndex + 1]
  ) %>%
  select(-StandardIndex)

standards_processed <- standard_map %>%
  mutate(
    duplicate_summary = pmap(
      list(Well1, Well2, `Plate Number`),
      summarize_duplicate_wells,
      plate_long = plate_long,
      excluded_well_table = excluded_well_table
    )
  ) %>%
  unnest(duplicate_summary)

# B0 is the maximum binding reference from S0. It is calculated separately for each plate.
B0_by_plate <- standards_processed %>%
  filter(Standard == "S0") %>%
  transmute(`Plate Number`, B0 = AverageBlankCorrectedOD)

if (any(is.na(B0_by_plate$B0) | B0_by_plate$B0 == 0)) {
  bad_plates <- B0_by_plate %>% filter(is.na(B0) | B0 == 0) %>% pull(`Plate Number`)
  stop("B0 could not be calculated or is zero for plate(s): ", paste(bad_plates, collapse = ", "), ". Check standards S0 wells A1 and A2, including any manually excluded wells.")
}

standards_processed <- standards_processed %>%
  left_join(B0_by_plate, by = "Plate Number") %>%
  mutate(Percent_B_B0 = (AverageBlankCorrectedOD / B0) * 100)

standards_processed
# A tibble: 8 × 21
  `Plate Number` Standard BindingName PlateRow Well1 Well2
  <chr>          <chr>    <chr>       <chr>    <chr> <chr>
1 1              S0       B0          A        A1    A2   
2 1              S1       B1          B        B1    B2   
3 1              S2       B2          C        C1    C2   
4 1              S3       B3          D        D1    D2   
5 1              S4       B4          E        E1    E2   
6 1              S5       B5          F        F1    F2   
7 1              S6       B6          G        G1    G2   
8 1              S7       B7          H        H1    H2   
# ℹ 15 more variables: KnownConcentration_ng_mL <dbl>, RawOD_1 <dbl>,
#   RawOD_2 <dbl>, BlankCorrectedOD_1 <dbl>, BlankCorrectedOD_2 <dbl>,
#   Well1_Excluded <lgl>, Well2_Excluded <lgl>, ExcludedWellsForSample <chr>,
#   WellsUsedForAverage <chr>, NumberOfWellsUsedForAverage <int>,
#   AverageRawOD <dbl>, AverageBlankCorrectedOD <dbl>, ExcludedWellNote <chr>,
#   B0 <dbl>, Percent_B_B0 <dbl>

13 Fit the 4-parameter standard curve

S0/B0 is used as the 100% maximum binding reference. The curve itself is fit with S1-S7 because S0 has a known concentration of zero and the model is concentration-based.

standards_for_curve <- standards_processed %>%
  filter(Standard != "S0")

fit_one_plate_curve <- function(plate_standard_data) {
  # Cortisol assays are competitive ELISAs, so the standard curve should be
  # monotonic decreasing: higher %B/B0 = lower cortisol concentration.
  #
  # This re-parameterized 4PL enforces:
  #   Top > Bottom, IC50 > 0, and Hill > 0
  # so the curve cannot accidentally flip direction during fitting.
  response_span_start <- max(plate_standard_data$Percent_B_B0, na.rm = TRUE) -
    min(plate_standard_data$Percent_B_B0, na.rm = TRUE)

  if (!is.finite(response_span_start) || response_span_start <= 0) {
    stop("The standard responses do not span a usable range for this plate. Check standards and blank correction.")
  }

  fit_4pl <- nlsLM(
    Percent_B_B0 ~ Bottom +
      (exp(LogSpan) / (1 + (KnownConcentration_ng_mL / exp(LogIC50))^exp(LogHill))),
    data = plate_standard_data,
    start = list(
      Bottom = min(plate_standard_data$Percent_B_B0, na.rm = TRUE),
      LogSpan = log(response_span_start),
      LogIC50 = log(stats::median(plate_standard_data$KnownConcentration_ng_mL, na.rm = TRUE)),
      LogHill = log(1)
    ),
    control = nls.lm.control(maxiter = 1000)
  )

  curve_parameters <- coef(fit_4pl)

  tibble(
    Bottom = unname(curve_parameters["Bottom"]),
    Top = unname(curve_parameters["Bottom"] + exp(curve_parameters["LogSpan"])),
    IC50 = unname(exp(curve_parameters["LogIC50"])),
    Hill = unname(exp(curve_parameters["LogHill"])),
    CurveDirection = "Competitive/decreasing: higher %B/B0 corresponds to lower cortisol"
  )
}

curve_parameters_by_plate <- standards_for_curve %>%
  group_by(`Plate Number`) %>%
  group_modify(~ fit_one_plate_curve(.x)) %>%
  ungroup()

cat("4-parameter logistic equation used for each plate:\n")
4-parameter logistic equation used for each plate:
cat("%B/B0 = Bottom + ((Top - Bottom) / (1 + (Concentration / IC50)^Hill))\n\n")
%B/B0 = Bottom + ((Top - Bottom) / (1 + (Concentration / IC50)^Hill))
cat("Fitted parameters by plate:\n")
Fitted parameters by plate:
print(curve_parameters_by_plate)
# A tibble: 1 × 6
  `Plate Number` Bottom   Top   IC50  Hill CurveDirection                       
  <chr>           <dbl> <dbl>  <dbl> <dbl> <chr>                                
1 1                4.14  120. 0.0189 0.579 Competitive/decreasing: higher %B/B0…
# Observed standard curve range by plate.
# This is used to flag concentrations that are extrapolated beyond the measured
# standards but still calculable from the fitted 4PL equation.
tail_coefficients_by_plate <- standards_for_curve %>%
  group_by(`Plate Number`) %>%
  group_modify(~ {
    d <- .x %>%
      filter(KnownConcentration_ng_mL > 0, is.finite(Percent_B_B0)) %>%
      arrange(KnownConcentration_ng_mL)

    if (nrow(d) < 2) {
      return(tibble(
        LowResponseTailSlope = NA_real_,
        LowResponseTailIntercept = NA_real_,
        HighResponseTailSlope = NA_real_,
        HighResponseTailIntercept = NA_real_,
        TailExtrapolationDirectionRule = "Insufficient standards for monotonic tail extrapolation"
      ))
    }

    # Competitive cortisol ELISA direction:
    #   higher %B/B0  -> lower cortisol concentration
    #   lower %B/B0   -> higher cortisol concentration
    #
    # To avoid the extrapolated tail flipping direction because of noisy standards,
    # each tail line is anchored at the nearest measured standard boundary and the
    # slope is forced to be negative when enforce_competitive_monotonic_tail_extrapolation = TRUE.
    # This ensures that samples above the high-binding end of the curve get lower,
    # not higher, concentrations as %B/B0 increases.

    high_response_boundary <- d %>% slice_max(order_by = Percent_B_B0, n = 1, with_ties = FALSE)
    high_response_neighbor <- d %>%
      filter(row_number() != match(high_response_boundary$KnownConcentration_ng_mL, KnownConcentration_ng_mL)) %>%
      slice_max(order_by = Percent_B_B0, n = 1, with_ties = FALSE)

    low_response_boundary <- d %>% slice_min(order_by = Percent_B_B0, n = 1, with_ties = FALSE)
    low_response_neighbor <- d %>%
      filter(row_number() != match(low_response_boundary$KnownConcentration_ng_mL, KnownConcentration_ng_mL)) %>%
      slice_min(order_by = Percent_B_B0, n = 1, with_ties = FALSE)

    calculate_tail_line <- function(boundary, neighbor) {
      dy <- neighbor$Percent_B_B0 - boundary$Percent_B_B0
      dx <- log10(neighbor$KnownConcentration_ng_mL) - log10(boundary$KnownConcentration_ng_mL)
      slope <- dx / dy

      if (!is.finite(slope) || dy == 0) {
        slope <- NA_real_
      }

      if (enforce_competitive_monotonic_tail_extrapolation && is.finite(slope)) {
        slope <- -abs(slope)
      }

      intercept <- log10(boundary$KnownConcentration_ng_mL) - slope * boundary$Percent_B_B0
      tibble(Slope = slope, Intercept = intercept)
    }

    high_tail <- calculate_tail_line(high_response_boundary, high_response_neighbor)
    low_tail <- calculate_tail_line(low_response_boundary, low_response_neighbor)

    tibble(
      LowResponseTailSlope = low_tail$Slope,
      LowResponseTailIntercept = low_tail$Intercept,
      HighResponseTailSlope = high_tail$Slope,
      HighResponseTailIntercept = high_tail$Intercept,
      TailExtrapolationDirectionRule = if_else(
        enforce_competitive_monotonic_tail_extrapolation,
        "Strict competitive monotonic tail enforced: higher %B/B0 gives lower concentration",
        "Tail slopes based only on nearest standards; monotonic direction not forced"
      )
    )
  }) %>%
  ungroup()

standard_curve_range_by_plate <- standards_for_curve %>%
  group_by(`Plate Number`) %>%
  summarize(
    MinStandardConcentration_ng_mL = min(KnownConcentration_ng_mL, na.rm = TRUE),
    MaxStandardConcentration_ng_mL = max(KnownConcentration_ng_mL, na.rm = TRUE),
    MinStandardPercent_B_B0 = min(Percent_B_B0, na.rm = TRUE),
    MaxStandardPercent_B_B0 = max(Percent_B_B0, na.rm = TRUE),
    # These two columns make boundary estimates direction-safe.
    # For a competitive ELISA, the lowest %B/B0 usually corresponds to the highest
    # concentration, but this is calculated from the standards rather than assumed.
    BoundaryConcentrationAtMinPercent_ng_mL = KnownConcentration_ng_mL[which.min(Percent_B_B0)],
    BoundaryConcentrationAtMaxPercent_ng_mL = KnownConcentration_ng_mL[which.max(Percent_B_B0)],
    .groups = "drop"
  ) %>%
  left_join(tail_coefficients_by_plate, by = "Plate Number")

cat("Observed standard curve range by plate:\n")
Observed standard curve range by plate:
print(standard_curve_range_by_plate)
# A tibble: 1 × 12
  `Plate Number` MinStandardConcentration_ng_mL MaxStandardConcentration_ng_mL
  <chr>                                   <dbl>                          <dbl>
1 1                                       0.005                             50
# ℹ 9 more variables: MinStandardPercent_B_B0 <dbl>,
#   MaxStandardPercent_B_B0 <dbl>,
#   BoundaryConcentrationAtMinPercent_ng_mL <dbl>,
#   BoundaryConcentrationAtMaxPercent_ng_mL <dbl>, LowResponseTailSlope <dbl>,
#   LowResponseTailIntercept <dbl>, HighResponseTailSlope <dbl>,
#   HighResponseTailIntercept <dbl>, TailExtrapolationDirectionRule <chr>

14 Plate-specific processing summary

This table is a quick check that each plate was processed independently. It prints and saves the blank average, B0, and 4-parameter curve coefficients for every plate.

plate_processing_summary <- blank_summary %>%
  select(`Plate Number`, AverageBlankRawOD, UsedBlankWells, ExcludedBlankWells) %>%
  left_join(B0_by_plate, by = "Plate Number") %>%
  left_join(curve_parameters_by_plate, by = "Plate Number") %>%
  left_join(standard_curve_range_by_plate, by = "Plate Number") %>%
  arrange(`Plate Number`)

cat("Plate-specific processing summary. Each row represents an independently processed plate:\n")
Plate-specific processing summary. Each row represents an independently processed plate:
print(plate_processing_summary)
# A tibble: 1 × 21
  `Plate Number` AverageBlankRawOD UsedBlankWells ExcludedBlankWells    B0
  <chr>                      <dbl> <chr>          <chr>              <dbl>
1 1                         0.0541 1:H11;1:H12    ""                 0.392
# ℹ 16 more variables: Bottom <dbl>, Top <dbl>, IC50 <dbl>, Hill <dbl>,
#   CurveDirection <chr>, MinStandardConcentration_ng_mL <dbl>,
#   MaxStandardConcentration_ng_mL <dbl>, MinStandardPercent_B_B0 <dbl>,
#   MaxStandardPercent_B_B0 <dbl>,
#   BoundaryConcentrationAtMinPercent_ng_mL <dbl>,
#   BoundaryConcentrationAtMaxPercent_ng_mL <dbl>, LowResponseTailSlope <dbl>,
#   LowResponseTailIntercept <dbl>, HighResponseTailSlope <dbl>, …
write_csv(plate_processing_summary, plate_processing_summary_output)
cat("Plate processing summary written to: ", plate_processing_summary_output, "\n", sep = "")
Plate processing summary written to: cortisol_plate_processing_summary.csv

15 Plot the standard curve

make_standard_curve_plot <- function(plate_number) {
  this_standards_all <- standards_processed %>%
    filter(`Plate Number` == plate_number)
  
  this_standards_curve <- standards_for_curve %>%
    filter(`Plate Number` == plate_number)
  
  this_params <- curve_parameters_by_plate %>%
    filter(`Plate Number` == plate_number)
  
  curve_data <- tibble(
    KnownConcentration_ng_mL = exp(seq(
      log(min(this_standards_curve$KnownConcentration_ng_mL)),
      log(max(this_standards_curve$KnownConcentration_ng_mL)),
      length.out = 300
    ))
  ) %>%
    mutate(
      Predicted_Percent_B_B0 = four_parameter_logistic(
        KnownConcentration_ng_mL,
        Bottom = this_params$Bottom,
        Top = this_params$Top,
        IC50 = this_params$IC50,
        Hill = this_params$Hill
      )
    )
  
  ggplot() +
    geom_point(
      data = this_standards_all,
      aes(x = KnownConcentration_ng_mL, y = Percent_B_B0),
      size = 3
    ) +
    geom_line(
      data = curve_data,
      aes(x = KnownConcentration_ng_mL, y = Predicted_Percent_B_B0),
      linewidth = 1
    ) +
    scale_x_continuous(
      trans = scales::pseudo_log_trans(sigma = 0.001),
      breaks = c(0, 0.005, 0.02, 0.1, 0.5, 2, 10, 50)
    ) +
    labs(
      title = paste("Cortisol ELISA 4-Parameter Standard Curve - Plate", plate_number),
      x = "Known cortisol concentration (ng/mL)",
      y = "%B/B0"
    ) +
    theme_bw()
}

standard_curve_plots <- map(plate_file_map$`Plate Number`, make_standard_curve_plot)

# Print plots in the knitted report.
walk(standard_curve_plots, print)

# Save all standard curves into one multi-page PDF.
pdf(standard_curve_plot_output, width = 7, height = 5)
walk(standard_curve_plots, print)
dev.off()
quartz_off_screen 
                2 
cat("Standard curve plot(s) written to: ", standard_curve_plot_output, "\n", sep = "")
Standard curve plot(s) written to: cortisol_standard_curve.pdf

16 Process samples

For each sample listed in the metadata file, the two duplicate wells are extracted from the correct plate, blank-corrected with that plate’s blank average, averaged, converted to %B/B0 using that plate’s B0, and back-calculated to a cortisol concentration from that plate’s own 4-parameter standard curve.

The blank row is retained in the full output, but concentration calculations are only meaningful for actual samples.

processed_metadata <- metadata %>%
  mutate(
    duplicate_summary = pmap(
      list(`Cell Position 1`, `Cell Position 2`, `Plate Number`),
      summarize_duplicate_wells,
      plate_long = plate_long,
      excluded_well_table = excluded_well_table
    )
  ) %>%
  unnest(duplicate_summary) %>%
  left_join(B0_by_plate, by = "Plate Number") %>%
  left_join(curve_parameters_by_plate, by = "Plate Number") %>%
  left_join(standard_curve_range_by_plate, by = "Plate Number") %>%
  mutate(
    Percent_B_B0 = (AverageBlankCorrectedOD / B0) * 100
  )

# Add sample-specific or default dilution factors.
# If these columns already exist in the metadata file, keep and use them.
# If not, create them using the default values specified above.
if (!("DilutionFactor" %in% colnames(processed_metadata))) {
  processed_metadata <- processed_metadata %>%
    mutate(DilutionFactor = default_sample_dilution_factor)
}

if (!("ExtractionCorrectionFactor" %in% colnames(processed_metadata))) {
  processed_metadata <- processed_metadata %>%
    mutate(ExtractionCorrectionFactor = default_extraction_correction_factor)
}

processed_metadata <- processed_metadata %>%
  mutate(
    IsBlank = str_to_upper(SampleID) == "BLK",
    IsStandard = str_detect(str_to_upper(SampleID), "^S[0-7]$"),
    Percent_B_B0_Warning = case_when(
      IsBlank ~ NA_character_,
      Percent_B_B0 < 20 ~ "Below recommended 20% B/B0 range; dilute/rerun or interpret cautiously.",
      Percent_B_B0 > 80 ~ "Above recommended 80% B/B0 range; dilute/rerun or interpret cautiously.",
      TRUE ~ NA_character_
    ),
    # Direction-independent check for whether a sample falls within the observed
    # standard %B/B0 range on its own plate. Competitive ELISAs often decrease
    # with concentration, so use min/max rather than assuming a direction.
    WithinObservedStandardPercentRange = !is.na(Percent_B_B0) &
      Percent_B_B0 >= MinStandardPercent_B_B0 &
      Percent_B_B0 <= MaxStandardPercent_B_B0,
    RawCurveConcentration_ng_mL = if_else(
      IsBlank,
      NA_real_,
      invert_four_parameter_logistic(
        Percent_B_B0,
        Bottom = Bottom,
        Top = Top,
        IC50 = IC50,
        Hill = Hill
      )
    ),
    # The fitted 4PL equation has upper and lower asymptotes. If a sample %B/B0
    # is beyond those fitted asymptotes, the exact 4PL inverse cannot return a
    # real concentration even if the sample is only slightly above the preferred
    # 20-80% range. In that case, optional log-linear tail extrapolation is used.
    FittedMinAsymptotePercent_B_B0 = pmin(Bottom, Top),
    FittedMaxAsymptotePercent_B_B0 = pmax(Bottom, Top),
    NeedsHighResponseTailExtrapolation = is.na(RawCurveConcentration_ng_mL) &
      !is.na(Percent_B_B0) &
      Percent_B_B0 > FittedMaxAsymptotePercent_B_B0,
    NeedsLowResponseTailExtrapolation = is.na(RawCurveConcentration_ng_mL) &
      !is.na(Percent_B_B0) &
      Percent_B_B0 < FittedMinAsymptotePercent_B_B0,
    # Fallback: if the 4PL inverse returns NA for another numerical reason, use
    # the nearest observed standard-response tail. This prevents samples such as
    # high-%B/B0 values above the preferred 80% range from being left as NA when
    # exploratory extrapolation is enabled.
    NearestTailForUnsolvedSample = case_when(
      is.na(RawCurveConcentration_ng_mL) & !is.na(Percent_B_B0) &
        abs(Percent_B_B0 - MaxStandardPercent_B_B0) <= abs(Percent_B_B0 - MinStandardPercent_B_B0) ~ "high_response",
      is.na(RawCurveConcentration_ng_mL) & !is.na(Percent_B_B0) ~ "low_response",
      TRUE ~ NA_character_
    ),
    CurveBackCalculationType = case_when(
      IsBlank ~ "Blank - not calculated",
      is.na(Percent_B_B0) ~ "No %B/B0 available",
      !is.na(RawCurveConcentration_ng_mL) & WithinObservedStandardPercentRange ~ "Interpolated within observed standard range",
      !is.na(RawCurveConcentration_ng_mL) & !WithinObservedStandardPercentRange ~ "4PL equation extrapolated outside observed standard range - interpret cautiously",
      use_tail_extrapolation_outside_model_asymptotes & NeedsHighResponseTailExtrapolation ~ "Outside fitted 4PL upper asymptote; log-linear low-concentration/high-binding tail extrapolation - exploratory only",
      use_tail_extrapolation_outside_model_asymptotes & NeedsLowResponseTailExtrapolation ~ "Outside fitted 4PL lower asymptote; log-linear high-concentration/low-binding tail extrapolation - exploratory only",
      use_tail_extrapolation_outside_model_asymptotes & is.na(RawCurveConcentration_ng_mL) & NearestTailForUnsolvedSample == "high_response" ~ "4PL inverse failed numerically; log-linear low-concentration/high-binding tail extrapolation - exploratory only",
      use_tail_extrapolation_outside_model_asymptotes & is.na(RawCurveConcentration_ng_mL) & NearestTailForUnsolvedSample == "low_response" ~ "4PL inverse failed numerically; log-linear high-concentration/low-binding tail extrapolation - exploratory only",
      is.na(RawCurveConcentration_ng_mL) & use_boundary_estimates_outside_model_asymptotes & Percent_B_B0 > MaxStandardPercent_B_B0 ~ "Outside fitted curve/asymptote; boundary estimate at low concentration limit",
      is.na(RawCurveConcentration_ng_mL) & use_boundary_estimates_outside_model_asymptotes & Percent_B_B0 < MinStandardPercent_B_B0 ~ "Outside fitted curve/asymptote; boundary estimate at high concentration limit",
      is.na(RawCurveConcentration_ng_mL) ~ "Outside fitted curve/asymptote; no exact concentration solution",
      TRUE ~ "Unclassified curve result"
    ),
    TailExtrapolatedConcentration_ng_mL = case_when(
      IsBlank ~ NA_real_,
      !use_tail_extrapolation_outside_model_asymptotes ~ NA_real_,
      is.na(Percent_B_B0) ~ NA_real_,
      is.finite(RawCurveConcentration_ng_mL) ~ NA_real_,
      NeedsHighResponseTailExtrapolation | NearestTailForUnsolvedSample == "high_response" ~ tail_extrapolate_log_linear(
        Percent_B_B0,
        slope = HighResponseTailSlope,
        intercept = HighResponseTailIntercept
      ),
      NeedsLowResponseTailExtrapolation | NearestTailForUnsolvedSample == "low_response" ~ tail_extrapolate_log_linear(
        Percent_B_B0,
        slope = LowResponseTailSlope,
        intercept = LowResponseTailIntercept
      ),
      TRUE ~ NA_real_
    ),
    BoundaryEstimateConcentration_ng_mL = case_when(
      IsBlank ~ NA_real_,
      !use_boundary_estimates_outside_model_asymptotes ~ NA_real_,
      is.na(Percent_B_B0) ~ NA_real_,
      # If the sample response is above the observed standard-response range,
      # assign the concentration of the standard at the high-response boundary.
      Percent_B_B0 > MaxStandardPercent_B_B0 ~ BoundaryConcentrationAtMaxPercent_ng_mL,
      # If the sample response is below the observed standard-response range,
      # assign the concentration of the standard at the low-response boundary.
      # In a typical competitive cortisol ELISA this is the highest standard, S7.
      Percent_B_B0 < MinStandardPercent_B_B0 ~ BoundaryConcentrationAtMinPercent_ng_mL,
      TRUE ~ NA_real_
    ),
    Concentration_ng_mL = case_when(
      IsBlank ~ NA_real_,
      !allow_curve_extrapolation & !WithinObservedStandardPercentRange ~ NA_real_,
      is.finite(RawCurveConcentration_ng_mL) ~ RawCurveConcentration_ng_mL,
      use_tail_extrapolation_outside_model_asymptotes & is.finite(TailExtrapolatedConcentration_ng_mL) ~ TailExtrapolatedConcentration_ng_mL,
      use_boundary_estimates_outside_model_asymptotes & is.finite(BoundaryEstimateConcentration_ng_mL) ~ BoundaryEstimateConcentration_ng_mL,
      TRUE ~ NA_real_
    ),
    UsesExtrapolatedConcentration = str_detect(CurveBackCalculationType, "extrapolat"),
    UsesTailExtrapolation = str_detect(CurveBackCalculationType, "tail extrapolation"),
    UsesBoundaryEstimate = str_detect(CurveBackCalculationType, "boundary estimate"),
    QuantificationConfidence = case_when(
      IsBlank ~ "Blank/control - not quantified",
      IsStandard ~ "Standard/control - not a research sample",
      is.na(Concentration_ng_mL) ~ "No reportable concentration under current settings",
      UsesBoundaryEstimate ~ "Censored boundary estimate - interpret with extreme caution",
      UsesTailExtrapolation ~ "Very low confidence - log-linear tail extrapolation outside fitted 4PL asymptote",
      UsesExtrapolatedConcentration ~ "Low confidence - 4PL extrapolation outside measured standard range",
      WithinObservedStandardPercentRange &
        Percent_B_B0 >= preferred_percent_B_B0_lower &
        Percent_B_B0 <= preferred_percent_B_B0_upper ~ "High confidence - within measured standard range and preferred 20-80% B/B0 range",
      WithinObservedStandardPercentRange ~ "Moderate confidence - within measured standard range but outside preferred 20-80% B/B0 range",
      TRUE ~ "Unclassified quantification confidence"
    ),
    ConcentrationInterpretationNote = case_when(
      IsBlank ~ NA_character_,
      CurveBackCalculationType == "Interpolated within observed standard range" ~ "Best-supported value: sample %B/B0 is within the measured standard curve range.",
      UsesTailExtrapolation ~ "Exploratory-only value: sample is outside the fitted 4PL asymptote, so the exact 4PL equation cannot be inverted. This value uses a local log-linear tail extrapolation from the nearest two standards; rerun at a better dilution for final interpretation.",
      UsesExtrapolatedConcentration ~ "Cautious/preliminary value: sample is outside the measured standard range but still has a finite 4PL back-calculated concentration.",
      UsesBoundaryEstimate ~ "Censored/boundary estimate: sample response is outside the fitted curve/asymptote, so this is not an exact concentration. Treat as approximate and rerun at a better dilution.",
      is.na(Concentration_ng_mL) ~ "No concentration reported because the value cannot be reliably solved from the fitted curve under current settings.",
      TRUE ~ NA_character_
    ),
    AssayRerunRecommendation = case_when(
      IsBlank ~ "Blank/control - no rerun recommendation",
      IsStandard ~ "Standard/control - no rerun recommendation",
      is.na(Percent_B_B0) ~ "Cannot evaluate - missing %B/B0",
      Percent_B_B0 > preferred_percent_B_B0_upper ~ "Signal too high / cortisol likely low: rerun less diluted or concentrate sample if accurate quantification is needed",
      Percent_B_B0 < preferred_percent_B_B0_lower ~ "Signal too low / cortisol likely high: dilute sample more and rerun if accurate quantification is needed",
      TRUE ~ "Within preferred 20-80% B/B0 working range"
    ),
    CompetitiveELISAInterpretation = case_when(
      IsBlank ~ "Blank/control",
      IsStandard ~ "Standard/control",
      is.na(Percent_B_B0) ~ "Missing %B/B0",
      Percent_B_B0 > preferred_percent_B_B0_upper ~ "High OD/%B/B0 = low cortisol; near top plateau",
      Percent_B_B0 < preferred_percent_B_B0_lower ~ "Low OD/%B/B0 = high cortisol; near bottom plateau",
      TRUE ~ "Preferred working range"
    ),
    FinalConcentration_ng_mL = Concentration_ng_mL * DilutionFactor * ExtractionCorrectionFactor,
    CompetitiveELISADirectionCheck = case_when(
      IsBlank ~ "Blank - not evaluated",
      IsStandard ~ "Standard/control - not evaluated as sample",
      is.na(Percent_B_B0) | is.na(Concentration_ng_mL) ~ "Cannot evaluate direction - missing value",
      Percent_B_B0 > preferred_percent_B_B0_upper & Concentration_ng_mL > IC50 ~
        "CHECK: high %B/B0 should indicate low cortisol; concentration is above plate IC50",
      Percent_B_B0 < preferred_percent_B_B0_lower & Concentration_ng_mL < IC50 ~
        "CHECK: low %B/B0 should indicate high cortisol; concentration is below plate IC50",
      TRUE ~ "Direction consistent with competitive ELISA"
    )
  ) %>%
  mutate(
    ProcessingStatus = case_when(
      IsBlank ~ "Blank control; excluded from concentration statistics.",
      IsStandard ~ "Standard/control well; used for curve fitting but excluded from sample statistics.",
      NumberOfWellsUsedForAverage == 0 ~ "No usable duplicate wells after exclusions or missing OD values.",
      is.na(AverageBlankCorrectedOD) ~ "Missing average blank-corrected OD.",
      is.na(B0) ~ "Missing B0 for this plate; check standards on this plate.",
      is.na(Percent_B_B0) ~ "Could not calculate %B/B0.",
      is.na(Concentration_ng_mL) ~ "Could not back-calculate concentration from the plate-specific 4PL curve under current settings.",
      !is.finite(Concentration_ng_mL) ~ "Back-calculated concentration was not finite.",
      is.na(DilutionFactor) | is.na(ExtractionCorrectionFactor) ~ "Missing dilution or extraction correction factor.",
      is.na(FinalConcentration_ng_mL) ~ "Final concentration is missing after corrections.",
      UsesBoundaryEstimate ~ "Usable only as a cautious boundary estimate; rerun at a better dilution for final interpretation.",
      UsesExtrapolatedConcentration ~ "Usable as cautious extrapolated concentration; rerun if this will be used for final interpretation.",
      TRUE ~ "Usable sample concentration."
    ),
    IncludedInPlotsAndStatistics = !IsBlank &
      (!exclude_standards_from_plots_and_statistics | !IsStandard) &
      !(SampleID %in% exclude_sample_ids_from_plots_and_statistics) &
      !is.na(FinalConcentration_ng_mL) &
      (include_extrapolated_concentrations_in_statistics | (!UsesExtrapolatedConcentration & !UsesBoundaryEstimate))
  )

processed_metadata
# A tibble: 40 × 67
   SampleID `Plate Number` `Cell Position 1` `Cell Position 2` `Dilution Factor`
   <chr>    <chr>          <chr>             <chr>                         <dbl>
 1 D-19_06… 1              A3                A4                              0  
 2 D-19_06… 1              B3                B4                             35  
 3 D-19_06… 1              C3                C4                             30  
 4 D-48_06… 1              D3                D4                              0  
 5 D-48_06… 1              E3                E4                             26.5
 6 D-48_06… 1              F3                F4                              0  
 7 D-CTL-4… 1              G3                G4                             50  
 8 D-43_06… 1              H3                H4                              0  
 9 D-CTL-1… 1              A5                A6                             50  
10 D-CTL-1… 1              B5                B6                             45.5
# ℹ 30 more rows
# ℹ 62 more variables:
#   SampleTissue_MassorVolume_ForCortisolExtraction_mg_or_ul <dbl>,
#   ExtractionCorrectionFactor <dbl>, Tissue <chr>, Ploidy <chr>,
#   SampleTime <chr>, Treatment <chr>, RawOD_1 <dbl>, RawOD_2 <dbl>,
#   BlankCorrectedOD_1 <dbl>, BlankCorrectedOD_2 <dbl>, Well1_Excluded <lgl>,
#   Well2_Excluded <lgl>, ExcludedWellsForSample <chr>, …

16.1 Interpretation of extrapolated or boundary-estimated concentrations

The column CurveBackCalculationType separates concentration values into practical calculation categories:

  • Interpolated within observed standard range: strongest values; the sample %B/B0 falls within the measured standard curve range for that plate.
  • 4PL equation extrapolated outside observed standard range: calculated from the fitted 4PL equation, but beyond the measured standards; useful for preliminary screening, but interpret cautiously and rerun at a better dilution for final work.
  • Outside fitted 4PL asymptote; log-linear tail extrapolation: the exact 4PL inverse has no real-valued solution for that %B/B0, so the script uses the nearest two nonzero standards at that end of the curve to preserve variation for exploratory visualization. This is not a true 4PL concentration and should be rerun at a better dilution for final interpretation.
  • Outside fitted curve/asymptote; boundary estimate: optional fallback only. The script reports the nearest standard-boundary estimate if use_boundary_estimates_outside_model_asymptotes <- TRUE. These are censored/approximate values, not exact concentrations, and are off by default because they can collapse many samples to identical values.

The column QuantificationConfidence gives a quick user-facing interpretation of those categories:

  • High confidence: within the measured standard range and the preferred 20-80% B/B0 range.
  • Moderate confidence: within the measured standard range but outside the preferred 20-80% B/B0 range.
  • Low confidence: 4PL extrapolation outside the measured standard range.
  • Very low confidence: log-linear tail extrapolation outside the fitted 4PL asymptote.
  • Censored boundary estimate: nearest standard-limit estimate only, used only if boundary estimates are enabled.
curve_backcalculation_summary <- processed_metadata %>%
  filter(!IsBlank, !IsStandard) %>%
  count(`Plate Number`, CurveBackCalculationType, name = "n_samples")

curve_backcalculation_summary
# A tibble: 3 × 3
  `Plate Number` CurveBackCalculationType                              n_samples
  <chr>          <chr>                                                     <int>
1 1              4PL equation extrapolated outside observed standard …        11
2 1              Interpolated within observed standard range                   2
3 1              Outside fitted 4PL upper asymptote; log-linear low-c…        26

17 Visualize samples relative to the standard curve

This plot overlays research samples onto the plate-specific standard curve. The shaded region marks the measured standard concentration range. Samples outside this shaded region were not interpolated between standards and should be interpreted cautiously.

make_standard_curve_with_samples_plot <- function(plate_number) {
  this_standards_all <- standards_processed %>%
    filter(`Plate Number` == plate_number)

  this_standards_curve <- standards_for_curve %>%
    filter(`Plate Number` == plate_number)

  this_samples <- processed_metadata %>%
    filter(
      `Plate Number` == plate_number,
      !IsBlank,
      !IsStandard,
      !is.na(Percent_B_B0),
      !is.na(Concentration_ng_mL),
      is.finite(Concentration_ng_mL),
      Concentration_ng_mL > 0
    )

  this_params <- curve_parameters_by_plate %>%
    filter(`Plate Number` == plate_number)

  if (nrow(this_params) == 0 || nrow(this_standards_curve) == 0) {
    return(
      ggplot() +
        annotate("text", x = 1, y = 1, label = paste("No standard curve available for Plate", plate_number)) +
        theme_void()
    )
  }

  min_standard_conc <- min(this_standards_curve$KnownConcentration_ng_mL, na.rm = TRUE)
  max_standard_conc <- max(this_standards_curve$KnownConcentration_ng_mL, na.rm = TRUE)

  # Extend the displayed x-axis to include samples when extrapolated values exist.
  min_plot_conc <- min(c(min_standard_conc, this_samples$Concentration_ng_mL), na.rm = TRUE)
  max_plot_conc <- max(c(max_standard_conc, this_samples$Concentration_ng_mL), na.rm = TRUE)

  # Keep plotting limits finite and positive.
  if (!is.finite(min_plot_conc) || min_plot_conc <= 0) min_plot_conc <- min_standard_conc
  if (!is.finite(max_plot_conc) || max_plot_conc <= 0) max_plot_conc <- max_standard_conc

  curve_data <- tibble(
    KnownConcentration_ng_mL = exp(seq(
      log(min_plot_conc),
      log(max_plot_conc),
      length.out = 500
    ))
  ) %>%
    mutate(
      Predicted_Percent_B_B0 = four_parameter_logistic(
        KnownConcentration_ng_mL,
        Bottom = this_params$Bottom,
        Top = this_params$Top,
        IC50 = this_params$IC50,
        Hill = this_params$Hill
      )
    )

  p <- ggplot()

  if (isTRUE(show_standard_range_shading)) {
    p <- p +
      annotate(
        "rect",
        xmin = min_standard_conc, xmax = max_standard_conc,
        ymin = -Inf, ymax = Inf,
        alpha = 0.10
      )
  }

  p +
    geom_line(
      data = curve_data,
      aes(x = KnownConcentration_ng_mL, y = Predicted_Percent_B_B0),
      linewidth = 1
    ) +
    geom_point(
      data = this_standards_all,
      aes(x = KnownConcentration_ng_mL, y = Percent_B_B0),
      size = 3,
      shape = 17
    ) +
    geom_point(
      data = this_samples,
      aes(x = Concentration_ng_mL, y = Percent_B_B0, color = QuantificationConfidence),
      size = 2.3,
      alpha = 0.85
    ) +
    scale_x_continuous(
      trans = scales::pseudo_log_trans(sigma = 0.001),
      breaks = sort(unique(c(0.005, 0.02, 0.1, 0.5, 2, 10, 50, min_standard_conc, max_standard_conc)))
    ) +
    labs(
      title = paste("Cortisol ELISA Standard Curve with Samples - Plate", plate_number),
      subtitle = "Shaded region = measured standard concentration range; sample points outside this region are extrapolated/cautious.",
      x = "Cortisol concentration from curve (ng/mL before dilution/extraction correction)",
      y = "%B/B0",
      color = "Quantification confidence"
    ) +
    theme_bw() +
    guides(
    color = guide_legend(
      nrow = 4,
      byrow = TRUE,
      title.position = "top"
    )
  ) +
  theme(
    legend.position = "bottom",
    legend.text = element_text(size = 8),
    legend.title = element_text(size = 9)
  )
}

standard_curve_with_samples_plots <- map(plate_file_map$`Plate Number`, make_standard_curve_with_samples_plot)

walk(standard_curve_with_samples_plots, print)

pdf(standard_curve_with_samples_plot_output, width = 8, height = 5)
walk(standard_curve_with_samples_plots, print)
dev.off()
quartz_off_screen 
                2 
cat("Standard curve with sample overlay plot(s) written to: ", standard_curve_with_samples_plot_output, "\n", sep = "")
Standard curve with sample overlay plot(s) written to: cortisol_standard_curve_with_samples.pdf

18 Competitive ELISA interpretation plot with actual samples

This plot is a teaching/QC aid for competitive ELISA interpretation. Cortisol is measured by a competitive assay, so the signal is inverse to analyte concentration: higher OD or %B/B0 means lower cortisol, while lower OD or %B/B0 means higher cortisol.

The plot overlays actual samples onto a conceptual inverse curve and labels the assay action that would usually improve quantification:

make_competitive_interpretation_plot_one_plate <- function(plate_number) {
  this_samples <- processed_metadata %>%
    filter(
      `Plate Number` == plate_number,
      !IsBlank,
      !IsStandard,
      !is.na(Percent_B_B0),
      !is.na(Concentration_ng_mL),
      is.finite(Concentration_ng_mL),
      Concentration_ng_mL > 0
    ) %>%
    mutate(
      AssayActionShort = case_when(
        Percent_B_B0 > preferred_percent_B_B0_upper ~ "Less diluted / concentrate",
        Percent_B_B0 < preferred_percent_B_B0_lower ~ "Dilute more",
        TRUE ~ "Preferred range"
      )
    )

  if (nrow(this_samples) == 0) {
    return(
      ggplot() +
        annotate("text", x = 1, y = 1, label = paste("No quantified samples available for Plate", plate_number)) +
        theme_void()
    )
  }

  min_sample_conc <- min(this_samples$Concentration_ng_mL, na.rm = TRUE)
  max_sample_conc <- max(this_samples$Concentration_ng_mL, na.rm = TRUE)
  min_x <- min(c(0.001, min_sample_conc), na.rm = TRUE)
  max_x <- max(c(100, max_sample_conc), na.rm = TRUE)

  if (!is.finite(min_x) || min_x <= 0) min_x <- 0.001
  if (!is.finite(max_x) || max_x <= min_x) max_x <- 100

  conceptual_curve <- tibble(
    Concentration_ng_mL = exp(seq(log(min_x), log(max_x), length.out = 500))
  ) %>%
    mutate(
      Percent_B_B0 = 100 / (1 + (Concentration_ng_mL / 1)^1.5)
    )

  label_data <- this_samples %>%
    mutate(
      DistanceFromPreferredRange = case_when(
        Percent_B_B0 > preferred_percent_B_B0_upper ~ Percent_B_B0 - preferred_percent_B_B0_upper,
        Percent_B_B0 < preferred_percent_B_B0_lower ~ preferred_percent_B_B0_lower - Percent_B_B0,
        TRUE ~ 0
      )
    ) %>%
    arrange(desc(DistanceFromPreferredRange)) %>%
    slice_head(n = number_of_extreme_samples_to_label) %>%
    filter(label_extreme_samples_on_interpretation_plot, DistanceFromPreferredRange > 0)

  ggplot(conceptual_curve, aes(x = Concentration_ng_mL, y = Percent_B_B0)) +
    annotate("rect", xmin = -Inf, xmax = Inf, ymin = preferred_percent_B_B0_upper, ymax = 100,
             alpha = 0.08) +
    annotate("rect", xmin = -Inf, xmax = Inf, ymin = preferred_percent_B_B0_lower, ymax = preferred_percent_B_B0_upper,
             alpha = 0.08) +
    annotate("rect", xmin = -Inf, xmax = Inf, ymin = 0, ymax = preferred_percent_B_B0_lower,
             alpha = 0.08) +
    geom_hline(yintercept = c(preferred_percent_B_B0_lower, preferred_percent_B_B0_upper), linetype = "dashed") +
    geom_line(linewidth = 1) +
    geom_point(
      data = this_samples,
      aes(x = Concentration_ng_mL, y = Percent_B_B0, color = AssayActionShort),
      inherit.aes = FALSE,
      size = 2.4,
      alpha = 0.85
    ) +
    geom_text(
      data = label_data,
      aes(x = Concentration_ng_mL, y = Percent_B_B0, label = SampleID),
      inherit.aes = FALSE,
      size = 2.8,
      vjust = -0.7,
      check_overlap = TRUE
    ) +
    annotate("text", x = min_x, y = 92, hjust = 0,
             label = "High OD / high %B/B0\nLow cortisol\nLess diluted or concentrate") +
    annotate("text", x = sqrt(min_x * max_x), y = 50, hjust = 0.5,
             label = "Preferred 20-80% B/B0\nBest quantification zone") +
    annotate("text", x = max_x, y = 8, hjust = 1,
             label = "Low OD / low %B/B0\nHigh cortisol\nDilute more") +
    scale_x_continuous(trans = scales::pseudo_log_trans(sigma = 0.001)) +
    labs(
      title = paste("Competitive ELISA Interpretation with Actual Samples - Plate", plate_number),
      subtitle = "Cortisol concentration is inversely related to OD/%B/B0; labels show rerun direction for out-of-range samples.",
      x = "Back-calculated cortisol concentration (ng/mL before dilution/extraction correction)",
      y = "%B/B0 (signal)",
      color = "Assay action"
    ) +
    theme_bw() +
    theme(legend.position = "bottom")
}

if (isTRUE(make_competitive_elisa_interpretation_plot)) {
  competitive_interpretation_plots <- map(plate_file_map$`Plate Number`, make_competitive_interpretation_plot_one_plate)

  walk(competitive_interpretation_plots, print)

  pdf(competitive_elisa_interpretation_plot_output, width = 8, height = 5)
  walk(competitive_interpretation_plots, print)
  dev.off()

  cat("Competitive ELISA interpretation plot(s) written to: ", competitive_elisa_interpretation_plot_output, "\n", sep = "")
}

Competitive ELISA interpretation plot(s) written to: cortisol_competitive_elisa_interpretation_with_samples.pdf

20 Export updated metadata CSV and excluded-well notes

The processed metadata output includes ExcludedWellNote, WellsUsedForAverage, and NumberOfWellsUsedForAverage. These columns document cases where a duplicate average was calculated from only one usable well.

write_csv(processed_metadata, processed_metadata_output)
cat("Processed metadata written to: ", processed_metadata_output, "\n", sep = "")
Processed metadata written to: processed_cortisol_metadata.csv
excluded_well_notes <- processed_metadata %>%
  filter(!is.na(ExcludedWellNote)) %>%
  select(
    `Plate Number`,
    SampleID,
    `Cell Position 1`,
    `Cell Position 2`,
    ExcludedWellsForSample,
    WellsUsedForAverage,
    NumberOfWellsUsedForAverage,
    ExcludedWellNote
  )

if (nrow(excluded_well_notes) == 0) {
  cat("No sample duplicate pairs were affected by manually excluded wells.\n")
} else {
  write_csv(excluded_well_notes, excluded_well_notes_output)
  cat("Excluded-well notes written to: ", excluded_well_notes_output, "\n", sep = "")
  print(excluded_well_notes)
}
No sample duplicate pairs were affected by manually excluded wells.

21 Make jittered boxplots by metadata grouping variables

This section uses the metadata CSV as the guide for group names and comparisons. Only these metadata columns are assumed to be fixed: SampleID, Cell Position 1, Cell Position 2, and Plate Number. Any additional metadata columns are treated as possible grouping variables, such as Treatment, SampleDay, Timepoint, Site, Tank, Exposure, etc.

You can also define optional combined grouping variables in the user-input section, such as TreatmentDay = Treatment + SampleDay. If those columns are not present in a future metadata file, that combined grouping variable is simply skipped.

Important plotting details:

plot_data <- processed_metadata %>%
  # Do not include blanks in any downstream plots or grouping statistics.
  filter(!IsBlank) %>%
  # Standards are curve-calibration wells, not research samples.
  filter(!exclude_standards_from_plots_and_statistics | !IsStandard) %>%
  # Optional manual sample-level exclusions for plotting/statistics only.
  filter(!(SampleID %in% exclude_sample_ids_from_plots_and_statistics)) %>%
  # Statistical tests require a calculated concentration.
  filter(!is.na(FinalConcentration_ng_mL))

# Identify the core metadata columns that are always present and should not be treated
# as experimental grouping variables.
required_metadata_columns <- c("SampleID", "Cell Position 1", "Cell Position 2", "Plate Number")

# Technical columns are generated by this workflow and should not be used as experimental groups.
technical_columns <- c(
  required_metadata_columns,
  "RawOD_1", "RawOD_2", "AverageRawOD",
  "BlankCorrectedOD_1", "BlankCorrectedOD_2", "AverageBlankCorrectedOD",
  "Well1_Excluded", "Well2_Excluded", "ExcludedWellsForSample", "WellsUsedForAverage",
  "NumberOfWellsUsedForAverage", "ExcludedWellNote",
  "Percent_B_B0", "Percent_B_B0_Warning",
  "DilutionFactor", "ExtractionCorrectionFactor",
  "Concentration_ng_mL", "FinalConcentration_ng_mL", "IsBlank", "IsStandard"
)

# Candidate grouping variables come from any metadata columns beyond the required first three columns.
# This keeps the workflow general: Treatment/SampleDay can be used when present, but future metadata
# files can instead use columns such as Site, Timepoint, Sex, Tank, Exposure, Batch, etc.
metadata_grouping_columns <- setdiff(colnames(metadata), required_metadata_columns)
metadata_grouping_columns <- setdiff(metadata_grouping_columns, c("DilutionFactor", "ExtractionCorrectionFactor"))

# If the researcher requested specific grouping variables, use only those columns if they exist.
if (!is.null(grouping_variables_to_analyze)) {
  missing_requested_groups <- setdiff(grouping_variables_to_analyze, colnames(plot_data))
  if (length(missing_requested_groups) > 0) {
    warning(
      "These requested grouping variable(s) were not found in the metadata and will be skipped: ",
      paste(missing_requested_groups, collapse = ", ")
    )
  }
  metadata_grouping_columns <- intersect(grouping_variables_to_analyze, colnames(plot_data))
}

# Create any requested combined grouping variables, but only when all component columns exist.
created_combined_grouping_columns <- character()
if (length(combined_grouping_variables) > 0) {
  for (new_group_name in names(combined_grouping_variables)) {
    component_columns <- combined_grouping_variables[[new_group_name]]
    missing_components <- setdiff(component_columns, colnames(plot_data))

    if (length(missing_components) == 0) {
      plot_data <- plot_data %>%
        mutate(
          !!new_group_name := pmap_chr(
            select(., all_of(component_columns)),
            ~ paste(c(...), collapse = "_")
          )
        )
      created_combined_grouping_columns <- c(created_combined_grouping_columns, new_group_name)
    } else {
      message(
        "Skipping combined grouping variable '", new_group_name,
        "' because these metadata column(s) were not found: ",
        paste(missing_components, collapse = ", ")
      )
    }
  }
}

# Analyze combined variables first, followed by individual metadata variables.
grouping_columns <- c(created_combined_grouping_columns, metadata_grouping_columns)
grouping_columns <- unique(setdiff(grouping_columns, technical_columns))

# Keep grouping columns that have at least two unique non-missing values.
grouping_columns <- grouping_columns[map_int(plot_data[grouping_columns], ~ n_distinct(.x, na.rm = TRUE)) >= 2]

cat("Grouping variables that will be plotted/analyzed:\n")
Grouping variables that will be plotted/analyzed:
if (length(grouping_columns) == 0) {
  cat("None found. Add metadata columns beyond SampleID, Cell Position 1, Cell Position 2, and Plate Number, or check grouping_variables_to_analyze.\n")
} else {
  print(grouping_columns)
}
[1] "Ploidy"                                                  
[2] "Dilution Factor"                                         
[3] "SampleTissue_MassorVolume_ForCortisolExtraction_mg_or_ul"
[4] "SampleTime"                                              
[5] "Treatment"                                               
# Create compact letter displays for a grouping variable.
# Interpretation: groups that share at least one letter are not significantly different.
# Groups with no letters in common are significantly different at the selected alpha level.
make_compact_letters <- function(data, group_col, method = "tukey", alpha = 0.05) {
  # Build a clean data frame for the selected grouping variable.
  # Important: group names can sometimes contain spaces, hyphens, slashes, or other
  # characters that confuse multcompView::multcompLetters(), because that function
  # expects comparison names formatted like "GroupA-GroupB". To avoid this, the
  # statistical model uses safe temporary group names, then maps the letters back
  # onto the original metadata group names for plotting/output.
  stats_data <- data %>%
    filter(!is.na(.data[[group_col]]), !is.na(FinalConcentration_ng_mL)) %>%
    mutate(GroupOriginal = as.character(.data[[group_col]])) %>%
    droplevels()

  if (n_distinct(stats_data$GroupOriginal) < 2 || nrow(stats_data) < 3) {
    return(tibble())
  }

  original_levels <- sort(unique(stats_data$GroupOriginal))
  safe_levels <- make.names(original_levels, unique = TRUE)
  original_to_safe <- setNames(safe_levels, original_levels)
  safe_to_original <- setNames(original_levels, safe_levels)

  stats_data <- stats_data %>%
    mutate(GroupSafe = factor(original_to_safe[GroupOriginal], levels = safe_levels))

  method <- str_to_lower(method)

  p_values <- tryCatch({
    if (method == "tukey") {
      anova_model <- aov(FinalConcentration_ng_mL ~ GroupSafe, data = stats_data)
      tukey_result <- TukeyHSD(anova_model)$GroupSafe
      p <- tukey_result[, "p adj"]
      names(p) <- rownames(tukey_result)
      p
    } else if (method == "wilcoxon") {
      wilcox_result <- pairwise.wilcox.test(
        x = stats_data$FinalConcentration_ng_mL,
        g = stats_data$GroupSafe,
        p.adjust.method = "BH",
        exact = FALSE
      )
      p_table <- as.data.frame(as.table(wilcox_result$p.value), stringsAsFactors = FALSE) %>%
        filter(!is.na(Freq))

      p <- p_table$Freq
      names(p) <- paste(p_table$Var1, p_table$Var2, sep = "-")
      p
    } else {
      stop("significance_letter_method must be either 'tukey' or 'wilcoxon'.")
    }
  }, error = function(e) {
    warning("Could not calculate post-hoc p-values for compact letters for ", group_col, ": ", e$message)
    NULL
  })

  if (is.null(p_values) || length(p_values) == 0) {
    return(tibble())
  }

  letters <- tryCatch({
    multcompView::multcompLetters(p_values, threshold = alpha)$Letters
  }, error = function(e) {
    warning(
      "Post-hoc tests were calculated for ", group_col,
      ", but compact letters could not be generated. This usually happens when group names ",
      "or pairwise-comparison names are difficult to parse. The post-hoc CSV outputs are still valid. ",
      "Original error: ", e$message
    )
    NULL
  })

  if (is.null(letters) || length(letters) == 0) {
    return(tibble())
  }

  tibble(
    GroupingVariable = group_col,
    Group = unname(safe_to_original[names(letters)]),
    SignificanceLetters = unname(letters),
    LetterMethod = method,
    Alpha = alpha
  )
}

# Pre-calculate letters for every grouping variable so they can be printed, saved, and added to plots.
if (show_significance_letters && length(grouping_columns) > 0) {
  significance_letters <- map_dfr(
    grouping_columns,
    ~ make_compact_letters(
      data = plot_data,
      group_col = .x,
      method = significance_letter_method,
      alpha = posthoc_significance_alpha
    )
  )

  if (nrow(significance_letters) > 0) {
    write_csv(significance_letters, significance_letters_output)
    cat("\nCompact letter displays for plots:\n")
    print(significance_letters)
    cat("Significance letters written to: ", significance_letters_output, "\n", sep = "")
  } else {
    cat("\nNo compact letter displays could be calculated for the available grouping variables.\n")
  }
} else {
  significance_letters <- tibble()
}

Compact letter displays for plots:
# A tibble: 38 × 5
   GroupingVariable Group               SignificanceLetters LetterMethod Alpha
   <chr>            <chr>               <chr>               <chr>        <dbl>
 1 Ploidy           Triploid_Homogenate a                   wilcoxon      0.05
 2 Ploidy           Diploid_Homogenate  a                   wilcoxon      0.05
 3 Dilution Factor  25                  a                   wilcoxon      0.05
 4 Dilution Factor  26.5                a                   wilcoxon      0.05
 5 Dilution Factor  28                  a                   wilcoxon      0.05
 6 Dilution Factor  29                  a                   wilcoxon      0.05
 7 Dilution Factor  30                  a                   wilcoxon      0.05
 8 Dilution Factor  34                  a                   wilcoxon      0.05
 9 Dilution Factor  35                  a                   wilcoxon      0.05
10 Dilution Factor  45                  a                   wilcoxon      0.05
# ℹ 28 more rows
Significance letters written to: cortisol_significance_letters.csv
make_grouped_concentration_plot <- function(data, group_col) {
  # Build a plotting-specific data frame so the x-axis values and the letter-table
  # values are guaranteed to use the same variable name and type. This avoids a
  # common issue where significance letters are calculated correctly but fail to
  # appear because the letter table does not join cleanly to the plotted factor.
  plot_df <- data %>%
    filter(!is.na(.data[[group_col]]), !is.na(FinalConcentration_ng_mL)) %>%
    mutate(PlotGroup = as.factor(as.character(.data[[group_col]]))) %>%
    droplevels()

  p <- ggplot(
    plot_df,
    aes(
      x = PlotGroup,
      y = FinalConcentration_ng_mL,
      color = PlotGroup
    )
  ) +
    geom_boxplot(outlier.shape = NA, alpha = 0.35) +
    geom_point(
      position = position_jitter(width = 0.15, height = 0),
      alpha = 0.8,
      size = 2
    ) +
    labs(
      title = paste("Final cortisol concentration by", group_col),
      subtitle = if (show_significance_letters) {
        paste0(
          "Letters show ", significance_letter_method,
          " post-hoc groupings. Groups sharing a letter are not significantly different."
        )
      } else {
        NULL
      },
      x = group_col,
      y = "Final cortisol concentration (ng/mL)",
      color = group_col
    ) +
    theme_bw() +
    theme(
      axis.text.x = element_text(angle = 45, hjust = 1),
      plot.margin = margin(t = 10, r = 10, b = 10, l = 10)
    )

  if (show_significance_letters && exists("significance_letters") && nrow(significance_letters) > 0) {
    letters_for_plot <- significance_letters %>%
      filter(GroupingVariable == group_col) %>%
      transmute(
        PlotGroup = as.factor(as.character(Group)),
        SignificanceLetters = SignificanceLetters
      )

    if (nrow(letters_for_plot) > 0) {
      label_positions <- plot_df %>%
        group_by(PlotGroup) %>%
        summarize(group_max = max(FinalConcentration_ng_mL, na.rm = TRUE), .groups = "drop")

      y_range <- range(plot_df$FinalConcentration_ng_mL, na.rm = TRUE)
      y_padding <- ifelse(diff(y_range) == 0, abs(max(y_range, na.rm = TRUE)) * 0.1, diff(y_range) * 0.10)
      if (is.na(y_padding) || y_padding == 0) y_padding <- 0.1

      letters_for_plot <- letters_for_plot %>%
        left_join(label_positions, by = "PlotGroup") %>%
        filter(!is.na(group_max)) %>%
        mutate(label_y = if (use_log10_y_axis) group_max * 1.25 else group_max + y_padding)

      if (nrow(letters_for_plot) > 0) {
        # Add the letters and explicitly expand the y-axis so labels are not clipped.
        p <- p +
          geom_text(
            data = letters_for_plot,
            aes(x = PlotGroup, y = label_y, label = SignificanceLetters),
            inherit.aes = FALSE,
            color = "black",
            size = 5,
            fontface = "bold",
            vjust = 0
          ) +
          coord_cartesian(clip = "off")

        if (!use_log10_y_axis) {
          p <- p + expand_limits(y = max(letters_for_plot$label_y, na.rm = TRUE) + y_padding)
        }
      } else {
        message("Significance letters were calculated for ", group_col,
                ", but none matched the plotted groups. Check group names in the metadata.")
      }
    } else {
      message("No significance letters available for plotting variable: ", group_col)
    }
  }

  if (use_log10_y_axis) {
    p <- p +
      scale_y_log10(expand = expansion(mult = c(0.05, 0.25))) +
      labs(y = "Final cortisol concentration (ng/mL), log10 scale")
  } else {
    p <- p + scale_y_continuous(expand = expansion(mult = c(0.05, 0.20)))
  }

  p
}

if (length(grouping_columns) == 0) {
  cat("No grouping columns with at least two groups were found. No grouped boxplots were created.\n")
} else {
  pdf(grouped_boxplots_output, width = 9, height = 6)
  
  for (group_col in grouping_columns) {
    print(make_grouped_concentration_plot(plot_data, group_col))
  }
  
  dev.off()
  cat("Grouped boxplots written to: ", grouped_boxplots_output, "\n", sep = "")
}
Grouped boxplots written to: cortisol_grouped_boxplots.pdf
# Extreme outlier screening for QC review.
#
# This section is intentionally separate from the primary analysis. The primary
# figures/statistics above still include all processed samples unless you manually
# exclude SampleIDs in the user-input section. This chunk only identifies samples
# that look unusually high/low relative to their group and, if requested, creates
# a secondary exploratory data set with those rows removed.

make_outlier_key <- function(data) {
  key_cols <- c("Plate Number", "SampleID", "Cell Position 1", "Cell Position 2")
  key_cols <- key_cols[key_cols %in% names(data)]
  if (length(key_cols) == 0) {
    return(as.character(seq_len(nrow(data))))
  }
  do.call(paste, c(data[key_cols], sep = "||"))
}

screen_extreme_outliers_for_group <- function(data, group_col, iqr_multiplier = 3, min_group_n = 6) {
  if (!(group_col %in% names(data))) return(tibble())

  data %>%
    mutate(
      OutlierKey = make_outlier_key(.),
      GroupingVariable = group_col,
      Group = as.character(.data[[group_col]])
    ) %>%
    filter(!is.na(Group), !is.na(FinalConcentration_ng_mL)) %>%
    group_by(GroupingVariable, Group) %>%
    mutate(
      GroupN = n(),
      Q1 = quantile(FinalConcentration_ng_mL, 0.25, na.rm = TRUE),
      Q3 = quantile(FinalConcentration_ng_mL, 0.75, na.rm = TRUE),
      IQR_Value = IQR(FinalConcentration_ng_mL, na.rm = TRUE),
      LowerExtremeOutlierCutoff = Q1 - iqr_multiplier * IQR_Value,
      UpperExtremeOutlierCutoff = Q3 + iqr_multiplier * IQR_Value,
      IsExtremeOutlier = GroupN >= min_group_n & IQR_Value > 0 &
        (FinalConcentration_ng_mL < LowerExtremeOutlierCutoff |
           FinalConcentration_ng_mL > UpperExtremeOutlierCutoff),
      OutlierDirection = case_when(
        IsExtremeOutlier & FinalConcentration_ng_mL < LowerExtremeOutlierCutoff ~ "Low extreme outlier",
        IsExtremeOutlier & FinalConcentration_ng_mL > UpperExtremeOutlierCutoff ~ "High extreme outlier",
        TRUE ~ NA_character_
      )
    ) %>%
    ungroup() %>%
    filter(IsExtremeOutlier) %>%
    select(
      OutlierKey,
      GroupingVariable, Group, GroupN,
      any_of(c("Plate Number", "SampleID", biological_individual_id_column, "Tissue", "HatcheryLocation")),
      FinalConcentration_ng_mL,
      Q1, Q3, IQR_Value, LowerExtremeOutlierCutoff, UpperExtremeOutlierCutoff,
      OutlierDirection,
      any_of(c("ConcentrationBackCalculationStatus", "ConcentrationCaution", "Percent_B_B0", "Percent_B_B0_Warning", "ExcludedWellNote"))
    ) %>%
    distinct()
}

if (isTRUE(screen_extreme_outliers) && exists("plot_data") && nrow(plot_data) > 0 && length(grouping_columns) > 0) {
  extreme_outlier_screening <- map_dfr(
    grouping_columns,
    ~ screen_extreme_outliers_for_group(
      data = plot_data,
      group_col = .x,
      iqr_multiplier = extreme_outlier_iqr_multiplier,
      min_group_n = minimum_group_n_for_outlier_screening
    )
  )

  write_csv(extreme_outlier_screening, outlier_screening_output)

  cat("\nExtreme outlier screening used the within-group Q1 - ", extreme_outlier_iqr_multiplier,
      "*IQR and Q3 + ", extreme_outlier_iqr_multiplier,
      "*IQR rule. This is a QC screen, not an automatic reason to delete data.\n", sep = "")

  if (nrow(extreme_outlier_screening) == 0) {
    cat("No extreme outliers were flagged using the current settings.\n")
  } else {
    cat("\nPotential extreme outliers flagged for researcher review:\n")
    print(extreme_outlier_screening)
    cat("\nExtreme outlier screening table written to: ", outlier_screening_output, "\n", sep = "")
  }

  if (isTRUE(produce_secondary_outputs_without_extreme_outliers) && nrow(extreme_outlier_screening) > 0) {
    outlier_keys_to_exclude <- if (isTRUE(exclude_if_extreme_outlier_in_any_grouping)) {
      unique(extreme_outlier_screening$OutlierKey)
    } else {
      character(0)
    }

    plot_data_no_extreme_outliers <- plot_data %>%
      mutate(OutlierKey = make_outlier_key(.)) %>%
      filter(!(OutlierKey %in% outlier_keys_to_exclude)) %>%
      select(-OutlierKey)

    cat("\nSecondary exploratory data set created with ", length(outlier_keys_to_exclude),
        " row(s) removed because they were flagged as extreme outliers.\n", sep = "")
  } else {
    plot_data_no_extreme_outliers <- plot_data
  }
} else {
  extreme_outlier_screening <- tibble()
  write_csv(extreme_outlier_screening, outlier_screening_output)
  plot_data_no_extreme_outliers <- plot_data
  cat("Extreme outlier screening was skipped.\n")
}

Extreme outlier screening used the within-group Q1 - 3*IQR and Q3 + 3*IQR rule. This is a QC screen, not an automatic reason to delete data.

Potential extreme outliers flagged for researcher review:
# A tibble: 7 × 17
  OutlierKey        GroupingVariable Group GroupN `Plate Number` SampleID Tissue
  <chr>             <chr>            <chr>  <int> <chr>          <chr>    <chr> 
1 1||T-48_0610_P2|… Ploidy           Trip…     20 1              T-48_06… Homog…
2 1||T-48_0610_P4|… Ploidy           Trip…     20 1              T-48_06… Homog…
3 1||T-48_0610_P2|… SampleTime       48        15 1              T-48_06… Homog…
4 1||T-48_0610_P4|… SampleTime       48        15 1              T-48_06… Homog…
5 1||D48-RECOVER24… SampleTime       24         7 1              D48-REC… Homog…
6 1||T-48_0610_P4|… Treatment        EXP       15 1              T-48_06… Homog…
7 1||D48-RECOVER24… Treatment        REC        7 1              D48-REC… Homog…
# ℹ 10 more variables: FinalConcentration_ng_mL <dbl>, Q1 <dbl>, Q3 <dbl>,
#   IQR_Value <dbl>, LowerExtremeOutlierCutoff <dbl>,
#   UpperExtremeOutlierCutoff <dbl>, OutlierDirection <chr>,
#   Percent_B_B0 <dbl>, Percent_B_B0_Warning <chr>, ExcludedWellNote <chr>

Extreme outlier screening table written to: cortisol_extreme_outlier_screening.csv

Secondary exploratory data set created with 3 row(s) removed because they were flagged as extreme outliers.

22 Secondary exploratory figures with extreme outliers excluded

If extreme outliers were flagged above, this section creates a second set of jittered boxplots with those rows removed. These are exploratory QC plots and should be interpreted alongside the primary plots, not as a replacement unless the researcher decides the flagged samples should be excluded.

if (isTRUE(produce_secondary_outputs_without_extreme_outliers) &&
    exists("plot_data_no_extreme_outliers") &&
    exists("extreme_outlier_screening") &&
    nrow(extreme_outlier_screening) > 0 &&
    length(grouping_columns) > 0) {

  pdf(grouped_boxplots_extreme_outliers_removed_output, width = 9, height = 6)

  for (group_col in grouping_columns) {
    p_no_outliers <- make_grouped_concentration_plot(plot_data_no_extreme_outliers, group_col) +
      labs(
        title = paste("Final cortisol concentration by", group_col, "(extreme outliers excluded)"),
        subtitle = "Exploratory QC version: rows flagged by the extreme outlier screen were removed."
      )
    print(p_no_outliers)
  }

  dev.off()
  cat("Secondary grouped boxplots with extreme outliers excluded written to: ",
      grouped_boxplots_extreme_outliers_removed_output, "\n", sep = "")
} else {
  cat("No secondary outlier-excluded grouped boxplots were created because no extreme outliers were flagged or the option is turned off.\n")
}
Secondary grouped boxplots with extreme outliers excluded written to: cortisol_grouped_boxplots_extreme_outliers_removed.pdf

23 Diagnostic sample counts before statistics

This table helps troubleshoot unexpected group sizes. It reports how many non-blank, non-standard research-sample rows are present in the metadata, how many have a calculated final concentration, and how many unique research/sample IDs are represented. If group sizes look wrong, check this table first.

The diagnostic files have been written to the output folder. Check sample_processing_diagnostics.csv, samples_excluded_from_statistics_diagnostic.csv, and group_representation_diagnostics.csv if group sizes look unexpected.

24 Statistical comparisons across groups

This section runs the statistical comparisons you described: an overall test comparing all groups against one another, followed by post-hoc pairwise tests to identify which groups differ.

For each eligible metadata grouping variable, including any combined grouping variables requested in combined_grouping_variables, the code outputs:

  1. Group summary statistics: n, mean, SD, median, and IQR.
  2. Parametric overall test: one-way ANOVA.
  3. Parametric post-hoc test: Tukey HSD.
  4. Non-parametric overall test: Kruskal-Wallis.
  5. Non-parametric post-hoc test: pairwise Wilcoxon tests with Benjamini-Hochberg adjusted p-values.
  6. Quick-reference tables showing only significant post-hoc comparisons.

The code also prints a simple assumption check. If the ANOVA residuals look non-normal or the group variances look unequal, the Kruskal-Wallis/Wilcoxon results are usually safer to emphasize, especially with small sample sizes.

# This section is intentionally defensive because metadata-driven grouping variables
# can vary a lot from one run to the next. Some grouping variables may have too few
# samples, only one replicate per group, all-identical values, or otherwise fail one
# of the statistical tests. Those cases are now skipped with an explanatory note
# instead of stopping the whole R Markdown render.

# Helper: safely convert a value to TRUE/FALSE.
is_true <- function(x) {
  is.logical(x) && length(x) == 1 && !is.na(x) && x
}

# Helper: format p-values for assumption notes without creating "NA" warnings.
format_p <- function(x) {
  if (length(x) != 1 || is.na(x)) {
    return("not available")
  }
  as.character(round(x, 4))
}

# Helper to format pairwise.wilcox.test output into a tidy table.
tidy_pairwise_wilcox <- function(pairwise_result, grouping_variable) {
  if (is.null(pairwise_result) || is.null(pairwise_result$p.value)) {
    return(tibble())
  }

  p_matrix <- pairwise_result$p.value

  as.data.frame(as.table(p_matrix), stringsAsFactors = FALSE) %>%
    as_tibble() %>%
    filter(!is.na(Freq)) %>%
    transmute(
      GroupingVariable = grouping_variable,
      group1 = as.character(Var1),
      group2 = as.character(Var2),
      p_adjusted = as.numeric(Freq),
      p_adjust_method = pairwise_result$p.adjust.method
    )
}

# Helper to format TukeyHSD output into a tidy table.
tidy_tukey <- function(tukey_result, grouping_variable) {
  if (is.null(tukey_result) || is.null(tukey_result$Group)) {
    return(tibble())
  }

  tukey_table <- as.data.frame(tukey_result$Group)
  tukey_table$comparison <- rownames(tukey_table)
  rownames(tukey_table) <- NULL

  tukey_table %>%
    as_tibble() %>%
    separate(comparison, into = c("group1", "group2"), sep = "-", remove = FALSE, extra = "merge") %>%
    transmute(
      GroupingVariable = grouping_variable,
      comparison,
      group1,
      group2,
      mean_difference = diff,
      conf_low = lwr,
      conf_high = upr,
      p_adjusted = `p adj`
    )
}

# Helper: create a standardized "not run" row.
not_run_row <- function(group_col, test_name, note) {
  tibble(
    GroupingVariable = group_col,
    Test = test_name,
    Statistic = NA_real_,
    df = NA_real_,
    p_value = NA_real_,
    AssumptionNotes = note,
    RecommendedInterpretation = "Do not interpret this test for this grouping variable."
  )
}

if (is_true(run_group_statistics) && length(grouping_columns) > 0) {
  group_summary_statistics <- list()
  overall_test_results <- list()
  tukey_posthoc_results <- list()
  wilcoxon_posthoc_results <- list()

  for (group_col in grouping_columns) {
    stats_data <- plot_data %>%
      filter(!is.na(.data[[group_col]]), !is.na(FinalConcentration_ng_mL)) %>%
      mutate(Group = as.factor(.data[[group_col]])) %>%
      droplevels()

    number_of_groups <- n_distinct(stats_data$Group)
    total_n <- nrow(stats_data)

    # Summary statistics for each group.
    # n_rows_used_in_stats counts concentration rows used in the tests.
    # n_unique_ids_used_in_stats counts unique SampleID values by default, or the
    # column named in statistics_unique_id_column if you set one. This prevents
    # confusing row counts with biological/sample counts.
    this_group_summary <- stats_data %>%
      group_by(Group) %>%
      summarize(
        n_rows_used_in_stats = n(),
        n_unique_ids_used_in_stats = n_distinct(.data[[statistics_unique_id_column]], na.rm = TRUE),
        mean_ng_mL = mean(FinalConcentration_ng_mL, na.rm = TRUE),
        sd_ng_mL = sd(FinalConcentration_ng_mL, na.rm = TRUE),
        median_ng_mL = median(FinalConcentration_ng_mL, na.rm = TRUE),
        iqr_ng_mL = IQR(FinalConcentration_ng_mL, na.rm = TRUE),
        .groups = "drop"
      ) %>%
      mutate(GroupingVariable = group_col) %>%
      select(GroupingVariable, Group, n_rows_used_in_stats, n_unique_ids_used_in_stats, mean_ng_mL, sd_ng_mL, median_ng_mL, iqr_ng_mL)

    group_summary_statistics[[group_col]] <- this_group_summary

    # For ANOVA/Tukey/Kruskal/Wilcoxon, at least 2 groups and at least 3 total samples are needed.
    # Tukey and ANOVA also require residual degrees of freedom, which fails when every group has n = 1.
    group_counts <- stats_data %>% count(Group, name = "n")
    residual_df <- total_n - number_of_groups

    if (number_of_groups < 2 || total_n < 3) {
      note <- "Fewer than 2 groups or too few total samples."
      overall_test_results[[paste0(group_col, "_NotRun")]] <- not_run_row(group_col, "Not run", note)
      next
    }

    if (residual_df < 1) {
      note <- "No residual degrees of freedom because each group has too few observations, often n = 1 per group."
      overall_test_results[[paste0(group_col, "_ANOVA_NotRun")]] <- not_run_row(group_col, "One-way ANOVA", note)
      overall_test_results[[paste0(group_col, "_Kruskal_NotRun")]] <- not_run_row(group_col, "Kruskal-Wallis", note)
      next
    }

    # One-way ANOVA across all groups.
    anova_model <- tryCatch(
      aov(FinalConcentration_ng_mL ~ Group, data = stats_data),
      error = function(e) {
        message("ANOVA skipped for ", group_col, ": ", e$message)
        NULL
      }
    )

    if (!is.null(anova_model)) {
      anova_table <- summary(anova_model)[[1]]
      anova_p <- anova_table[["Pr(>F)"]][1]
      anova_f <- anova_table[["F value"]][1]
      anova_df <- anova_table[["Df"]][1]

      residuals_anova <- residuals(anova_model)
      shapiro_p <- tryCatch(
        if (length(residuals_anova) >= 3 && length(residuals_anova) <= 5000 && sd(residuals_anova, na.rm = TRUE) > 0) {
          shapiro.test(residuals_anova)$p.value
        } else {
          NA_real_
        },
        error = function(e) NA_real_
      )

      bartlett_p <- tryCatch(
        bartlett.test(FinalConcentration_ng_mL ~ Group, data = stats_data)$p.value,
        error = function(e) NA_real_
      )

      assumption_notes <- paste0(
        "ANOVA residual Shapiro-Wilk p = ", format_p(shapiro_p),
        "; Bartlett equal-variance p = ", format_p(bartlett_p), "."
      )

      recommended_interpretation <- if (
        (!is.na(shapiro_p) && shapiro_p < 0.05) ||
          (!is.na(bartlett_p) && bartlett_p < 0.05)
      ) {
        "Assumptions may be weak; emphasize Kruskal-Wallis and pairwise Wilcoxon results."
      } else {
        "ANOVA assumptions are not strongly contradicted; ANOVA/Tukey are reasonable, but still review sample size and plots."
      }

      overall_test_results[[paste0(group_col, "_ANOVA")]] <- tibble(
        GroupingVariable = group_col,
        Test = "One-way ANOVA",
        Statistic = anova_f,
        df = anova_df,
        p_value = anova_p,
        AssumptionNotes = assumption_notes,
        RecommendedInterpretation = recommended_interpretation
      )

      # Tukey post-hoc pairwise comparisons for ANOVA.
      # This can fail for some metadata groupings, so it is wrapped in tryCatch.
      tukey_result <- tryCatch(
        {
          tr <- TukeyHSD(anova_model)
          names(tr) <- "Group"
          tr
        },
        error = function(e) {
          message("Tukey post-hoc skipped for ", group_col, ": ", e$message)
          NULL
        }
      )
      tukey_posthoc_results[[group_col]] <- tidy_tukey(tukey_result, group_col)
    } else {
      overall_test_results[[paste0(group_col, "_ANOVA_NotRun")]] <- not_run_row(
        group_col,
        "One-way ANOVA",
        "ANOVA model could not be fit for this grouping variable."
      )
    }

    # Kruskal-Wallis non-parametric overall test across all groups.
    kruskal_result <- tryCatch(
      kruskal.test(FinalConcentration_ng_mL ~ Group, data = stats_data),
      error = function(e) {
        message("Kruskal-Wallis skipped for ", group_col, ": ", e$message)
        NULL
      }
    )

    if (!is.null(kruskal_result)) {
      overall_test_results[[paste0(group_col, "_Kruskal")]] <- tibble(
        GroupingVariable = group_col,
        Test = "Kruskal-Wallis",
        Statistic = unname(kruskal_result$statistic),
        df = unname(kruskal_result$parameter),
        p_value = kruskal_result$p.value,
        AssumptionNotes = "Non-parametric comparison of group distributions; useful when ANOVA assumptions are questionable.",
        RecommendedInterpretation = "Use with the pairwise Wilcoxon table to identify which groups differ."
      )
    } else {
      overall_test_results[[paste0(group_col, "_Kruskal_NotRun")]] <- not_run_row(
        group_col,
        "Kruskal-Wallis",
        "Kruskal-Wallis test could not be run for this grouping variable."
      )
    }

    # Pairwise Wilcoxon post-hoc comparisons with multiple-test correction.
    wilcox_result <- tryCatch(
      pairwise.wilcox.test(
        x = stats_data$FinalConcentration_ng_mL,
        g = stats_data$Group,
        p.adjust.method = "BH",
        exact = FALSE
      ),
      error = function(e) {
        message("Pairwise Wilcoxon skipped for ", group_col, ": ", e$message)
        NULL
      }
    )
    wilcoxon_posthoc_results[[group_col]] <- tidy_pairwise_wilcox(wilcox_result, group_col)
  }

  group_statistics <- bind_rows(group_summary_statistics)
  overall_statistics <- bind_rows(overall_test_results)
  tukey_posthoc_statistics <- bind_rows(tukey_posthoc_results)
  wilcoxon_posthoc_statistics <- bind_rows(wilcoxon_posthoc_results)

  # Ensure empty post-hoc objects still have the expected columns so downstream filtering/writing does not fail.
  if (nrow(tukey_posthoc_statistics) == 0) {
    tukey_posthoc_statistics <- tibble(
      GroupingVariable = character(), comparison = character(), group1 = character(), group2 = character(),
      mean_difference = numeric(), conf_low = numeric(), conf_high = numeric(), p_adjusted = numeric()
    )
  }

  if (nrow(wilcoxon_posthoc_statistics) == 0) {
    wilcoxon_posthoc_statistics <- tibble(
      GroupingVariable = character(), group1 = character(), group2 = character(),
      p_adjusted = numeric(), p_adjust_method = character()
    )
  }

  cat("\nGroup summary statistics:\n")
  print(group_statistics)

  cat("\nOverall group comparison tests. These p-values test whether there is evidence of any difference among groups:\n")
  print(overall_statistics)

  cat("\nTukey post-hoc pairwise comparisons. These adjusted p-values identify which specific groups differ after ANOVA:\n")
  print(tukey_posthoc_statistics)

  cat("\nPairwise Wilcoxon post-hoc comparisons. These adjusted p-values identify which specific groups differ using a non-parametric approach:\n")
  print(wilcoxon_posthoc_statistics)

  significant_tukey_posthoc <- tukey_posthoc_statistics %>%
    filter(!is.na(p_adjusted), p_adjusted < posthoc_significance_alpha) %>%
    arrange(GroupingVariable, p_adjusted)

  significant_wilcoxon_posthoc <- wilcoxon_posthoc_statistics %>%
    filter(!is.na(p_adjusted), p_adjusted < posthoc_significance_alpha) %>%
    arrange(GroupingVariable, p_adjusted)

  cat("\nQuick reference: significant Tukey post-hoc comparisons only ",
      "(adjusted p < ", posthoc_significance_alpha, "):\n", sep = "")
  if (nrow(significant_tukey_posthoc) == 0) {
    cat("No significant Tukey post-hoc comparisons found at this threshold.\n")
  } else {
    print(significant_tukey_posthoc)
  }

  cat("\nQuick reference: significant pairwise Wilcoxon post-hoc comparisons only ",
      "(BH-adjusted p < ", posthoc_significance_alpha, "):\n", sep = "")
  if (nrow(significant_wilcoxon_posthoc) == 0) {
    cat("No significant pairwise Wilcoxon post-hoc comparisons found at this threshold.\n")
  } else {
    print(significant_wilcoxon_posthoc)
  }

  write_csv(group_statistics, statistical_summary_output)
  write_csv(overall_statistics, anova_kruskal_summary_output)
  write_csv(tukey_posthoc_statistics, posthoc_tukey_output)
  write_csv(wilcoxon_posthoc_statistics, posthoc_wilcoxon_output)
  write_csv(significant_tukey_posthoc, significant_tukey_output)
  write_csv(significant_wilcoxon_posthoc, significant_wilcoxon_output)

  cat("\nGroup summary statistics written to: ", statistical_summary_output, "\n", sep = "")
  cat("Overall ANOVA/Kruskal-Wallis results written to: ", anova_kruskal_summary_output, "\n", sep = "")
  cat("Tukey post-hoc results written to: ", posthoc_tukey_output, "\n", sep = "")
  cat("Pairwise Wilcoxon post-hoc results written to: ", posthoc_wilcoxon_output, "\n", sep = "")
  cat("Significant Tukey post-hoc quick reference written to: ", significant_tukey_output, "\n", sep = "")
  cat("Significant Wilcoxon post-hoc quick reference written to: ", significant_wilcoxon_output, "\n", sep = "")
} else {
  cat("Optional group statistics were skipped.\n")
}

Group summary statistics:
# A tibble: 38 × 8
   GroupingVariable Group n_rows_used_in_stats n_unique_ids_used_in…¹ mean_ng_mL
   <chr>            <fct>                <int>                  <int>      <dbl>
 1 Ploidy           Dipl…                   19                     19    0.0131 
 2 Ploidy           Trip…                   20                     20    0.0163 
 3 Dilution Factor  0                        9                      9    0.0248 
 4 Dilution Factor  25                       1                      1    0.00504
 5 Dilution Factor  26.5                     1                      1    0.00239
 6 Dilution Factor  28                       1                      1    0.00633
 7 Dilution Factor  29                       1                      1    0.0129 
 8 Dilution Factor  30                       1                      1    0.00143
 9 Dilution Factor  34                       1                      1    0.0123 
10 Dilution Factor  35                       1                      1    0.0118 
# ℹ 28 more rows
# ℹ abbreviated name: ¹​n_unique_ids_used_in_stats
# ℹ 3 more variables: sd_ng_mL <dbl>, median_ng_mL <dbl>, iqr_ng_mL <dbl>

Overall group comparison tests. These p-values test whether there is evidence of any difference among groups:
# A tibble: 10 × 7
   GroupingVariable                Test  Statistic    df p_value AssumptionNotes
   <chr>                           <chr>     <dbl> <dbl>   <dbl> <chr>          
 1 Ploidy                          One-…    0.400      1 5.31e-1 ANOVA residual…
 2 Ploidy                          Krus…    0.0639     1 8.00e-1 Non-parametric…
 3 Dilution Factor                 One-…    0.616     10 7.87e-1 ANOVA residual…
 4 Dilution Factor                 Krus…    7.23      10 7.04e-1 Non-parametric…
 5 SampleTissue_MassorVolume_ForC… One-…   14.0       17 7.41e-8 ANOVA residual…
 6 SampleTissue_MassorVolume_ForC… Krus…   25.7       17 7.92e-2 Non-parametric…
 7 SampleTime                      One-…    0.881      3 4.60e-1 ANOVA residual…
 8 SampleTime                      Krus…    0.816      3 8.46e-1 Non-parametric…
 9 Treatment                       One-…    0.445      2 6.44e-1 ANOVA residual…
10 Treatment                       Krus…    0.876      2 6.45e-1 Non-parametric…
# ℹ 1 more variable: RecommendedInterpretation <chr>

Tukey post-hoc pairwise comparisons. These adjusted p-values identify which specific groups differ after ANOVA:
# A tibble: 218 × 8
   GroupingVariable comparison  group1 group2 mean_difference conf_low conf_high
   <chr>            <chr>       <chr>  <chr>            <dbl>    <dbl>     <dbl>
 1 Ploidy           Triploid_H… Tripl… Diplo…         0.00310 -0.00684    0.0130
 2 Dilution Factor  25-0        25     0             -0.0198  -0.0788     0.0393
 3 Dilution Factor  26.5-0      26.5   0             -0.0224  -0.0815     0.0367
 4 Dilution Factor  28-0        28     0             -0.0185  -0.0776     0.0406
 5 Dilution Factor  29-0        29     0             -0.0119  -0.0710     0.0472
 6 Dilution Factor  30-0        30     0             -0.0234  -0.0825     0.0357
 7 Dilution Factor  34-0        34     0             -0.0125  -0.0716     0.0466
 8 Dilution Factor  35-0        35     0             -0.0130  -0.0721     0.0460
 9 Dilution Factor  45-0        45     0             -0.00850 -0.0523     0.0353
10 Dilution Factor  45.5-0      45.5   0             -0.00404 -0.0631     0.0550
# ℹ 208 more rows
# ℹ 1 more variable: p_adjusted <dbl>

Pairwise Wilcoxon post-hoc comparisons. These adjusted p-values identify which specific groups differ using a non-parametric approach:
# A tibble: 218 × 5
   GroupingVariable group1              group2        p_adjusted p_adjust_method
   <chr>            <chr>               <chr>              <dbl> <chr>          
 1 Ploidy           Triploid_Homogenate Diploid_Homo…      0.811 BH             
 2 Dilution Factor  25                  0                  1     BH             
 3 Dilution Factor  26.5                0                  1     BH             
 4 Dilution Factor  28                  0                  1     BH             
 5 Dilution Factor  29                  0                  1     BH             
 6 Dilution Factor  30                  0                  1     BH             
 7 Dilution Factor  34                  0                  1     BH             
 8 Dilution Factor  35                  0                  1     BH             
 9 Dilution Factor  45                  0                  1     BH             
10 Dilution Factor  45.5                0                  1     BH             
# ℹ 208 more rows

Quick reference: significant Tukey post-hoc comparisons only (adjusted p < 0.05):
# A tibble: 44 × 8
   GroupingVariable  comparison group1 group2 mean_difference conf_low conf_high
   <chr>             <chr>      <chr>  <chr>            <dbl>    <dbl>     <dbl>
 1 SampleTissue_Mas… 100-25     100    25             -0.0594  -0.0829   -0.0359
 2 SampleTissue_Mas… 25-10      25     10              0.0705   0.0424    0.0985
 3 SampleTissue_Mas… 40-25      40     25             -0.0717  -0.104    -0.0392
 4 SampleTissue_Mas… 43-10      43     10              0.0612   0.0331    0.0893
 5 SampleTissue_Mas… 60-25      60     25             -0.0702  -0.103    -0.0378
 6 SampleTissue_Mas… 53-25      53     25             -0.0693  -0.102    -0.0368
 7 SampleTissue_Mas… 25-11      25     11              0.0693   0.0368    0.102 
 8 SampleTissue_Mas… 100-43     100    43             -0.0501  -0.0736   -0.0266
 9 SampleTissue_Mas… 50-25      50     25             -0.0666  -0.0991   -0.0342
10 SampleTissue_Mas… 56-25      56     25             -0.0653  -0.0978   -0.0329
# ℹ 34 more rows
# ℹ 1 more variable: p_adjusted <dbl>

Quick reference: significant pairwise Wilcoxon post-hoc comparisons only (BH-adjusted p < 0.05):
No significant pairwise Wilcoxon post-hoc comparisons found at this threshold.

Group summary statistics written to: cortisol_group_statistics.csv
Overall ANOVA/Kruskal-Wallis results written to: cortisol_overall_group_tests.csv
Tukey post-hoc results written to: cortisol_tukey_posthoc_tests.csv
Pairwise Wilcoxon post-hoc results written to: cortisol_wilcoxon_posthoc_tests.csv
Significant Tukey post-hoc quick reference written to: cortisol_significant_tukey_posthoc_tests.csv
Significant Wilcoxon post-hoc quick reference written to: cortisol_significant_wilcoxon_posthoc_tests.csv
# Secondary exploratory statistics with extreme outliers excluded.
# This reuses the same grouping variables as the primary statistics section, but
# uses plot_data_no_extreme_outliers. The primary statistics are not modified.

run_group_statistics_simple <- function(data, grouping_columns, label_suffix = "") {
  group_summary_statistics <- list()
  overall_test_results <- list()
  tukey_posthoc_results <- list()
  wilcoxon_posthoc_results <- list()

  for (group_col in grouping_columns) {
    if (!(group_col %in% names(data))) next

    stats_data <- data %>%
      filter(!is.na(.data[[group_col]]), !is.na(FinalConcentration_ng_mL)) %>%
      mutate(Group = as.factor(as.character(.data[[group_col]]))) %>%
      droplevels()

    if (nrow(stats_data) == 0 || n_distinct(stats_data$Group) < 2) next

    group_summary_statistics[[group_col]] <- stats_data %>%
      group_by(Group) %>%
      summarize(
        GroupingVariable = group_col,
        n_rows_used_in_stats = n(),
        n_unique_ids_used_in_stats = if (statistics_unique_id_column %in% names(.)) {
          n_distinct(.data[[statistics_unique_id_column]], na.rm = TRUE)
        } else {
          n_distinct(SampleID, na.rm = TRUE)
        },
        mean_FinalConcentration_ng_mL = mean(FinalConcentration_ng_mL, na.rm = TRUE),
        sd_FinalConcentration_ng_mL = sd(FinalConcentration_ng_mL, na.rm = TRUE),
        median_FinalConcentration_ng_mL = median(FinalConcentration_ng_mL, na.rm = TRUE),
        IQR_FinalConcentration_ng_mL = IQR(FinalConcentration_ng_mL, na.rm = TRUE),
        .groups = "drop"
      ) %>%
      relocate(GroupingVariable, Group) %>%
      mutate(AnalysisDataset = label_suffix)

    if (nrow(stats_data) >= 3 && n_distinct(stats_data$Group) >= 2) {
      anova_result <- tryCatch({
        anova_model <- aov(FinalConcentration_ng_mL ~ Group, data = stats_data)
        anova_summary <- summary(anova_model)[[1]]
        tibble(
          GroupingVariable = group_col,
          Test = "ANOVA",
          Statistic = as.numeric(anova_summary$`F value`[1]),
          df = as.numeric(anova_summary$Df[1]),
          p_value = as.numeric(anova_summary$`Pr(>F)`[1]),
          AnalysisDataset = label_suffix
        )
      }, error = function(e) {
        tibble(GroupingVariable = group_col, Test = "ANOVA", Statistic = NA_real_, df = NA_real_, p_value = NA_real_, AnalysisDataset = label_suffix)
      })

      kruskal_result <- tryCatch({
        kt <- kruskal.test(FinalConcentration_ng_mL ~ Group, data = stats_data)
        tibble(
          GroupingVariable = group_col,
          Test = "Kruskal-Wallis",
          Statistic = as.numeric(kt$statistic),
          df = as.numeric(kt$parameter),
          p_value = as.numeric(kt$p.value),
          AnalysisDataset = label_suffix
        )
      }, error = function(e) {
        tibble(GroupingVariable = group_col, Test = "Kruskal-Wallis", Statistic = NA_real_, df = NA_real_, p_value = NA_real_, AnalysisDataset = label_suffix)
      })

      overall_test_results[[group_col]] <- bind_rows(anova_result, kruskal_result)

      tukey_posthoc_results[[group_col]] <- tryCatch({
        anova_model <- aov(FinalConcentration_ng_mL ~ Group, data = stats_data)
        tidy_tukey(TukeyHSD(anova_model), group_col) %>% mutate(AnalysisDataset = label_suffix)
      }, error = function(e) tibble())

      wilcoxon_posthoc_results[[group_col]] <- tryCatch({
        tidy_pairwise_wilcox(
          pairwise.wilcox.test(stats_data$FinalConcentration_ng_mL, stats_data$Group, p.adjust.method = "BH", exact = FALSE),
          group_col
        ) %>% mutate(AnalysisDataset = label_suffix)
      }, error = function(e) tibble())
    }
  }

  list(
    summary = bind_rows(group_summary_statistics),
    overall = bind_rows(overall_test_results),
    tukey = bind_rows(tukey_posthoc_results),
    wilcoxon = bind_rows(wilcoxon_posthoc_results)
  )
}

if (isTRUE(produce_secondary_outputs_without_extreme_outliers) &&
    exists("plot_data_no_extreme_outliers") &&
    exists("extreme_outlier_screening") &&
    nrow(extreme_outlier_screening) > 0 &&
    isTRUE(run_group_statistics) &&
    length(grouping_columns) > 0) {

  secondary_stats <- run_group_statistics_simple(
    data = plot_data_no_extreme_outliers,
    grouping_columns = grouping_columns,
    label_suffix = "Extreme outliers excluded"
  )

  secondary_significant_tukey <- secondary_stats$tukey %>%
    filter(!is.na(p_adjusted), p_adjusted < posthoc_significance_alpha) %>%
    arrange(GroupingVariable, p_adjusted)

  secondary_significant_wilcoxon <- secondary_stats$wilcoxon %>%
    filter(!is.na(p_adjusted), p_adjusted < posthoc_significance_alpha) %>%
    arrange(GroupingVariable, p_adjusted)

  cat("\nSecondary exploratory statistics with extreme outliers excluded:\n")
  cat("\nGroup summary statistics:\n")
  print(secondary_stats$summary)
  cat("\nOverall tests:\n")
  print(secondary_stats$overall)
  cat("\nSignificant Tukey post-hoc comparisons only:\n")
  if (nrow(secondary_significant_tukey) == 0) cat("No significant Tukey comparisons found.\n") else print(secondary_significant_tukey)
  cat("\nSignificant Wilcoxon post-hoc comparisons only:\n")
  if (nrow(secondary_significant_wilcoxon) == 0) cat("No significant Wilcoxon comparisons found.\n") else print(secondary_significant_wilcoxon)

  write_csv(secondary_stats$summary, statistical_summary_extreme_outliers_removed_output)
  write_csv(secondary_stats$overall, anova_kruskal_summary_extreme_outliers_removed_output)
  write_csv(secondary_stats$tukey, posthoc_tukey_extreme_outliers_removed_output)
  write_csv(secondary_stats$wilcoxon, posthoc_wilcoxon_extreme_outliers_removed_output)
  write_csv(secondary_significant_tukey, significant_tukey_extreme_outliers_removed_output)
  write_csv(secondary_significant_wilcoxon, significant_wilcoxon_extreme_outliers_removed_output)
} else {
  cat("Secondary outlier-excluded statistics were skipped because no extreme outliers were flagged, statistics are off, or the option is turned off.\n")
}

Secondary exploratory statistics with extreme outliers excluded:

Group summary statistics:
# A tibble: 35 × 9
   GroupingVariable Group            n_rows_used_in_stats n_unique_ids_used_in…¹
   <chr>            <fct>                           <int>                  <int>
 1 Ploidy           Diploid_Homogen…                   18                     18
 2 Ploidy           Triploid_Homoge…                   18                     18
 3 Dilution Factor  0                                   6                      6
 4 Dilution Factor  25                                  1                      1
 5 Dilution Factor  26.5                                1                      1
 6 Dilution Factor  28                                  1                      1
 7 Dilution Factor  29                                  1                      1
 8 Dilution Factor  30                                  1                      1
 9 Dilution Factor  34                                  1                      1
10 Dilution Factor  35                                  1                      1
# ℹ 25 more rows
# ℹ abbreviated name: ¹​n_unique_ids_used_in_stats
# ℹ 5 more variables: mean_FinalConcentration_ng_mL <dbl>,
#   sd_FinalConcentration_ng_mL <dbl>, median_FinalConcentration_ng_mL <dbl>,
#   IQR_FinalConcentration_ng_mL <dbl>, AnalysisDataset <chr>

Overall tests:
# A tibble: 10 × 6
   GroupingVariable                Test  Statistic    df p_value AnalysisDataset
   <chr>                           <chr>     <dbl> <dbl>   <dbl> <chr>          
 1 Ploidy                          ANOVA   0.168       1 0.685   Extreme outlie…
 2 Ploidy                          Krus…   0.00100     1 0.975   Extreme outlie…
 3 Dilution Factor                 ANOVA   0.649      10 0.759   Extreme outlie…
 4 Dilution Factor                 Krus…  12.6        10 0.247   Extreme outlie…
 5 SampleTissue_MassorVolume_ForC… ANOVA   3.17       14 0.00839 Extreme outlie…
 6 SampleTissue_MassorVolume_ForC… Krus…  20.7        14 0.111   Extreme outlie…
 7 SampleTime                      ANOVA   0.0779      3 0.972   Extreme outlie…
 8 SampleTime                      Krus…   0.146       3 0.986   Extreme outlie…
 9 Treatment                       ANOVA   0.382       2 0.686   Extreme outlie…
10 Treatment                       Krus…   2.60        2 0.273   Extreme outlie…

Significant Tukey post-hoc comparisons only:
# A tibble: 7 × 9
  GroupingVariable   comparison group1 group2 mean_difference conf_low conf_high
  <chr>              <chr>      <chr>  <chr>            <dbl>    <dbl>     <dbl>
1 SampleTissue_Mass… 14-10      14     10              0.0361  0.00893  0.0633  
2 SampleTissue_Mass… 40-14      40     14             -0.0373 -0.0687  -0.00594 
3 SampleTissue_Mass… 60-14      60     14             -0.0359 -0.0673  -0.00450 
4 SampleTissue_Mass… 53-14      53     14             -0.0349 -0.0663  -0.00355 
5 SampleTissue_Mass… 14-11      14     11              0.0349  0.00353  0.0663  
6 SampleTissue_Mass… 14-100     14     100             0.0250  0.00228  0.0478  
7 SampleTissue_Mass… 50-14      50     14             -0.0323 -0.0637  -0.000903
# ℹ 2 more variables: p_adjusted <dbl>, AnalysisDataset <chr>

Significant Wilcoxon post-hoc comparisons only:
No significant Wilcoxon comparisons found.

25 Within-group dispersion analysis

This section summarizes within-group variation in final cortisol concentration. It reports robust dispersion metrics (IQR and MAD), plus SD, variance, and coefficient of variation. It also runs a Fligner-Killeen test for each grouping variable to test whether dispersion differs among groups. Violin plots with boxplot overlays are generated to emphasize distribution shape and within-group spread; each group is labeled with its IQR. The same analysis is run for the full dataset and, when available, the exploratory dataset with extreme outliers removed.

# -----------------------------------------------------------------------------
# Within-group dispersion analysis
# -----------------------------------------------------------------------------
# This section is intended to complement the group-comparison statistics above.
# Kruskal-Wallis / ANOVA ask whether group central tendencies differ, whereas
# this section asks how variable the samples are within each group.

run_dispersion_analysis <- function(data, grouping_columns, dataset_label,
                                    summary_output_file, fligner_output_file,
                                    plot_output_file) {
  if (!is.data.frame(data) || nrow(data) == 0) {
    cat("Dispersion analysis skipped for ", dataset_label, ": no data rows available.\n", sep = "")
    write_csv(tibble(), summary_output_file)
    write_csv(tibble(), fligner_output_file)
    return(invisible(NULL))
  }

  if (length(grouping_columns) == 0) {
    cat("Dispersion analysis skipped for ", dataset_label, ": no grouping columns available.\n", sep = "")
    write_csv(tibble(), summary_output_file)
    write_csv(tibble(), fligner_output_file)
    return(invisible(NULL))
  }

  dispersion_data <- data %>%
    filter(!is.na(FinalConcentration_ng_mL))

  dispersion_summary <- map_dfr(grouping_columns, function(group_col) {
    if (!(group_col %in% names(dispersion_data))) return(tibble())

    dispersion_data %>%
      mutate(Group = as.character(.data[[group_col]])) %>%
      filter(!is.na(Group), Group != "") %>%
      group_by(Group) %>%
      summarize(
        AnalysisDataset = dataset_label,
        GroupingVariable = group_col,
        n = n(),
        median_ng_mL = median(FinalConcentration_ng_mL, na.rm = TRUE),
        mean_ng_mL = mean(FinalConcentration_ng_mL, na.rm = TRUE),
        sd_ng_mL = sd(FinalConcentration_ng_mL, na.rm = TRUE),
        variance_ng_mL = var(FinalConcentration_ng_mL, na.rm = TRUE),
        IQR_ng_mL = IQR(FinalConcentration_ng_mL, na.rm = TRUE),
        MAD_ng_mL = mad(FinalConcentration_ng_mL, na.rm = TRUE),
        min_ng_mL = min(FinalConcentration_ng_mL, na.rm = TRUE),
        max_ng_mL = max(FinalConcentration_ng_mL, na.rm = TRUE),
        coefficient_of_variation_percent = ifelse(
          is.na(mean_ng_mL) | mean_ng_mL == 0,
          NA_real_,
          100 * sd_ng_mL / mean_ng_mL
        ),
        .groups = "drop"
      ) %>%
      select(AnalysisDataset, GroupingVariable, Group, everything())
  })

  fligner_results <- map_dfr(grouping_columns, function(group_col) {
    if (!(group_col %in% names(dispersion_data))) return(tibble())

    test_data <- dispersion_data %>%
      mutate(Group = as.factor(.data[[group_col]])) %>%
      filter(!is.na(Group), !is.na(FinalConcentration_ng_mL)) %>%
      group_by(Group) %>%
      filter(n() >= 2) %>%
      ungroup()

    if (n_distinct(test_data$Group) < 2 || nrow(test_data) < 3) {
      return(tibble(
        AnalysisDataset = dataset_label,
        GroupingVariable = group_col,
        Test = "Fligner-Killeen test",
        Statistic = NA_real_,
        P_Value = NA_real_,
        Note = "Skipped: fewer than two groups with at least two observations"
      ))
    }

    tryCatch({
      fligner_result <- fligner.test(FinalConcentration_ng_mL ~ Group, data = test_data)
      tibble(
        AnalysisDataset = dataset_label,
        GroupingVariable = group_col,
        Test = "Fligner-Killeen test",
        Statistic = unname(fligner_result$statistic),
        P_Value = fligner_result$p.value,
        Note = ifelse(fligner_result$p.value < 0.05,
                      "Significant evidence that dispersion differs among groups",
                      "No significant evidence that dispersion differs among groups")
      )
    }, error = function(e) {
      tibble(
        AnalysisDataset = dataset_label,
        GroupingVariable = group_col,
        Test = "Fligner-Killeen test",
        Statistic = NA_real_,
        P_Value = NA_real_,
        Note = paste("Skipped/error:", e$message)
      )
    })
  })

  write_csv(dispersion_summary, summary_output_file)
  write_csv(fligner_results, fligner_output_file)

  cat("\nDispersion summary written to: ", summary_output_file, "\n", sep = "")
  cat("Fligner-Killeen dispersion test results written to: ", fligner_output_file, "\n", sep = "")
  cat("\nDispersion summary for ", dataset_label, ":\n", sep = "")
  print(dispersion_summary)
  cat("\nFligner-Killeen tests for ", dataset_label, ":\n", sep = "")
  print(fligner_results)

  # Visualize within-group dispersion for each grouping variable.
  pdf(plot_output_file, width = 9, height = 6)
  for (group_col in grouping_columns) {
    if (!(group_col %in% names(dispersion_data))) next

    plot_df <- dispersion_data %>%
      mutate(Group = as.factor(.data[[group_col]])) %>%
      filter(!is.na(Group), !is.na(FinalConcentration_ng_mL))

    if (n_distinct(plot_df$Group) < 2) next

    # Pull IQR values from the summary table for this grouping variable so they
    # can be displayed directly on the dispersion plots.
    label_df <- dispersion_summary %>%
      filter(
        AnalysisDataset == dataset_label,
        GroupingVariable == group_col,
        Group %in% as.character(plot_df$Group)
      ) %>%
      mutate(
        Group = factor(Group, levels = levels(plot_df$Group)),
        IQR_Label = paste0("IQR = ", signif(IQR_ng_mL, 3))
      )

    # Put labels slightly above the observed data for each group. This is done
    # group-by-group so labels remain readable even if some groups have much
    # larger concentrations than others.
    label_positions <- plot_df %>%
      group_by(Group) %>%
      summarize(
        LabelY = max(FinalConcentration_ng_mL, na.rm = TRUE),
        .groups = "drop"
      ) %>%
      mutate(
        LabelY = ifelse(is.finite(LabelY) & LabelY > 0, LabelY * 1.20, NA_real_)
      )

    label_df <- label_df %>%
      left_join(label_positions, by = "Group") %>%
      filter(!is.na(LabelY))

    p <- ggplot(plot_df, aes(x = Group, y = FinalConcentration_ng_mL, fill = Group)) +
      geom_violin(trim = FALSE, alpha = 0.35, color = NA) +
      geom_boxplot(width = 0.18, outlier.shape = NA, alpha = 0.75, color = "black") +
      geom_jitter(aes(color = Group), width = 0.12, alpha = 0.45, size = 1.6) +
      geom_text(
        data = label_df,
        aes(x = Group, y = LabelY, label = IQR_Label),
        inherit.aes = FALSE,
        size = 3,
        angle = 90,
        hjust = 0,
        vjust = 0.5
      ) +
      labs(
        title = paste("Within-group dispersion of cortisol concentrations:", group_col),
        subtitle = paste(dataset_label, "— violin width shows distribution shape; labels show within-group IQR"),
        x = group_col,
        y = "Final cortisol concentration (ng/mL)"
      ) +
      theme_bw() +
      theme(
        legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1),
        plot.margin = margin(t = 10, r = 35, b = 10, l = 10)
      ) +
      coord_cartesian(clip = "off")

    if (exists("use_log10_y_axis") && isTRUE(use_log10_y_axis)) {
      p <- p + scale_y_log10() + labs(y = "Final cortisol concentration (ng/mL), log10 scale")
    }

    print(p)
  }
  dev.off()
  cat("Dispersion plots written to: ", plot_output_file, "\n", sep = "")

  invisible(list(summary = dispersion_summary, fligner = fligner_results))
}

# Full dataset dispersion analysis.
if (exists("plot_data")) {
  dispersion_full_results <- run_dispersion_analysis(
    data = plot_data,
    grouping_columns = grouping_columns,
    dataset_label = "Full dataset",
    summary_output_file = "cortisol_dispersion_summary_full_dataset.csv",
    fligner_output_file = "cortisol_dispersion_fligner_tests_full_dataset.csv",
    plot_output_file = "cortisol_dispersion_plots_full_dataset.pdf"
  )
} else {
  cat("Full dataset dispersion analysis skipped because plot_data was not found.\n")
}

Dispersion summary written to: cortisol_dispersion_summary_full_dataset.csv
Fligner-Killeen dispersion test results written to: cortisol_dispersion_fligner_tests_full_dataset.csv

Dispersion summary for Full dataset:
# A tibble: 38 × 13
   AnalysisDataset GroupingVariable Group     n median_ng_mL mean_ng_mL sd_ng_mL
   <chr>           <chr>            <chr> <int>        <dbl>      <dbl>    <dbl>
 1 Full dataset    Ploidy           Dipl…    19      0.0120     0.0131    0.0113
 2 Full dataset    Ploidy           Trip…    20      0.0124     0.0163    0.0183
 3 Full dataset    Dilution Factor  0         9      0.00788    0.0248    0.0285
 4 Full dataset    Dilution Factor  25        1      0.00504    0.00504  NA     
 5 Full dataset    Dilution Factor  26.5      1      0.00239    0.00239  NA     
 6 Full dataset    Dilution Factor  28        1      0.00633    0.00633  NA     
 7 Full dataset    Dilution Factor  29        1      0.0129     0.0129   NA     
 8 Full dataset    Dilution Factor  30        1      0.00143    0.00143  NA     
 9 Full dataset    Dilution Factor  34        1      0.0123     0.0123   NA     
10 Full dataset    Dilution Factor  35        1      0.0118     0.0118   NA     
# ℹ 28 more rows
# ℹ 6 more variables: variance_ng_mL <dbl>, IQR_ng_mL <dbl>, MAD_ng_mL <dbl>,
#   min_ng_mL <dbl>, max_ng_mL <dbl>, coefficient_of_variation_percent <dbl>

Fligner-Killeen tests for Full dataset:
# A tibble: 5 × 6
  AnalysisDataset GroupingVariable                 Test  Statistic P_Value Note 
  <chr>           <chr>                            <chr>     <dbl>   <dbl> <chr>
1 Full dataset    Ploidy                           Flig…     0.573  0.449  No s…
2 Full dataset    Dilution Factor                  Flig…     7.63   0.0220 Sign…
3 Full dataset    SampleTissue_MassorVolume_ForCo… Flig…     2.78   0.249  No s…
4 Full dataset    SampleTime                       Flig…     3.70   0.295  No s…
5 Full dataset    Treatment                        Flig…     2.63   0.269  No s…
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
for 'Full dataset — violin width shows distribution shape; labels show
within-group IQR' in 'mbcsToSbcs': - substituted for — (U+2014)
Warning: Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
for 'Full dataset — violin width shows distribution shape; labels show
within-group IQR' in 'mbcsToSbcs': - substituted for — (U+2014)
Warning: Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
for 'Full dataset — violin width shows distribution shape; labels show
within-group IQR' in 'mbcsToSbcs': - substituted for — (U+2014)
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
for 'Full dataset — violin width shows distribution shape; labels show
within-group IQR' in 'mbcsToSbcs': - substituted for — (U+2014)
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
for 'Full dataset — violin width shows distribution shape; labels show
within-group IQR' in 'mbcsToSbcs': - substituted for — (U+2014)
Dispersion plots written to: cortisol_dispersion_plots_full_dataset.pdf
# Outlier-removed dispersion analysis.
if (exists("plot_data_no_extreme_outliers")) {
  dispersion_no_outlier_results <- run_dispersion_analysis(
    data = plot_data_no_extreme_outliers,
    grouping_columns = grouping_columns,
    dataset_label = "Extreme outliers excluded",
    summary_output_file = "cortisol_dispersion_summary_extreme_outliers_removed.csv",
    fligner_output_file = "cortisol_dispersion_fligner_tests_extreme_outliers_removed.csv",
    plot_output_file = "cortisol_dispersion_plots_extreme_outliers_removed.pdf"
  )
} else {
  cat("Outlier-removed dispersion analysis skipped because plot_data_no_extreme_outliers was not found.\n")
}

Dispersion summary written to: cortisol_dispersion_summary_extreme_outliers_removed.csv
Fligner-Killeen dispersion test results written to: cortisol_dispersion_fligner_tests_extreme_outliers_removed.csv

Dispersion summary for Extreme outliers excluded:
# A tibble: 35 × 13
   AnalysisDataset GroupingVariable Group     n median_ng_mL mean_ng_mL sd_ng_mL
   <chr>           <chr>            <chr> <int>        <dbl>      <dbl>    <dbl>
 1 Extreme outlie… Ploidy           Dipl…    18      0.0119     0.0117   0.00970
 2 Extreme outlie… Ploidy           Trip…    18      0.0107     0.0106   0.00595
 3 Extreme outlie… Dilution Factor  0         6      0.00191    0.00834  0.0145 
 4 Extreme outlie… Dilution Factor  25        1      0.00504    0.00504 NA      
 5 Extreme outlie… Dilution Factor  26.5      1      0.00239    0.00239 NA      
 6 Extreme outlie… Dilution Factor  28        1      0.00633    0.00633 NA      
 7 Extreme outlie… Dilution Factor  29        1      0.0129     0.0129  NA      
 8 Extreme outlie… Dilution Factor  30        1      0.00143    0.00143 NA      
 9 Extreme outlie… Dilution Factor  34        1      0.0123     0.0123  NA      
10 Extreme outlie… Dilution Factor  35        1      0.0118     0.0118  NA      
# ℹ 25 more rows
# ℹ 6 more variables: variance_ng_mL <dbl>, IQR_ng_mL <dbl>, MAD_ng_mL <dbl>,
#   min_ng_mL <dbl>, max_ng_mL <dbl>, coefficient_of_variation_percent <dbl>

Fligner-Killeen tests for Extreme outliers excluded:
# A tibble: 5 × 6
  AnalysisDataset           GroupingVariable       Test  Statistic P_Value Note 
  <chr>                     <chr>                  <chr>     <dbl>   <dbl> <chr>
1 Extreme outliers excluded Ploidy                 Flig…     3.16   0.0754 No s…
2 Extreme outliers excluded Dilution Factor        Flig…     0.148  0.929  No s…
3 Extreme outliers excluded SampleTissue_MassorVo… Flig…     2.78   0.249  No s…
4 Extreme outliers excluded SampleTime             Flig…     6.75   0.0803 No s…
5 Extreme outliers excluded Treatment              Flig…     1.17   0.558  No s…
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
for 'Extreme outliers excluded — violin width shows distribution shape; labels
show within-group IQR' in 'mbcsToSbcs': - substituted for — (U+2014)
Warning: Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
for 'Extreme outliers excluded — violin width shows distribution shape; labels
show within-group IQR' in 'mbcsToSbcs': - substituted for — (U+2014)
Warning: Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
for 'Extreme outliers excluded — violin width shows distribution shape; labels
show within-group IQR' in 'mbcsToSbcs': - substituted for — (U+2014)
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
for 'Extreme outliers excluded — violin width shows distribution shape; labels
show within-group IQR' in 'mbcsToSbcs': - substituted for — (U+2014)
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
for 'Extreme outliers excluded — violin width shows distribution shape; labels
show within-group IQR' in 'mbcsToSbcs': - substituted for — (U+2014)
Dispersion plots written to: cortisol_dispersion_plots_extreme_outliers_removed.pdf

26 Optional: focused two-variable plot

This optional plot is also metadata-driven. It uses focused_plot_x_variable for the x-axis and focused_plot_color_variable for color. For your example metadata, this creates a SampleDay-by-Treatment plot. For a different metadata file, change those two variables in the user-input section or set either one to NULL to skip this plot.

if (!is.null(focused_plot_x_variable) && !is.null(focused_plot_color_variable)) {
  if (all(c(focused_plot_x_variable, focused_plot_color_variable) %in% colnames(plot_data))) {
    focused_two_variable_plot <- ggplot(
      plot_data,
      aes(
        x = as.factor(.data[[focused_plot_x_variable]]),
        y = FinalConcentration_ng_mL,
        color = as.factor(.data[[focused_plot_color_variable]])
      )
    ) +
      geom_boxplot(
        aes(group = interaction(.data[[focused_plot_x_variable]], .data[[focused_plot_color_variable]])),
        outlier.shape = NA,
        alpha = 0.35,
        position = position_dodge(width = 0.8)
      ) +
      geom_point(
        position = position_jitterdodge(jitter.width = 0.15, jitter.height = 0, dodge.width = 0.8),
        alpha = 0.8,
        size = 2
      ) +
      labs(
        title = paste("Final cortisol concentration by", focused_plot_x_variable, "and", focused_plot_color_variable),
        x = focused_plot_x_variable,
        y = "Final cortisol concentration (ng/mL)",
        color = focused_plot_color_variable
      ) +
      theme_bw()

    if (use_log10_y_axis) {
      focused_two_variable_plot <- focused_two_variable_plot +
        scale_y_log10() +
        labs(y = "Final cortisol concentration (ng/mL), log10 scale")
    }

    focused_two_variable_plot
  } else {
    message(
      "Focused two-variable plot skipped because one or both requested columns were not found: ",
      focused_plot_x_variable, ", ", focused_plot_color_variable
    )
  }
}

27 Clean ranked cortisol summary output

This section creates a simplified CSV for quick review. It includes only the sample ID, calculated %B/B0, and the final dilution-corrected cortisol concentration. Rows are sorted from highest to lowest final cortisol concentration.

# This output intentionally uses the same working/output location as the rest of
# the script. If a future version defines an output folder, this chunk will use it;
# otherwise it writes to the current working directory just like the other CSVs.
clean_ranked_cortisol_summary_output <- if (exists("output_directory")) {
  file.path(output_directory, "clean_ranked_cortisol_summary.csv")
} else {
  "clean_ranked_cortisol_summary.csv"
}

# Use the main processed metadata table created earlier in the workflow.
# This table contains the final concentration columns after blank correction,
# plate-specific standard curve back-calculation, dilution correction, and QC flags.
clean_summary_source <- processed_metadata

# If IsBlank / IsStandard columns are missing for any reason, create safe defaults.
if (!("IsBlank" %in% colnames(clean_summary_source))) {
  clean_summary_source$IsBlank <- grepl("^BLK$|blank", clean_summary_source$SampleID, ignore.case = TRUE)
}
if (!("IsStandard" %in% colnames(clean_summary_source))) {
  clean_summary_source$IsStandard <- grepl("^S[0-7]$", clean_summary_source$SampleID, ignore.case = TRUE)
}

# Choose the final concentration column used by this version of the script.
final_conc_col <- dplyr::case_when(
  "FinalConcentration_ng_mL" %in% colnames(clean_summary_source) ~ "FinalConcentration_ng_mL",
  "Final_Concentration_ng_mL" %in% colnames(clean_summary_source) ~ "Final_Concentration_ng_mL",
  TRUE ~ NA_character_
)

if (is.na(final_conc_col)) {
  stop("Could not find a final concentration column for the clean ranked cortisol summary.")
}

clean_ranked_cortisol_summary <- clean_summary_source %>%
  dplyr::mutate(
    IsBlank_clean = as.logical(IsBlank),
    IsStandard_clean = as.logical(IsStandard)
  ) %>%
  dplyr::filter(
    !is.na(SampleID),
    is.na(IsBlank_clean) | IsBlank_clean == FALSE,
    is.na(IsStandard_clean) | IsStandard_clean == FALSE
  ) %>%
  dplyr::transmute(
    SampleID = as.character(SampleID),
    Percent_B_B0 = Percent_B_B0,
    FinalConcentration_ng_mL = .data[[final_conc_col]]
  ) %>%
  dplyr::arrange(dplyr::desc(FinalConcentration_ng_mL), SampleID)

readr::write_csv(clean_ranked_cortisol_summary, clean_ranked_cortisol_summary_output)

cat("Clean ranked cortisol summary saved to: ", clean_ranked_cortisol_summary_output, "
", sep = "")
Clean ranked cortisol summary saved to: clean_ranked_cortisol_summary.csv
knitr::kable(
  clean_ranked_cortisol_summary,
  digits = 4,
  caption = "Clean ranked cortisol summary: samples ordered from highest to lowest final cortisol concentration."
)
Clean ranked cortisol summary: samples ordered from highest to lowest final cortisol concentration.
SampleID Percent_B_B0 FinalConcentration_ng_mL
T-48_0610_P4 81.3737 0.0717
T-48_0610_P2 90.7954 0.0624
D48-RECOVER24_0611_P1 92.1486 0.0391
D-48_0610_P4 82.4461 0.0373
D-CTL-19_0609_P2 122.0478 0.0208
T48-RECOVER24_0611_P1 122.0094 0.0206
D-CTL-43_0610_P3 123.8095 0.0204
D-CTL-19_0609_P4 124.8564 0.0192
T-CTL-48_0610_P1 125.0606 0.0189
D-CTL-19_0609_P3 125.7883 0.0181
T-CTL-19_0609_P1 126.5033 0.0173
D-CTL-48_0610_P1 126.5926 0.0172
T-CTL-48_0610_P3 128.1374 0.0156
T-19_0609_P4 128.3927 0.0154
D-CTL-48_0610_P3 128.9417 0.0149
T-CTL-19_0609_P2 129.0310 0.0148
T-CTL-48_0610_P2 129.3757 0.0145
D48-RECOVER24_0611_P4 122.4180 0.0129
T48-RECOVER24_0611_P3 131.8141 0.0125
T-19_0609_P1 125.7883 0.0123
D48-RECOVER24_0611_P2 130.6779 0.0120
D-19_0609_P3 126.9884 0.0118
D-CTL-48_0610_P4 136.2441 0.0095
T-19_0609_P2 136.7931 0.0091
T-CTL-48_0610_P4 137.1633 0.0089
T48-RECOVER24_0611_P4 103.3321 0.0079
D-CTL-48_0610_P2 140.0102 0.0075
T-19_0609_P3 133.3844 0.0063
T-48_0610_P3 135.2355 0.0050
T-CTL-19_0609_P3 115.4092 0.0040
T48-RECOVER24_0611_P2 115.5113 0.0038
T-CTL-19_0609_P4 116.0858 0.0030
D-19_0609_P2 122.7244 0.0024
D-48_0610_P2 115.0645 0.0024
D-CTL-19_0609_P1 117.2986 0.0016
D-19_0609_P4 116.5582 0.0014
D-48_0610_P1 129.7332 0.0014
T-48_0610_P1 112.2814 0.0010
D-43_0610_P3 108.1323 0.0000

28 Sanity check that %B/B0 and final concentration are correlated correctly

You should see:

👉 downward trend

If not: 👉 something is off in conversion/extrapolation

# =========================
# Sanity Check Plot
# =========================

sanity_plot_data <- processed_metadata

# Create grouping variable if needed
if (!"HatcheryLocation_Tissue" %in% colnames(sanity_plot_data)) {
  sanity_plot_data$HatcheryLocation_Tissue <- paste(
    sanity_plot_data$HatcheryLocation,
    sanity_plot_data$Tissue,
    sep = "_"
  )
}
Warning: Unknown or uninitialised column: `HatcheryLocation`.
# Assign colors by group
group_colors <- as.factor(sanity_plot_data$HatcheryLocation_Tissue)
color_palette <- rainbow(length(levels(group_colors)))
point_colors <- color_palette[group_colors]

# =========================
# Define extrapolated values
# =========================
lower_bound <- 0.005
upper_bound <- 50

if ("Concentration_ng_mL" %in% colnames(sanity_plot_data)) {
  is_extrapolated <- sanity_plot_data$Concentration_ng_mL < lower_bound |
                     sanity_plot_data$Concentration_ng_mL > upper_bound
} else {
  is_extrapolated <- rep(FALSE, nrow(sanity_plot_data))
}

is_extrapolated[is.na(is_extrapolated)] <- FALSE

print(table(is_extrapolated))
is_extrapolated
FALSE  TRUE 
    3    37 
# =========================
# Plot with shaded region
# =========================

# Set up empty plot first
plot(sanity_plot_data$Percent_B_B0,
     sanity_plot_data$FinalConcentration_ng_mL,
     type = "n",  # important: don't plot points yet
     xlab = "%B/B0",
     ylab = "Final Concentration",
     main = "Sanity Check: Should be inverse relationship")

# Add shaded ranges for interpolation and extrapolation
usr <- par("usr")
# Add full theoretical quantification area (0–100% B/B0)
rect(
  xleft = 0, xright = 100,
  ybottom = usr[3], ytop = usr[4],
  col = adjustcolor("lightgray", alpha.f = 0.15),
  border = NA
)

# Add dashed boundary lines for full 0–100% B/B0 scale
abline(v = c(0, 100), lty = 2, col = "gray40", lwd = 1.5)
# Add shaded optimal range (20–80% B/B0)
rect(
  xleft = 20, xright = 80,
  ybottom = usr[3], ytop = usr[4],
  col = adjustcolor("lightgreen", alpha.f = 0.25),
  border = NA
)

# Add dashed boundary lines for optimal range
abline(v = c(20, 80), lty = 2, col = "darkgreen", lwd = 1.5)


# Add points
points(sanity_plot_data$Percent_B_B0,
       sanity_plot_data$FinalConcentration_ng_mL,
       col = point_colors,
       pch = ifelse(is_extrapolated, 17, 16),
       cex = ifelse(is_extrapolated, 1.5, 1.1))

# Add regression line
abline(lm(FinalConcentration_ng_mL ~ Percent_B_B0,
          data = sanity_plot_data),
       col = "red", lwd = 2)

# Legend: groups
legend("topright",
       legend = levels(group_colors),
       col = color_palette,
       pch = 16,
       cex = 0.7,
       title = "Group")

# Legend: extrapolation
legend("bottomleft",
       legend = c("Interpolated", "Extrapolated"),
       pch = c(16, 17),
       col = c("black", "black"),
       cex = 0.8,
       title = "Point Type")

# Get plot limits
usr <- par("usr")

# Y position near top of plot (slightly below top edge)
y_pos <- usr[4] - 0.12 * (usr[4] - usr[3])

# Label optimal range (centered in 20–80 region)
text(x = 50,
     y = y_pos,
     labels = "Optimal Range (20–80%)",
     col = "darkgreen",
     cex = 0.9,
     font = 2)

# Label full range (centered in 0–100 region, slightly higher)
text(x = 50,
     y = usr[4] - 0.07 * (usr[4] - usr[3]),
     labels = "Full Assay Range (0–100%)",
     col = "gray40",
     cex = 0.9,
     font = 2)

#if everything is flagged the concentrations are mostly outside of the standard curve, which is still biologically relevant information
summary(sanity_plot_data$Concentration_ng_mL)
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max.      NA's 
3.147e-05 1.953e-04 3.078e-04 6.829e-04 4.232e-04 5.733e-03         1

29 Oyster size comparisons, size-cortisol correlations, and size-adjusted models

This section evaluates whether oyster shell size differs between hatchery locations and whether shell size is associated with cortisol concentration. Because each oyster can contribute multiple tissue samples, all size-based analyses are summarized at the oyster level to avoid treating multiple tissues from the same oyster as independent biological replicates.

# ============================================================
# Oyster size, size-cortisol correlations, and size-adjusted models
# ============================================================
# This section is designed to run on both:
#   1. the full processed dataset
#   2. the secondary dataset with extreme outliers removed, if available
#
# Key design choice:
#   Size comparisons are oyster-level only.
#   Size-cortisol correlations and size-adjusted models use one consistent tissue
#   per oyster, selected as the most frequently represented tissue with usable data.
# ============================================================

library(dplyr)
library(ggplot2)
library(readr)
library(broom)

# ---------- Helper functions ----------

first_non_missing <- function(x) {
  x <- x[!is.na(x) & x != ""]
  if (length(x) == 0) return(NA)
  x[1]
}

make_oyster_id <- function(data) {
  if ("Oyster" %in% colnames(data)) {
    as.character(data$Oyster)
  } else {
    # Fallback: use the portion of SampleID before the final underscore-delimited tissue token.
    # If your SampleID format changes, it is better to include an explicit Oyster column.
    gsub("(_[^_]+)$", "", as.character(data$SampleID))
  }
}

clean_sample_rows <- function(data) {
  data %>%
    mutate(
      SampleID_chr = as.character(SampleID),
      IsBlank_safe = grepl("^BLK$|blank", SampleID_chr, ignore.case = TRUE),
      IsStandard_safe = grepl("^S[0-7]$", SampleID_chr, ignore.case = TRUE)
    ) %>%
    filter(!IsBlank_safe, !IsStandard_safe)
}

has_required_size_columns <- function(data) {
  all(c("SampleID", "HatcheryLocation", "Shell_Length_(mm)", "Shell_Width_(mm)") %in% colnames(data))
}

prepare_oyster_size_data <- function(data) {
  if (!has_required_size_columns(data)) {
    return(tibble(
      Note = "Size analysis skipped: required columns SampleID, HatcheryLocation, Shell_Length_(mm), Shell_Width_(mm) were not all present."
    ))
  }

  data %>%
    clean_sample_rows() %>%
    mutate(
      OysterID = make_oyster_id(.),
      Shell_Length_mm = suppressWarnings(as.numeric(.data[["Shell_Length_(mm)"]])),
      Shell_Width_mm = suppressWarnings(as.numeric(.data[["Shell_Width_(mm)"]])),
      HatcheryLocation = as.factor(HatcheryLocation)
    ) %>%
    filter(!is.na(OysterID), OysterID != "", !is.na(HatcheryLocation)) %>%
    group_by(OysterID, HatcheryLocation) %>%
    summarize(
      Shell_Length_mm = first_non_missing(Shell_Length_mm),
      Shell_Width_mm = first_non_missing(Shell_Width_mm),
      n_rows_for_oyster = n(),
      .groups = "drop"
    ) %>%
    filter(!is.na(Shell_Length_mm) | !is.na(Shell_Width_mm))
}

run_two_group_wilcox <- function(data, response_col) {
  if (!response_col %in% colnames(data)) return(tibble())
  test_data <- data %>% filter(!is.na(.data[[response_col]]), !is.na(HatcheryLocation))
  if (n_distinct(test_data$HatcheryLocation) != 2 || nrow(test_data) < 3) {
    return(tibble(
      response = response_col,
      method = "Wilcoxon rank-sum test",
      note = "Test skipped: exactly two HatcheryLocation groups with sufficient data are required."
    ))
  }

  wt <- wilcox.test(as.formula(paste(response_col, "~ HatcheryLocation")), data = test_data, exact = FALSE)
  broom::tidy(wt) %>% mutate(response = response_col, .before = 1)
}

plot_size_by_location <- function(oyster_size_data, response_col, y_label, plot_title) {
  ggplot(oyster_size_data,
         aes(x = HatcheryLocation, y = .data[[response_col]], color = HatcheryLocation, fill = HatcheryLocation)) +
    geom_violin(trim = FALSE, alpha = 0.25, na.rm = TRUE) +
    geom_boxplot(width = 0.20, outlier.shape = NA, alpha = 0.45, na.rm = TRUE) +
    geom_jitter(width = 0.12, size = 2, alpha = 0.75, na.rm = TRUE) +
    theme_minimal() +
    theme(legend.position = "none") +
    labs(title = plot_title, x = "Hatchery location", y = y_label)
}

choose_representative_tissue <- function(data) {
  if (!"Tissue" %in% colnames(data)) return(NA_character_)

  tissue_counts <- data %>%
    clean_sample_rows() %>%
    mutate(
      OysterID = make_oyster_id(.),
      FinalConcentration_ng_mL = suppressWarnings(as.numeric(FinalConcentration_ng_mL))
    ) %>%
    filter(!is.na(Tissue), Tissue != "", !is.na(OysterID), !is.na(FinalConcentration_ng_mL)) %>%
    distinct(OysterID, Tissue) %>%
    count(Tissue, sort = TRUE)

  if (nrow(tissue_counts) == 0) return(NA_character_)
  as.character(tissue_counts$Tissue[1])
}

prepare_oyster_cortisol_size_data <- function(data, selected_tissue = NULL) {
  if (!has_required_size_columns(data) || !"FinalConcentration_ng_mL" %in% colnames(data)) {
    return(tibble(
      Note = "Correlation/model analysis skipped: required size columns or FinalConcentration_ng_mL were not present."
    ))
  }

  if (is.null(selected_tissue) || is.na(selected_tissue)) {
    selected_tissue <- choose_representative_tissue(data)
  }

  if (is.na(selected_tissue)) {
    return(tibble(
      Note = "Correlation/model analysis skipped: no usable Tissue column/data were available."
    ))
  }

  data %>%
    clean_sample_rows() %>%
    mutate(
      OysterID = make_oyster_id(.),
      Shell_Length_mm = suppressWarnings(as.numeric(.data[["Shell_Length_(mm)"]])),
      Shell_Width_mm = suppressWarnings(as.numeric(.data[["Shell_Width_(mm)"]])),
      Cortisol_ng_mL = suppressWarnings(as.numeric(FinalConcentration_ng_mL)),
      HatcheryLocation = as.factor(HatcheryLocation),
      SelectedTissue = selected_tissue
    ) %>%
    filter(Tissue == selected_tissue,
           !is.na(OysterID), OysterID != "",
           !is.na(HatcheryLocation),
           !is.na(Cortisol_ng_mL), Cortisol_ng_mL > 0) %>%
    group_by(OysterID, HatcheryLocation, SelectedTissue) %>%
    summarize(
      Shell_Length_mm = first_non_missing(Shell_Length_mm),
      Shell_Width_mm = first_non_missing(Shell_Width_mm),
      # If more than one row exists for the selected tissue/dilution, use the median
      # to avoid allowing one oyster to contribute multiple values.
      Cortisol_ng_mL = median(Cortisol_ng_mL, na.rm = TRUE),
      n_rows_for_oyster_selected_tissue = n(),
      .groups = "drop"
    ) %>%
    filter(!is.na(Cortisol_ng_mL))
}

run_spearman_tests <- function(data) {
  out <- list()

  if ("Shell_Length_mm" %in% colnames(data) && sum(!is.na(data$Shell_Length_mm) & !is.na(data$Cortisol_ng_mL)) >= 4) {
    out[["Shell_Length_mm"]] <- broom::tidy(cor.test(data$Shell_Length_mm, data$Cortisol_ng_mL,
                                                     method = "spearman", exact = FALSE)) %>%
      mutate(size_metric = "Shell_Length_mm", .before = 1)
  }

  if ("Shell_Width_mm" %in% colnames(data) && sum(!is.na(data$Shell_Width_mm) & !is.na(data$Cortisol_ng_mL)) >= 4) {
    out[["Shell_Width_mm"]] <- broom::tidy(cor.test(data$Shell_Width_mm, data$Cortisol_ng_mL,
                                                    method = "spearman", exact = FALSE)) %>%
      mutate(size_metric = "Shell_Width_mm", .before = 1)
  }

  if (length(out) == 0) {
    return(tibble(note = "Spearman correlations skipped: insufficient complete observations."))
  }

  bind_rows(out)
}

plot_size_cortisol <- function(data, size_col, x_label, plot_title) {
  ggplot(data, aes(x = .data[[size_col]], y = Cortisol_ng_mL, color = HatcheryLocation)) +
    geom_point(size = 2.3, alpha = 0.85, na.rm = TRUE) +
    geom_smooth(method = "lm", se = FALSE, na.rm = TRUE) +
    scale_y_log10() +
    theme_minimal() +
    labs(title = plot_title, x = x_label, y = "Final cortisol concentration (ng/mL, log10 scale)", color = "Hatchery location")
}

run_size_adjusted_models <- function(data) {
  if (!"Cortisol_ng_mL" %in% colnames(data)) return(tibble(note = "Model skipped: no cortisol column."))

  model_data <- data %>%
    mutate(log10_Cortisol = log10(Cortisol_ng_mL)) %>%
    filter(is.finite(log10_Cortisol), !is.na(HatcheryLocation))

  model_outputs <- list()

  for (size_col in c("Shell_Length_mm", "Shell_Width_mm")) {
    if (!size_col %in% colnames(model_data)) next
    md <- model_data %>% filter(!is.na(.data[[size_col]]))

    if (nrow(md) < 6 || n_distinct(md$HatcheryLocation) < 2) {
      model_outputs[[size_col]] <- tibble(
        model = paste("log10_Cortisol ~ HatcheryLocation +", size_col),
        note = "Model skipped: insufficient complete observations or fewer than two hatchery locations."
      )
      next
    }

    # Main size-adjusted model: hatchery location + size covariate
    main_formula <- as.formula(paste("log10_Cortisol ~ HatcheryLocation +", size_col))
    main_model <- lm(main_formula, data = md)
    main_terms <- broom::tidy(main_model) %>%
      mutate(model = paste("log10_Cortisol ~ HatcheryLocation +", size_col), .before = 1)

    main_glance <- broom::glance(main_model) %>%
      mutate(model = paste("log10_Cortisol ~ HatcheryLocation +", size_col), .before = 1)

    # Interaction model: tests whether the size-cortisol slope differs by hatchery location
    interaction_formula <- as.formula(paste("log10_Cortisol ~ HatcheryLocation *", size_col))
    interaction_terms <- tibble()
    interaction_glance <- tibble()
    interaction_anova <- tibble()

    if (nrow(md) >= 8) {
      interaction_model <- lm(interaction_formula, data = md)
      interaction_terms <- broom::tidy(interaction_model) %>%
        mutate(model = paste("log10_Cortisol ~ HatcheryLocation *", size_col), .before = 1)
      interaction_glance <- broom::glance(interaction_model) %>%
        mutate(model = paste("log10_Cortisol ~ HatcheryLocation *", size_col), .before = 1)
      interaction_anova <- broom::tidy(anova(main_model, interaction_model)) %>%
        mutate(model_comparison = paste("Interaction test for", size_col), .before = 1)
    }

    model_outputs[[size_col]] <- list(
      terms = bind_rows(main_terms, interaction_terms),
      fit = bind_rows(main_glance, interaction_glance),
      interaction_test = interaction_anova
    )
  }

  model_outputs
}

save_model_outputs <- function(model_outputs, prefix) {
  if (is.data.frame(model_outputs)) {
    write_csv(model_outputs, paste0(prefix, "_size_adjusted_model_notes.csv"))
    return(invisible(NULL))
  }

  terms <- purrr::map_dfr(model_outputs, function(x) {
    if (is.list(x) && "terms" %in% names(x)) {
      return(x$terms)
    } else {
      return(tibble())
    }
  })

  fit <- purrr::map_dfr(model_outputs, function(x) {
    if (is.list(x) && "fit" %in% names(x)) {
      return(x$fit)
    } else {
      return(tibble())
    }
  })

  interaction <- purrr::map_dfr(model_outputs, function(x) {
    if (is.list(x) && "interaction_test" %in% names(x)) {
      return(x$interaction_test)
    } else {
      return(tibble())
    }
  })

  readr::write_csv(terms, paste0(prefix, "_size_adjusted_model_terms.csv"))
  readr::write_csv(fit, paste0(prefix, "_size_adjusted_model_fit.csv"))
  readr::write_csv(interaction, paste0(prefix, "_size_adjusted_interaction_tests.csv"))

  return(list(
    terms = terms,
    fit = fit,
    interaction = interaction
  ))
}
run_oyster_size_workflow <- function(data, dataset_label, output_prefix) {
  cat("\n==============================\n")
  cat("Oyster size workflow:", dataset_label, "\n")
  cat("==============================\n")

  if (!has_required_size_columns(data)) {
    cat("Skipping oyster size workflow because required shell size columns are missing.\n")
    return(invisible(NULL))
  }

  # 1. Oyster-level size comparisons between hatchery locations
  oyster_size_data <- prepare_oyster_size_data(data)
  write_csv(oyster_size_data, paste0(output_prefix, "_oyster_size_by_individual.csv"))

  size_summary <- oyster_size_data %>%
    group_by(HatcheryLocation) %>%
    summarize(
      n_oysters = n(),
      mean_length_mm = mean(Shell_Length_mm, na.rm = TRUE),
      median_length_mm = median(Shell_Length_mm, na.rm = TRUE),
      sd_length_mm = sd(Shell_Length_mm, na.rm = TRUE),
      mean_width_mm = mean(Shell_Width_mm, na.rm = TRUE),
      median_width_mm = median(Shell_Width_mm, na.rm = TRUE),
      sd_width_mm = sd(Shell_Width_mm, na.rm = TRUE),
      .groups = "drop"
    )
  write_csv(size_summary, paste0(output_prefix, "_oyster_size_summary_by_hatchery_location.csv"))
  print(size_summary)

  size_tests <- bind_rows(
    run_two_group_wilcox(oyster_size_data, "Shell_Length_mm"),
    run_two_group_wilcox(oyster_size_data, "Shell_Width_mm")
  )
  write_csv(size_tests, paste0(output_prefix, "_oyster_size_wilcoxon_tests.csv"))
  print(size_tests)

  p_size_length <- plot_size_by_location(
    oyster_size_data, "Shell_Length_mm", "Shell length (mm)",
    paste("Shell length by hatchery location -", dataset_label)
  )
  p_size_width <- plot_size_by_location(
    oyster_size_data, "Shell_Width_mm", "Shell width (mm)",
    paste("Shell width by hatchery location -", dataset_label)
  )

  pdf(paste0(output_prefix, "_oyster_size_by_hatchery_location.pdf"), width = 7, height = 5)
  print(p_size_length)
  print(p_size_width)
  dev.off()

  # 2. Choose the most consistently represented tissue for oyster-level cortisol-size correlations
  selected_tissue <- choose_representative_tissue(data)
  cat("Selected tissue for size-cortisol correlation/modeling:", selected_tissue, "\n")
  write_csv(tibble(Dataset = dataset_label, SelectedTissue = selected_tissue),
            paste0(output_prefix, "_selected_tissue_for_size_cortisol_analysis.csv"))

  corr_data <- prepare_oyster_cortisol_size_data(data, selected_tissue)
  write_csv(corr_data, paste0(output_prefix, "_oyster_size_cortisol_correlation_data.csv"))

  if ("Cortisol_ng_mL" %in% colnames(corr_data) && nrow(corr_data) >= 4) {
    corr_tests <- run_spearman_tests(corr_data)
    write_csv(corr_tests, paste0(output_prefix, "_spearman_size_cortisol_correlations.csv"))
    print(corr_tests)

    p_corr_length <- plot_size_cortisol(
      corr_data, "Shell_Length_mm", "Shell length (mm)",
      paste("Cortisol vs shell length -", selected_tissue, "-", dataset_label)
    )
    p_corr_width <- plot_size_cortisol(
      corr_data, "Shell_Width_mm", "Shell width (mm)",
      paste("Cortisol vs shell width -", selected_tissue, "-", dataset_label)
    )

    pdf(paste0(output_prefix, "_size_cortisol_correlation_plots.pdf"), width = 7, height = 5)
    print(p_corr_length)
    print(p_corr_width)
    dev.off()

    # 3. Size-adjusted models
    model_outputs <- run_size_adjusted_models(corr_data)
    saved_models <- save_model_outputs(model_outputs, output_prefix)
    if (!is.null(saved_models)) {
      cat("\nSize-adjusted model terms:\n")
      print(saved_models$terms)
      cat("\nInteraction tests comparing additive vs interaction models:\n")
      print(saved_models$interaction)
    }
  } else {
    cat("Size-cortisol correlations and size-adjusted models skipped because fewer than 4 complete oyster-level observations were available.\n")
  }

  invisible(TRUE)
}

# ---------- Run workflow for full data ----------
if (exists("plot_data") && nrow(plot_data) > 0) {
  run_oyster_size_workflow(
    data = plot_data,
    dataset_label = "Full dataset",
    output_prefix = "full_dataset"
  )
} else if (exists("processed_metadata") && nrow(processed_metadata) > 0) {
  run_oyster_size_workflow(
    data = processed_metadata,
    dataset_label = "Full dataset",
    output_prefix = "full_dataset"
  )
}

==============================
Oyster size workflow: Full dataset 
==============================
Skipping oyster size workflow because required shell size columns are missing.
# ---------- Run workflow for extreme-outlier-removed data ----------
if (exists("plot_data_no_extreme_outliers") && nrow(plot_data_no_extreme_outliers) > 0) {
  run_oyster_size_workflow(
    data = plot_data_no_extreme_outliers,
    dataset_label = "Extreme outliers removed",
    output_prefix = "outliers_removed_dataset"
  )
} else {
  cat("Extreme-outlier-removed oyster size workflow skipped because plot_data_no_extreme_outliers was not available.\n")
}

==============================
Oyster size workflow: Extreme outliers removed 
==============================
Skipping oyster size workflow because required shell size columns are missing.

30 Notes for adapting this workflow

getwd()

getwd()

list.files()

plate_file_map

file.exists(“Plate1_Tecan_RawData.xlsx”) read_tecan_plate(“Plate1_Tecan_RawData.xlsx”, plate_number = 1)

plate_data <- plate_file_map %>% pmap(read_tecan_plate) str(plate_file_map)

args(read_tecan_plate) names(plate_file_map)

plate_file_map <- plate_file_map %>% rename( tecan_excel_file = TecanExcelFile, plate_number = Plate Number )

names(plate_file_map)

plate_data <- plate_file_map %>% pmap(read_tecan_plate)

plate_data str(plate_data)

setwd(“/Users/sydneyirwin/Desktop/Cortisol Coding”)

getwd() list.files()

rstudioapi::getSourceEditorContext()$path

getwd()

metadata[41, ]

metadata <- read_csv(metadata_csv_file, show_col_types = FALSE)