Resazurin Assays - USDA M.gigas Families in Response to Temperature Stress

2026
Resazurin
Pacific oyster
Crassostrea gigas
Magellana gigas
Author

Sam White

Published

April 15, 2026

INTRO

Adult oysters from 9 families were submerged in 85 mL of room temperature resazurin working solution in clear plastic cups and held at 36°C. At each designated timepoint, 250 µL was taken from each cup and pipetted into a clear 96-well plate. Fluorescence was measured using a Synergy HTX (Agilent) plate reader. At an undetermined point during the experiment, the temperature was increased to 40°C to accelerate metabolism.

Important

There were some cups with duplicate numbers, so they have been excluded from the analysis. See the “Excluded samples” section below for details.

The post below was knitted from 01.00-resazurin-20260415-36C.Rmd (GitHub)

RESULTS

There were no significant differences in metabolism AUC (metabolism) between families (ANOVA p = 0.15).

1 Background

Adult oysters from 9 families were submerged in 85 mL of room temperature resazurin working solution in clear plastic cups and held at 36°C. At each designated timepoint, 250 µL was taken from each cup and pipetted into a clear 96-well plate. Fluorescence was measured using a Synergy HTX (Agilent) plate reader. At an undetermined point during the experiment, the temperature was increased to 40°C to accelerate metabolism.

See Resazurin/data/20260415-36C/README.md for full experimental notes.

1.1 Expected inputs

Path Description
Resazurin/data/20260415-36C/plate-*-T*.txt Plate reader fluorescence exports (one file per plate per timepoint)
Resazurin/data/20260415-36C/layout.csv Well metadata: plate ID, well ID, blank flag, family groups, sample IDs, area measurements (mm², from ImageJ)

1.2 Expected outputs

All outputs are written to Resazurin/outputs/01.00-resazurin-20260415-36C/.

File Description
figures/ All plots generated by this script
auc_all_metrics.csv Per-individual AUC values for every active measurement metric
auc_summary.csv Group-level AUC summary statistics (mean, SD, SE, median)
metabolism.csv Full per-well per-timepoint metabolism data frame
pairwise_stats.csv Tukey-adjusted pairwise comparisons from AUC linear models

2 Setup

2.1 Knitr options

knitr::opts_chunk$set(
  echo = TRUE,         # Display code chunks
  eval = TRUE,        # Evaluate code chunks
  warning = FALSE,     # Hide warnings
  message = FALSE,     # Hide messages
  comment = "",         # Prevents appending '##' to beginning of lines in code output
  results = 'hold'     # Holds output so it's all printed together after code chunk
)

2.2 Load libraries

library(tidyverse)
library(pracma)       # trapz()
library(lme4)
library(lmerTest)
library(emmeans)
library(multcompView)
library(cowplot)
library(colorspace)   # qualitative_hcl() for large palettes

3 Helper Functions

normalize_well_id <- function(x) {
  x <- toupper(trimws(x))
  valid <- str_detect(x, "^[A-Z]+[0-9]+$")
  out <- rep(NA_character_, length(x))
  if (!any(valid)) return(out)
  m <- str_match(x[valid], "^([A-Z]+)([0-9]+)$")
  out[valid] <- paste0(m[, 2], as.integer(m[, 3]))
  out
}

parse_time_hr <- function(path) {
  hit <- str_match(basename(path),
                   "(?i)-T([0-9]+(?:\\.[0-9]+)?)\\.txt$")
  as.numeric(hit[, 2])
}

parse_plate_id <- function(path) {
  hit <- str_match(basename(path),
    "(?i)^plate-([A-Za-z0-9-]+)-T[0-9]+(?:\\.[0-9]+)?\\.txt$")
  id <- hit[, 2]
  ifelse(is.na(id), "unknown", id)
}

extract_results_block <- function(lines) {
  results_idx <- which(trimws(lines) == "Results")
  if (length(results_idx) == 0) stop("No Results section found")
  idx <- results_idx[1]
  header_tokens <- str_split(lines[idx + 1], "\\t")[[1]] |> trimws()
  col_ids <- header_tokens[
    header_tokens != "" & str_detect(header_tokens, "^[0-9]+$")]
  j <- idx + 2
  data_lines <- character()
  while (j <= length(lines)) {
    line <- lines[j]
    if (trimws(line) == "") break
    if (!str_detect(line, "^[A-Za-z]\\t")) break
    data_lines <- c(data_lines, line)
    j <- j + 1
  }
  list(col_ids = col_ids, data_lines = data_lines)
}

parse_plate_export <- function(path) {
  lines <- readLines(path, warn = FALSE)
  res <- extract_results_block(lines)

  map_dfr(res$data_lines, function(line) {
    tokens <- str_split(line, "\\t")[[1]] |> trimws()
    tokens <- tokens[tokens != ""]
    row_letter <- tokens[1]
    nums <- suppressWarnings(as.numeric(tokens[-1]))
    valid_idx <- which(!is.na(nums))
    if (length(valid_idx) == 0) return(tibble())
    vals <- nums[valid_idx]
    n <- min(length(vals), length(res$col_ids))
    tibble(
      row_id  = toupper(row_letter),
      col_id  = as.integer(res$col_ids[seq_len(n)]),
      well_id = normalize_well_id(
        paste0(toupper(row_letter), res$col_ids[seq_len(n)])),
      value   = vals[seq_len(n)]
    )
  }) %>%
    mutate(
      plate_id = str_to_lower(parse_plate_id(path)),
      time_hr  = parse_time_hr(path)
    )
}

trapezoid_auc <- function(time_hr, value) {
  ok <- is.finite(time_hr) & is.finite(value)
  t <- time_hr[ok]
  v <- value[ok]
  if (length(t) < 2) return(NA_real_)
  ord <- order(t)
  t <- t[ord]; v <- v[ord]
  sum(diff(t) * (head(v, -1) + tail(v, -1)) / 2)
}

# Shared helper: extract display unit string from a measurement column name.
# e.g. "area_mm2_measurement" -> "mm²", "weight_mg_measurement" -> "mg"
parse_meas_unit <- function(col_name) {
  unit_raw <- col_name |>
    str_remove("^metabolism_per_") |>
    str_remove("_measurement$") |>
    str_extract("[^_]+$")
  case_when(
    unit_raw == "mm2" ~ "mm²",
    unit_raw == "cm2" ~ "cm²",
    unit_raw == "mm3" ~ "mm³",
    unit_raw == "cm3" ~ "cm³",
    TRUE              ~ unit_raw
  )
}

# y-axis label for metabolism line plots: "fold change/mm²"
metabolism_y_label <- function(col_name) {
  paste0("Metabolism (fold change/", parse_meas_unit(col_name), ")")
}

# y-axis label for AUC box plots: "Metabolism (AUC; mm²)"
auc_y_label <- function(metric_name) {
  paste0("Metabolism (AUC; ", parse_meas_unit(metric_name), ")")
}

4 Load Data

4.1 Plate export files

proj_root <- rprojroot::find_rstudio_root_file()
data_dir  <- file.path(proj_root, "Resazurin", "data", "20260415-36C")
out_dir   <- file.path(proj_root, "Resazurin", "outputs",
                        "01.00-resazurin-20260415-36C")
fig_dir   <- file.path(out_dir, "figures")

dir.create(fig_dir, recursive = TRUE, showWarnings = FALSE)
dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)

plate_files <- list.files(
  data_dir,
  pattern = "(?i)^plate-.*-T[0-9]+(?:\\.[0-9]+)?\\.txt$",
  full.names = TRUE
)

plate_raw <- map_dfr(plate_files, function(path) {
  tryCatch(parse_plate_export(path),
           error = function(e) {
             message("Parse error in ", basename(path), ": ", e$message)
             tibble()
           })
})

str(plate_raw)
tibble [960 × 6] (S3: tbl_df/tbl/data.frame)
 $ row_id  : chr [1:960] "A" "A" "A" "A" ...
 $ col_id  : int [1:960] 1 2 3 4 5 6 7 8 9 10 ...
 $ well_id : chr [1:960] "A1" "A2" "A3" "A4" ...
 $ value   : num [1:960] 727 599 672 685 639 720 683 686 794 975 ...
 $ plate_id: chr [1:960] "a" "a" "a" "a" ...
 $ time_hr : num [1:960] 1 1 1 1 1 1 1 1 1 1 ...

4.2 Plate consistency check

Checks that every plate has the same number of wells at every timepoint. The expected well count is the mode across all plate × timepoint reads. Any plate with at least one deviating read is flagged and dropped entirely before any further analysis — removing only the aberrant timepoint would break the fold-change baseline calculation.

well_counts <- plate_raw %>%
  group_by(plate_id, time_hr) %>%
  summarise(n_wells = n_distinct(well_id), .groups = "drop")

expected_n_wells <- as.integer(
  names(which.max(table(well_counts$n_wells)))
)

inconsistent_reads <- well_counts %>%
  filter(n_wells != expected_n_wells) %>%
  arrange(plate_id, time_hr)

inconsistent_plate_ids <- unique(inconsistent_reads$plate_id)

if (nrow(inconsistent_reads) > 0) {
  cat("**Plate consistency check FAILED.**",
      "Expected", expected_n_wells, "wells per plate-timepoint read.",
      length(inconsistent_plate_ids),
      "plate(s) have at least one deviating read and are excluded",
      "from all analyses:\n\n")
  cat(knitr::kable(
    inconsistent_reads,
    col.names = c("Plate", "Time (h)", "Wells read"),
    caption   = paste("Expected:", expected_n_wells, "wells per read")
  ), sep = "\n")
  cat("\n")
  plate_raw <- plate_raw %>%
    filter(!plate_id %in% inconsistent_plate_ids)
  message(length(inconsistent_plate_ids),
          " plate(s) removed from plate_raw: ",
          paste(inconsistent_plate_ids, collapse = ", "))
} else {
  cat("Plate consistency check passed: all",
      n_distinct(well_counts$plate_id), "plates have",
      expected_n_wells, "wells at every timepoint.\n")
}

Plate consistency check passed: all 10 plates have 96 wells at every timepoint.

4.3 Layout file

layout_path <- file.path(data_dir, "layout.csv")

layout_raw <- read_csv(layout_path,
                       col_types = cols(.default = "c"),
                       show_col_types = FALSE)

# Standardise column names to snake_case
names(layout_raw) <- names(layout_raw) |>
  str_to_lower() |>
  str_replace_all("[^a-z0-9]+", "_") |>
  str_replace_all("_+", "_") |>
  str_replace("_$", "")

# Normalise plate_id to match plate file ids (strip "plate-" prefix)
layout_clean <- layout_raw %>%
  mutate(
    plate_id = str_remove(str_to_lower(plate_id), "^plate-"),
    well_id  = normalize_well_id(plate_well),
    is_blank = if ("is_blank" %in% names(layout_raw))
      toupper(trimws(is_blank)) %in% c("TRUE", "T", "1", "YES", "Y")
    else
      FALSE
  )

found_exclude_col <- intersect(
  c("exclude_from_analysis", "exclude", "omit", "not_analyzed"),
  names(layout_clean)
)[1]
layout_clean <- layout_clean %>%
  mutate(
    exclude_from_analysis = if (!is.na(found_exclude_col))
      toupper(trimws(.data[[found_exclude_col]])) %in%
        c("TRUE", "T", "1", "YES", "Y")
    else
      FALSE
  )

# Identify measurement columns and group columns
measurement_cols <- names(layout_clean)[
  str_detect(names(layout_clean), "_measurement$")]
group_cols <- names(layout_clean)[
  str_detect(names(layout_clean), "_group$")]

# Cast measurement columns to numeric
layout_clean <- layout_clean %>%
  mutate(across(all_of(measurement_cols),
                ~ suppressWarnings(as.numeric(.x))))

# Determine which measurement columns actually contain finite data
active_meas_cols <- measurement_cols[
  sapply(measurement_cols, function(col)
    any(is.finite(layout_clean[[col]]), na.rm = TRUE))]

# Normalise group values to lowercase so they match colour scale definitions
layout_clean <- layout_clean %>%
  mutate(across(all_of(group_cols),
                ~ str_to_lower(trimws(as.character(.x)))))

message("Group columns: ", paste(group_cols, collapse = ", "))
message("Active measurement columns: ",
        paste(active_meas_cols, collapse = ", "))

str(layout_clean)
tibble [960 × 14] (S3: tbl_df/tbl/data.frame)
 $ plate_id             : chr [1:960] "a" "a" "a" "a" ...
 $ plate_well           : chr [1:960] "A01" "A02" "A03" "A04" ...
 $ is_blank             : logi [1:960] FALSE FALSE FALSE FALSE FALSE FALSE ...
 $ family_id_group      : chr [1:960] "1" "9" "9" "9" ...
 $ sample_id_group      : chr [1:960] "47" "38" "52" "57" ...
 $ exclude_from_analysis: logi [1:960] FALSE FALSE FALSE FALSE FALSE FALSE ...
 $ exclude_reason       : chr [1:960] NA NA NA NA ...
 $ weight_g_measurement : num [1:960] NA NA NA NA NA NA NA NA NA NA ...
 $ width_mm_measurement : num [1:960] NA NA NA NA NA NA NA NA NA NA ...
 $ length_mm_measurement: num [1:960] NA NA NA NA NA NA NA NA NA NA ...
 $ treatment_group      : chr [1:960] NA NA NA NA ...
 $ area_mm2_measurement : num [1:960] 1346 648 577 682 372 ...
 $ imagej_id            : chr [1:960] "3" "11" "7" "8" ...
 $ well_id              : chr [1:960] "A1" "A2" "A3" "A4" ...

5 Merge Plate Data with Layout

dat <- plate_raw %>%
  left_join(
    layout_clean %>%
      select(plate_id, well_id, is_blank, exclude_from_analysis,
             any_of("exclude_reason"),
             all_of(group_cols), all_of(measurement_cols)),
    by = c("plate_id", "well_id")
  ) %>%
  mutate(
    is_blank = replace_na(is_blank, FALSE),
    exclude_from_analysis = replace_na(exclude_from_analysis, FALSE)
  )

str(dat)
tibble [960 × 16] (S3: tbl_df/tbl/data.frame)
 $ row_id               : chr [1:960] "A" "A" "A" "A" ...
 $ col_id               : int [1:960] 1 2 3 4 5 6 7 8 9 10 ...
 $ well_id              : chr [1:960] "A1" "A2" "A3" "A4" ...
 $ value                : num [1:960] 727 599 672 685 639 720 683 686 794 975 ...
 $ plate_id             : chr [1:960] "a" "a" "a" "a" ...
 $ time_hr              : num [1:960] 1 1 1 1 1 1 1 1 1 1 ...
 $ is_blank             : logi [1:960] FALSE FALSE FALSE FALSE FALSE FALSE ...
 $ exclude_from_analysis: logi [1:960] FALSE FALSE FALSE FALSE FALSE FALSE ...
 $ exclude_reason       : chr [1:960] NA NA NA NA ...
 $ family_id_group      : chr [1:960] "1" "9" "9" "9" ...
 $ sample_id_group      : chr [1:960] "47" "38" "52" "57" ...
 $ treatment_group      : chr [1:960] NA NA NA NA ...
 $ weight_g_measurement : num [1:960] NA NA NA NA NA NA NA NA NA NA ...
 $ width_mm_measurement : num [1:960] NA NA NA NA NA NA NA NA NA NA ...
 $ length_mm_measurement: num [1:960] NA NA NA NA NA NA NA NA NA NA ...
 $ area_mm2_measurement : num [1:960] 1346 648 577 682 372 ...

6 Raw Fluorescence

6.1 Data frame

# Wells in the plate reader output that have no layout entry get all-NA group
# columns after the join. Keep only wells assigned to at least one group.
active_gc <- intersect(group_cols, names(dat))

raw_df <- dat %>%
  filter(
    !is_blank,
    if (length(active_gc) > 0)
      if_any(all_of(active_gc), ~ !is.na(.))
    else
      TRUE
  ) %>%
  mutate(
    trace_id = if_else(
      !is.na(sample_id_group) & trimws(as.character(sample_id_group)) != "",
      as.character(sample_id_group),
      paste(plate_id, well_id, sep = "_")
    )
  )

families   <- sort(unique(na.omit(raw_df$family_id_group)))
treatments <- sort(unique(na.omit(raw_df$treatment_group)))

n_fam <- length(families)
n_trt <- length(treatments)

# Palette strategy:
#   <= 7 groups : Okabe-Ito (gold standard for colorblind-safe figures).
#   >  7 groups : colorspace::qualitative_hcl("Dynamic") scales to any N
#                 using perceptually uniform HCL space — no colour collisions.
# Black (#000000) is excluded from both and reserved for blank wells.
okabe_ito_7 <- c(
  "#E69F00", "#56B4E9", "#009E73", "#F0E442",
  "#0072B2", "#D55E00", "#CC79A7"
)
make_palette <- function(n) {
  if (n == 0L) return(character(0))
  if (n <= length(okabe_ito_7)) return(okabe_ito_7[seq_len(n)])
  colorspace::qualitative_hcl(n, palette = "Dynamic")
}

all_colours   <- make_palette(n_fam + n_trt)
fam_colours   <- setNames(all_colours[seq_len(n_fam)], families)
trt_colours   <- setNames(all_colours[n_fam + seq_len(n_trt)], treatments)

lty_pool <- c("solid", "dashed", "dotted", "dotdash", "longdash")
trt_linetypes <- setNames(
  lty_pool[(seq_len(n_trt) - 1L) %% length(lty_pool) + 1L],
  treatments
)
plate_well_colours <- c(blank = "black", fam_colours)

has_trt <- n_trt > 0

str(raw_df)
tibble [700 × 17] (S3: tbl_df/tbl/data.frame)
 $ row_id               : chr [1:700] "A" "A" "A" "A" ...
 $ col_id               : int [1:700] 1 2 3 4 5 6 7 8 9 10 ...
 $ well_id              : chr [1:700] "A1" "A2" "A3" "A4" ...
 $ value                : num [1:700] 727 599 672 685 639 720 683 686 794 975 ...
 $ plate_id             : chr [1:700] "a" "a" "a" "a" ...
 $ time_hr              : num [1:700] 1 1 1 1 1 1 1 1 1 1 ...
 $ is_blank             : logi [1:700] FALSE FALSE FALSE FALSE FALSE FALSE ...
 $ exclude_from_analysis: logi [1:700] FALSE FALSE FALSE FALSE FALSE FALSE ...
 $ exclude_reason       : chr [1:700] NA NA NA NA ...
 $ family_id_group      : chr [1:700] "1" "9" "9" "9" ...
 $ sample_id_group      : chr [1:700] "47" "38" "52" "57" ...
 $ treatment_group      : chr [1:700] NA NA NA NA ...
 $ weight_g_measurement : num [1:700] NA NA NA NA NA NA NA NA NA NA ...
 $ width_mm_measurement : num [1:700] NA NA NA NA NA NA NA NA NA NA ...
 $ length_mm_measurement: num [1:700] NA NA NA NA NA NA NA NA NA NA ...
 $ area_mm2_measurement : num [1:700] 1346 648 577 682 372 ...
 $ trace_id             : chr [1:700] "47" "38" "52" "57" ...

6.2 Raw fluorescence by plate (including blanks)

p_raw_plates <- dat %>%
  filter(is.finite(time_hr), is.finite(value)) %>%
  mutate(
    colour_group = if_else(is_blank, "blank",
                           coalesce(family_id_group, "sample")),
    trace_id     = paste(plate_id, well_id, sep = "_")
  ) %>%
  ggplot(aes(x = time_hr, y = value,
             group = trace_id, colour = colour_group)) +
  geom_line(alpha = 0.6) +
  geom_point(size = 1, alpha = 0.7) +
  facet_wrap(~ plate_id) +
  scale_colour_manual(
    values   = plate_well_colours,
    name     = "Group",
    breaks   = names(plate_well_colours),
    na.value = "grey80"
  ) +
  labs(x = "Time (h)", y = "Raw fluorescence (RFU)") +
  theme_classic(base_size = 12) +
  theme(strip.background = element_blank(),
        strip.text       = element_text(face = "bold"))

p_raw_plates

ggsave(file.path(fig_dir, "raw_fluor_by_plate.png"),
       p_raw_plates, width = 10, height = 8)

6.3 Mean raw fluorescence by family

raw_family_summary <- raw_df %>%
  filter(!is.na(family_id_group), !exclude_from_analysis) %>%
  group_by(family_id_group, treatment_group, time_hr) %>%
  summarise(
    mean_fluor = mean(value, na.rm = TRUE),
    se_fluor   = sd(value, na.rm = TRUE) /
      sqrt(sum(!is.na(value))),
    n          = sum(!is.na(value)),
    .groups    = "drop"
  ) %>%
  mutate(group_var = if (has_trt)
    paste(family_id_group, treatment_group, sep = ".")
  else
    family_id_group)

p_raw_mean <- ggplot(raw_family_summary,
    aes(x = time_hr, y = mean_fluor,
        colour = family_id_group,
        group = group_var)) +
  geom_ribbon(aes(ymin = mean_fluor - se_fluor,
                  ymax = mean_fluor + se_fluor,
                  fill = family_id_group),
              alpha = 0.15, colour = NA) +
  geom_line(
    mapping   = if (has_trt) aes(linetype = treatment_group) else NULL,
    linewidth = 1) +
  geom_point(size = 2) +
  scale_colour_manual(values = fam_colours, name = "Family") +
  scale_fill_manual(values = fam_colours, name = "Family") +
  labs(x = "Time (h)", y = "Mean raw fluorescence (RFU ± SE)") +
  theme_classic(base_size = 13) +
  if (has_trt) scale_linetype_manual(values = trt_linetypes, name = "Treatment") else NULL

p_raw_mean

ggsave(file.path(fig_dir, "raw_mean_by_family.png"),
       p_raw_mean, width = 8, height = 5)

6.4 Individual raw fluorescence traces by family

p_raw_by_family <- raw_df %>%
  filter(!is.na(family_id_group)) %>%
  ggplot(aes(x = time_hr, y = value, group = trace_id,
             colour = .data[[if (has_trt) "treatment_group" else "family_id_group"]])) +
  geom_line(alpha = 0.6) +
  geom_point(size = 1.2, alpha = 0.7) +
  facet_wrap(~ family_id_group) +
  scale_colour_manual(
    values = if (has_trt) trt_colours else fam_colours,
    name   = if (has_trt) "Treatment" else "Family") +
  labs(x = "Time (h)", y = "Raw fluorescence (RFU)") +
  theme_classic(base_size = 12) +
  theme(strip.background = element_blank(),
        strip.text       = element_text(face = "bold"))

p_raw_by_family

ggsave(file.path(fig_dir, "raw_individual_by_family.png"),
       p_raw_by_family, width = 10, height = 5)

6.5 Individual raw fluorescence traces by treatment

if (has_trt) {
  p_raw_by_treatment <- raw_df %>%
    ggplot(aes(x = time_hr, y = value,
               group = trace_id, colour = family_id_group)) +
    geom_line(alpha = 0.6) +
    geom_point(size = 1.2, alpha = 0.7) +
    facet_wrap(~ treatment_group) +
    scale_colour_manual(values = fam_colours, name = "Family") +
    labs(x = "Time (h)", y = "Raw fluorescence (RFU)") +
    theme_classic(base_size = 12) +
    theme(strip.background = element_blank(),
          strip.text       = element_text(face = "bold"))

  p_raw_by_treatment
  ggsave(file.path(fig_dir, "raw_individual_by_treatment.png"),
         p_raw_by_treatment, width = 10, height = 5)
}

6.6 Excluded samples

Wells flagged exclude_from_analysis = TRUE appear in the raw fluorescence plots above but are omitted from all analyses that follow.

excluded_wells <- dat %>%
  filter(!is_blank, exclude_from_analysis) %>%
  mutate(
    sample = if_else(
      !is.na(sample_id_group) & trimws(as.character(sample_id_group)) != "",
      as.character(sample_id_group),
      paste(plate_id, well_id, sep = "_")
    )
  ) %>%
  select(plate_id, well_id, sample, family_id_group, treatment_group,
         any_of("exclude_reason")) %>%
  distinct() %>%
  arrange(plate_id, well_id)

if (nrow(excluded_wells) > 0) {
  col_names <- c("Plate", "Well", "Sample", "Family", "Treatment")
  if ("exclude_reason" %in% names(excluded_wells))
    col_names <- c(col_names, "Reason")
  cat(knitr::kable(excluded_wells, col.names = col_names), sep = "\n")
} else {
  cat("No wells are excluded from analysis.\n")
}
Plate Well Sample Family Treatment Reason
a E11 134 NA NA Duplicate sample ID
a F1 139 NA NA Duplicate sample ID
a F2 138 NA NA Duplicate sample ID
a F5 140 NA NA Duplicate sample ID
a F6 136 NA NA Duplicate sample ID
a F7 137 NA NA Duplicate sample ID
a F8 135 NA NA Duplicate sample ID
b A11 140 NA NA Duplicate sample ID
b A2 139 NA NA Duplicate sample ID
b A3 136 NA NA Duplicate sample ID
b A4 134 NA NA Duplicate sample ID
b A7 137 NA NA Duplicate sample ID
b A8 138 NA NA Duplicate sample ID
b A9 135 NA NA Duplicate sample ID
c D12 138 NA NA Duplicate sample ID
c E5 139 NA NA Duplicate sample ID
c E6 122 NA NA Duplicate sample ID
c F2 137 NA NA Duplicate sample ID
c F3 135 NA NA Duplicate sample ID
c F4 134 NA NA Duplicate sample ID
c F5 136 NA NA Duplicate sample ID
c F6 140 NA NA Duplicate sample ID
c H10 139 NA NA Duplicate sample ID
d A11 135 NA NA Duplicate sample ID
d A3 137 NA NA Duplicate sample ID
d A4 138 NA NA Duplicate sample ID
d A5 136 NA NA Duplicate sample ID
d A6 134 NA NA Duplicate sample ID
d A7 140 NA NA Duplicate sample ID
e D12 140 NA NA Duplicate sample ID
e E1 136 NA NA Duplicate sample ID
e E10 139 NA NA Duplicate sample ID
e E5 138 NA NA Duplicate sample ID
e F2 135 NA NA Duplicate sample ID
e F3 134 NA NA Duplicate sample ID
e F5 137 NA NA Duplicate sample ID
e H10 139 NA NA Duplicate sample ID
e H9 138 NA NA Duplicate sample ID
f A10 136 NA NA Duplicate sample ID
f A11 134 NA NA Duplicate sample ID
f A12 135 NA NA Duplicate sample ID
f A6 137 NA NA Duplicate sample ID
f A8 140 NA NA Duplicate sample ID
g D10 140 NA NA Duplicate sample ID
g E1 139 NA NA Duplicate sample ID
g E11 135 NA NA Duplicate sample ID
g E6 136 NA NA Duplicate sample ID
g E8 138 NA NA Duplicate sample ID
g F3 137 NA NA Duplicate sample ID
g F4 134 NA NA Duplicate sample ID
g H7 138 NA NA Duplicate sample ID
g H8 139 NA NA Duplicate sample ID
h A11 137 NA NA Duplicate sample ID
h A6 140 NA NA Duplicate sample ID
h A7 134 NA NA Duplicate sample ID
h A8 136 NA NA Duplicate sample ID
h A9 135 NA NA Duplicate sample ID
i D12 139 NA NA Duplicate sample ID
i D3 138 NA NA Duplicate sample ID
i D4 140 NA NA Duplicate sample ID
i E2 135 NA NA Duplicate sample ID
i F2 137 NA NA Duplicate sample ID
i F3 134 NA NA Duplicate sample ID
i F4 136 NA NA Duplicate sample ID
i H6 136 NA NA Duplicate sample ID
i H7 139 NA NA Duplicate sample ID
i H9 138 NA NA Duplicate sample ID
j A3 135 NA NA Duplicate sample ID
j A6 137 NA NA Duplicate sample ID
j A8 140 NA NA Duplicate sample ID
j A9 134 NA NA Duplicate sample ID

7 Blank Correction via Fold-Change Normalization

T0 is the earliest timepoint present in the dataset (not necessarily 0 hr). Sample fold-change is expressed relative to each individual’s T0 reading, resolved by sample_id_group when that column is populated — allowing the same animal to be tracked across plates — or by plate_id + well_id when no sample IDs exist (backward-compatible with single-plate, multi-timepoint designs). Blank fold-change is the per-plate mean blank RFU at each timepoint divided by the pooled mean blank RFU at T0. Subtracting blank fold-change from sample fold-change removes background fluorescence drift; all samples start at exactly 0 at T0 by construction.

7.1 Step 1 – Identify T0 and compute per-sample fold-change

# T0 = earliest timepoint present in the dataset
t0_time <- min(dat$time_hr[is.finite(dat$time_hr)], na.rm = TRUE)
message("T0 timepoint: ", t0_time, " hr")

# T0 reference value per individual.
# Resolved by sample_id_group (cross-plate tracking) when available;
# falls back to plate+well for layouts without explicit sample IDs.
t0_all <- dat %>%
  filter(time_hr == t0_time, !is_blank, is.finite(value)) %>%
  mutate(sample_key = if_else(
    !is.na(sample_id_group) & trimws(as.character(sample_id_group)) != "",
    as.character(sample_id_group),
    paste(plate_id, well_id, sep = "_")
  )) %>%
  group_by(sample_key) %>%
  summarise(value_t0 = mean(value, na.rm = TRUE), .groups = "drop")

dat_fc <- dat %>%
  mutate(sample_key = if_else(
    !is_blank &
      !is.na(sample_id_group) & trimws(as.character(sample_id_group)) != "",
    as.character(sample_id_group),
    paste(plate_id, well_id, sep = "_")
  )) %>%
  left_join(t0_all, by = "sample_key") %>%
  mutate(fold_change = if_else(
    !is_blank & is.finite(value_t0) & value_t0 > 0,
    value / value_t0,
    NA_real_
  ))

str(dat_fc)
tibble [960 × 19] (S3: tbl_df/tbl/data.frame)
 $ row_id               : chr [1:960] "A" "A" "A" "A" ...
 $ col_id               : int [1:960] 1 2 3 4 5 6 7 8 9 10 ...
 $ well_id              : chr [1:960] "A1" "A2" "A3" "A4" ...
 $ value                : num [1:960] 727 599 672 685 639 720 683 686 794 975 ...
 $ plate_id             : chr [1:960] "a" "a" "a" "a" ...
 $ time_hr              : num [1:960] 1 1 1 1 1 1 1 1 1 1 ...
 $ is_blank             : logi [1:960] FALSE FALSE FALSE FALSE FALSE FALSE ...
 $ exclude_from_analysis: logi [1:960] FALSE FALSE FALSE FALSE FALSE FALSE ...
 $ exclude_reason       : chr [1:960] NA NA NA NA ...
 $ family_id_group      : chr [1:960] "1" "9" "9" "9" ...
 $ sample_id_group      : chr [1:960] "47" "38" "52" "57" ...
 $ treatment_group      : chr [1:960] NA NA NA NA ...
 $ weight_g_measurement : num [1:960] NA NA NA NA NA NA NA NA NA NA ...
 $ width_mm_measurement : num [1:960] NA NA NA NA NA NA NA NA NA NA ...
 $ length_mm_measurement: num [1:960] NA NA NA NA NA NA NA NA NA NA ...
 $ area_mm2_measurement : num [1:960] 1346 648 577 682 372 ...
 $ sample_key           : chr [1:960] "47" "38" "52" "57" ...
 $ value_t0             : num [1:960] 727 599 672 685 639 720 683 686 794 975 ...
 $ fold_change          : num [1:960] 1 1 1 1 1 1 1 1 1 1 ...

7.2 Step 2 – Blank fold-change reference per plate per timepoint

# Pooled mean blank RFU at T0 across all T0 plates
mean_blank_t0 <- dat %>%
  filter(is_blank, time_hr == t0_time, is.finite(value)) %>%
  pull(value) %>%
  mean(na.rm = TRUE)

if (!is.finite(mean_blank_t0))
  message("No blank readings found at T0 (", t0_time,
          " hr); blank correction will produce NA.")

# Per-plate per-timepoint mean blank expressed as fold-change relative to T0
blank_fc_ref <- dat %>%
  filter(is_blank, is.finite(value)) %>%
  group_by(plate_id, time_hr) %>%
  summarise(mean_blank_rfu = mean(value, na.rm = TRUE), .groups = "drop") %>%
  mutate(mean_blank_fc = mean_blank_rfu / mean_blank_t0)

str(blank_fc_ref)
tibble [10 × 4] (S3: tbl_df/tbl/data.frame)
 $ plate_id      : chr [1:10] "a" "b" "c" "d" ...
 $ time_hr       : num [1:10] 1 1 2 2 3 3 5 5 7 7
 $ mean_blank_rfu: num [1:10] 590 545 622 621 602 604 638 687 955 916
 $ mean_blank_fc : num [1:10] 1.04 0.96 1.1 1.09 1.06 ...

7.3 Step 3 – Subtract blank fold-change from sample fold-change

samples <- dat_fc %>%
  filter(!is_blank, !exclude_from_analysis) %>%
  mutate(
    trace_id = if_else(
      !is.na(sample_id_group) & trimws(as.character(sample_id_group)) != "",
      as.character(sample_id_group),
      paste(plate_id, well_id, sep = "_")
    )
  ) %>%
  left_join(blank_fc_ref, by = c("plate_id", "time_hr")) %>%
  mutate(corrected_fc = fold_change - mean_blank_fc)

str(samples)
tibble [879 × 23] (S3: tbl_df/tbl/data.frame)
 $ row_id               : chr [1:879] "A" "A" "A" "A" ...
 $ col_id               : int [1:879] 1 2 3 4 5 6 7 8 9 10 ...
 $ well_id              : chr [1:879] "A1" "A2" "A3" "A4" ...
 $ value                : num [1:879] 727 599 672 685 639 720 683 686 794 975 ...
 $ plate_id             : chr [1:879] "a" "a" "a" "a" ...
 $ time_hr              : num [1:879] 1 1 1 1 1 1 1 1 1 1 ...
 $ is_blank             : logi [1:879] FALSE FALSE FALSE FALSE FALSE FALSE ...
 $ exclude_from_analysis: logi [1:879] FALSE FALSE FALSE FALSE FALSE FALSE ...
 $ exclude_reason       : chr [1:879] NA NA NA NA ...
 $ family_id_group      : chr [1:879] "1" "9" "9" "9" ...
 $ sample_id_group      : chr [1:879] "47" "38" "52" "57" ...
 $ treatment_group      : chr [1:879] NA NA NA NA ...
 $ weight_g_measurement : num [1:879] NA NA NA NA NA NA NA NA NA NA ...
 $ width_mm_measurement : num [1:879] NA NA NA NA NA NA NA NA NA NA ...
 $ length_mm_measurement: num [1:879] NA NA NA NA NA NA NA NA NA NA ...
 $ area_mm2_measurement : num [1:879] 1346 648 577 682 372 ...
 $ sample_key           : chr [1:879] "47" "38" "52" "57" ...
 $ value_t0             : num [1:879] 727 599 672 685 639 720 683 686 794 975 ...
 $ fold_change          : num [1:879] 1 1 1 1 1 1 1 1 1 1 ...
 $ trace_id             : chr [1:879] "47" "38" "52" "57" ...
 $ mean_blank_rfu       : num [1:879] 590 590 590 590 590 590 590 590 590 590 ...
 $ mean_blank_fc        : num [1:879] 1.04 1.04 1.04 1.04 1.04 ...
 $ corrected_fc         : num [1:879] -0.0396 -0.0396 -0.0396 -0.0396 -0.0396 ...

8 Blank-Corrected Fold-Change

8.1 Mean by family

bc_fc_summary <- samples %>%
  filter(!is.na(family_id_group), !exclude_from_analysis) %>%
  group_by(family_id_group, treatment_group, time_hr) %>%
  summarise(
    mean_val = mean(corrected_fc, na.rm = TRUE),
    se_val   = sd(corrected_fc, na.rm = TRUE) /
      sqrt(sum(!is.na(corrected_fc))),
    n        = sum(!is.na(corrected_fc)),
    .groups  = "drop"
  ) %>%
  mutate(group_var = if (has_trt)
    paste(family_id_group, treatment_group, sep = ".")
  else
    family_id_group)

p_bc_fc_mean <- ggplot(bc_fc_summary,
    aes(x = time_hr, y = mean_val,
        colour = family_id_group,
        group = group_var)) +
  geom_ribbon(aes(ymin = mean_val - se_val,
                  ymax = mean_val + se_val,
                  fill = family_id_group),
              alpha = 0.15, colour = NA) +
  geom_line(
    mapping   = if (has_trt) aes(linetype = treatment_group) else NULL,
    linewidth = 1) +
  geom_point(size = 2) +
  scale_colour_manual(values = fam_colours, name = "Family") +
  scale_fill_manual(values = fam_colours, name = "Family") +
  labs(x = "Time (h)",
       y = "Mean blank-corrected fold-change (± SE)") +
  theme_classic(base_size = 13) +
  if (has_trt) scale_linetype_manual(values = trt_linetypes, name = "Treatment") else NULL

p_bc_fc_mean

ggsave(file.path(fig_dir, "blank_corrected_fc_mean_by_family.png"),
       p_bc_fc_mean, width = 8, height = 5)

8.2 Individual traces by family

p_bc_fc_by_family <- samples %>%
  filter(!is.na(family_id_group)) %>%
  ggplot(aes(x = time_hr, y = corrected_fc, group = trace_id,
             colour = .data[[if (has_trt) "treatment_group" else "family_id_group"]])) +
  geom_line(alpha = 0.6) +
  geom_point(size = 1.2, alpha = 0.7) +
  facet_wrap(~ family_id_group) +
  scale_colour_manual(
    values = if (has_trt) trt_colours else fam_colours,
    name   = if (has_trt) "Treatment" else "Family") +
  labs(x = "Time (h)", y = "Blank-corrected fold-change") +
  theme_classic(base_size = 12) +
  theme(strip.background = element_blank(),
        strip.text       = element_text(face = "bold"))

p_bc_fc_by_family

ggsave(file.path(fig_dir, "blank_corrected_fc_by_family.png"),
       p_bc_fc_by_family, width = 10, height = 5)

8.3 Individual blank-corrected fold-change traces by treatment

if (has_trt) {
  p_bc_fc_by_treatment <- samples %>%
    ggplot(aes(x = time_hr, y = corrected_fc,
               group = trace_id, colour = family_id_group)) +
    geom_line(alpha = 0.6) +
    geom_point(size = 1.2, alpha = 0.7) +
    facet_wrap(~ treatment_group) +
    scale_colour_manual(values = fam_colours, name = "Family") +
    labs(x = "Time (h)", y = "Blank-corrected fold-change") +
    theme_classic(base_size = 12) +
    theme(strip.background = element_blank(),
          strip.text       = element_text(face = "bold"))

  p_bc_fc_by_treatment
  ggsave(file.path(fig_dir, "blank_corrected_fc_by_treatment.png"),
         p_bc_fc_by_treatment, width = 10, height = 5)
}

9 Metabolism (Size-Normalised Fold-Change)

Blank-corrected fold-change divided by each active measurement column. This is “metabolism” as defined in Huffmyer et al.

if (length(active_meas_cols) == 0) {
  message("No active measurement columns: skipping metabolism calculation.")
  metabolism_df <- tibble()
} else {
  metabolism_df <- samples
  for (mc in active_meas_cols) {
    out_col <- paste0("metabolism_per_", mc)
    metabolism_df <- metabolism_df %>%
      mutate(!!out_col := if_else(
        is.finite(.data[[mc]]) & .data[[mc]] > 0 &
          is.finite(corrected_fc),
        corrected_fc / .data[[mc]],
        NA_real_
      ))
  }
}

str(metabolism_df)
tibble [879 × 24] (S3: tbl_df/tbl/data.frame)
 $ row_id                             : chr [1:879] "A" "A" "A" "A" ...
 $ col_id                             : int [1:879] 1 2 3 4 5 6 7 8 9 10 ...
 $ well_id                            : chr [1:879] "A1" "A2" "A3" "A4" ...
 $ value                              : num [1:879] 727 599 672 685 639 720 683 686 794 975 ...
 $ plate_id                           : chr [1:879] "a" "a" "a" "a" ...
 $ time_hr                            : num [1:879] 1 1 1 1 1 1 1 1 1 1 ...
 $ is_blank                           : logi [1:879] FALSE FALSE FALSE FALSE FALSE FALSE ...
 $ exclude_from_analysis              : logi [1:879] FALSE FALSE FALSE FALSE FALSE FALSE ...
 $ exclude_reason                     : chr [1:879] NA NA NA NA ...
 $ family_id_group                    : chr [1:879] "1" "9" "9" "9" ...
 $ sample_id_group                    : chr [1:879] "47" "38" "52" "57" ...
 $ treatment_group                    : chr [1:879] NA NA NA NA ...
 $ weight_g_measurement               : num [1:879] NA NA NA NA NA NA NA NA NA NA ...
 $ width_mm_measurement               : num [1:879] NA NA NA NA NA NA NA NA NA NA ...
 $ length_mm_measurement              : num [1:879] NA NA NA NA NA NA NA NA NA NA ...
 $ area_mm2_measurement               : num [1:879] 1346 648 577 682 372 ...
 $ sample_key                         : chr [1:879] "47" "38" "52" "57" ...
 $ value_t0                           : num [1:879] 727 599 672 685 639 720 683 686 794 975 ...
 $ fold_change                        : num [1:879] 1 1 1 1 1 1 1 1 1 1 ...
 $ trace_id                           : chr [1:879] "47" "38" "52" "57" ...
 $ mean_blank_rfu                     : num [1:879] 590 590 590 590 590 590 590 590 590 590 ...
 $ mean_blank_fc                      : num [1:879] 1.04 1.04 1.04 1.04 1.04 ...
 $ corrected_fc                       : num [1:879] -0.0396 -0.0396 -0.0396 -0.0396 -0.0396 ...
 $ metabolism_per_area_mm2_measurement: num [1:879] -2.95e-05 -6.12e-05 -6.87e-05 -5.81e-05 -1.06e-04 ...

9.1 Mean metabolism by family

if (nrow(metabolism_df) > 0) {

  metab_cols <- paste0("metabolism_per_", active_meas_cols)

  for (col in metab_cols) {
    if (!col %in% names(metabolism_df)) next
    mc_label <- str_remove(col, "^metabolism_per_")

    metab_summary <- metabolism_df %>%
      filter(!is.na(family_id_group), !exclude_from_analysis) %>%
      group_by(family_id_group, treatment_group, time_hr) %>%
      summarise(
        mean_val = mean(.data[[col]], na.rm = TRUE),
        se_val   = sd(.data[[col]], na.rm = TRUE) /
          sqrt(sum(!is.na(.data[[col]]))),
        .groups  = "drop"
      ) %>%
      mutate(group_var = if (has_trt)
        paste(family_id_group, treatment_group, sep = ".")
      else
        family_id_group)

    p_metab_mean <- ggplot(metab_summary,
        aes(x = time_hr, y = mean_val,
            colour = family_id_group,
            group = group_var)) +
      geom_ribbon(aes(ymin = mean_val - se_val,
                      ymax = mean_val + se_val,
                      fill = family_id_group),
                  alpha = 0.15, colour = NA) +
      geom_line(
        mapping   = if (has_trt) aes(linetype = treatment_group) else NULL,
        linewidth = 1) +
      geom_point(size = 2) +
      scale_colour_manual(values = fam_colours, name = "Family") +
      scale_fill_manual(values = fam_colours, name = "Family") +
      labs(x = "Time (h)",
           y = paste0(metabolism_y_label(col), " (± SE)")) +
      theme_classic(base_size = 13) +
      if (has_trt) scale_linetype_manual(values = trt_linetypes, name = "Treatment") else NULL

    print(p_metab_mean)
    ggsave(
      file.path(fig_dir,
                paste0("metabolism_mean_", mc_label, ".png")),
      p_metab_mean, width = 8, height = 5)
  }
}

9.2 Individual metabolism traces by family

if (nrow(metabolism_df) > 0) {

  for (col in metab_cols) {
    if (!col %in% names(metabolism_df)) next
    mc_label <- str_remove(col, "^metabolism_per_")

    p_metab_by_family <- metabolism_df %>%
      filter(!is.na(family_id_group)) %>%
      ggplot(aes(x = time_hr, y = .data[[col]], group = trace_id,
                 colour = .data[[if (has_trt) "treatment_group" else "family_id_group"]])) +
      geom_line(alpha = 0.6) +
      geom_point(size = 1.2, alpha = 0.7) +
      facet_wrap(~ family_id_group) +
      scale_colour_manual(
        values = if (has_trt) trt_colours else fam_colours,
        name   = if (has_trt) "Treatment" else "Family") +
      labs(x = "Time (h)", y = metabolism_y_label(col)) +
      theme_classic(base_size = 12) +
      theme(strip.background = element_blank(),
            strip.text       = element_text(face = "bold"))

    print(p_metab_by_family)
    ggsave(
      file.path(fig_dir,
                paste0("metabolism_individual_", mc_label, "_by_family.png")),
      p_metab_by_family, width = 10, height = 5)

    if (has_trt) {
      p_metab_by_treatment <- ggplot(metabolism_df,
          aes(x = time_hr, y = .data[[col]],
              group = trace_id, colour = family_id_group)) +
        geom_line(alpha = 0.6) +
        geom_point(size = 1.2, alpha = 0.7) +
        facet_wrap(~ treatment_group) +
        scale_colour_manual(values = fam_colours, name = "Family") +
        labs(x = "Time (h)", y = metabolism_y_label(col)) +
        theme_classic(base_size = 12) +
        theme(strip.background = element_blank(),
              strip.text       = element_text(face = "bold"))

      print(p_metab_by_treatment)
      ggsave(
        file.path(fig_dir,
                  paste0("metabolism_individual_", mc_label, "_by_treatment.png")),
        p_metab_by_treatment, width = 10, height = 5)
    }
  }
}

10 Time-Series Statistical Analysis

Linear mixed effects models test the effect of experimental variables on metabolism over time. Individual (sample_id_group) is included as a random intercept to account for repeated measures across timepoints. Type III ANOVA with Satterthwaite’s approximation (lmerTest) assesses significance; post-hoc pairwise comparisons use estimated marginal means (emmeans, Tukey adjustment).

run_ts_stats <- function(df, value_col) {
  has_family    <- "family_id_group" %in% names(df) &&
    length(unique(na.omit(df$family_id_group))) > 1
  has_treatment <- "treatment_group" %in% names(df) &&
    length(unique(na.omit(df$treatment_group))) > 1

  if (!has_family && !has_treatment) return(NULL)

  df <- df %>%
    filter(is.finite(.data[[value_col]]), is.finite(time_hr)) %>%
    mutate(
      time_f     = factor(time_hr),
      individual = factor(trace_id)
    )

  if (nrow(df) == 0) return(NULL)

  if (has_family)    df <- df %>% mutate(family    = factor(family_id_group))
  if (has_treatment) df <- df %>% mutate(treatment = factor(treatment_group))

  if (has_family    && length(unique(na.omit(df$family)))    < 2) return(NULL)
  if (has_treatment && length(unique(na.omit(df$treatment))) < 2) return(NULL)

  fixed <- if (has_family && has_treatment) {
    paste0(value_col, " ~ time_f * family * treatment")
  } else if (has_family) {
    paste0(value_col, " ~ time_f * family")
  } else {
    paste0(value_col, " ~ time_f * treatment")
  }

  model <- lmer(
    as.formula(paste0(fixed, " + (1 | individual)")),
    data = df
  )

  anova_res <- anova(model, type = 3, ddf = "Satterthwaite")

  # Pairwise comparisons of group combinations at each timepoint
  emm_spec <- if (has_family && has_treatment) {
    ~ family * treatment | time_f
  } else if (has_family) {
    ~ family | time_f
  } else {
    ~ treatment | time_f
  }

  emm       <- emmeans(model, emm_spec)
  pairs_res <- as.data.frame(pairs(emm, adjust = "tukey"))

  # Main-effect marginal means (collapsed across time)
  emm_main <- if (has_family && has_treatment) {
    emmeans(model, ~ family * treatment)
  } else if (has_family) {
    emmeans(model, ~ family)
  } else {
    emmeans(model, ~ treatment)
  }

  pairs_main <- as.data.frame(pairs(emm_main, adjust = "tukey"))

  list(
    model         = model,
    anova         = anova_res,
    pairs_by_time = pairs_res,
    pairs_main    = pairs_main,
    has_family    = has_family,
    has_treatment = has_treatment
  )
}

ts_stats <- list()
if (nrow(metabolism_df) > 0) {
  for (mc in active_meas_cols) {
    col <- paste0("metabolism_per_", mc)
    if (col %in% names(metabolism_df))
      ts_stats[[col]] <- run_ts_stats(metabolism_df, col)
  }
}

10.1 Results

for (col in names(ts_stats)) {
  res <- ts_stats[[col]]
  if (is.null(res)) next

  cat("\n\n### Metric:", col, "\n\n")

  cat("**Type III ANOVA (Satterthwaite approximation):**\n\n")
  cat(knitr::kable(as.data.frame(res$anova), digits = 4, format = "pipe"), sep = "\n")
  cat("\n")

  cat("**Marginal means – main effects (collapsed across time):**\n\n")
  cat(knitr::kable(as.data.frame(res$pairs_main), digits = 4, format = "pipe"), sep = "\n")
  cat("\n")

  cat("**Pairwise comparisons by timepoint (Tukey):**\n\n")
  cat(knitr::kable(as.data.frame(res$pairs_by_time), digits = 4, format = "pipe"), sep = "\n")
  cat("\n")
}

10.1.1 Metric: metabolism_per_area_mm2_measurement

Type III ANOVA (Satterthwaite approximation):

Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
time_f 0.0065 0.0016 4 466.5627 111.4948 0.0000
family 0.0001 0.0000 8 116.4978 0.7997 0.6040
time_f:family 0.0004 0.0000 32 466.5881 0.9526 0.5441

Marginal means – main effects (collapsed across time):

contrast estimate SE df t.ratio p.value
family1 - family10 0.0000 0.0015 116.9082 0.0206 1.0000
family1 - family2 0.0004 0.0015 116.9082 0.2823 1.0000
family1 - family3 -0.0018 0.0012 116.9082 -1.4421 0.8790
family1 - family5 0.0000 0.0012 116.9082 -0.0184 1.0000
family1 - family6 0.0024 0.0024 116.9082 1.0081 0.9844
family1 - family7 -0.0011 0.0014 116.9082 -0.7674 0.9975
family1 - family8 -0.0003 0.0015 116.9082 -0.1924 1.0000
family1 - family9 0.0005 0.0012 117.1840 0.4427 1.0000
family10 - family2 0.0004 0.0017 116.9082 0.2267 1.0000
family10 - family3 -0.0018 0.0015 116.9082 -1.1980 0.9555
family10 - family5 -0.0001 0.0015 116.9082 -0.0355 1.0000
family10 - family6 0.0024 0.0025 116.9082 0.9361 0.9904
family10 - family7 -0.0011 0.0016 116.9082 -0.6689 0.9991
family10 - family8 -0.0003 0.0017 116.9082 -0.1844 1.0000
family10 - family9 0.0005 0.0015 117.0920 0.3410 1.0000
family2 - family3 -0.0022 0.0015 116.9082 -1.4598 0.8715
family2 - family5 -0.0004 0.0015 116.9082 -0.2973 1.0000
family2 - family6 0.0020 0.0025 116.9082 0.7821 0.9972
family2 - family7 -0.0015 0.0016 116.9082 -0.9100 0.9920
family2 - family8 -0.0007 0.0017 116.9082 -0.4111 1.0000
family2 - family9 0.0001 0.0015 117.0920 0.0793 1.0000
family3 - family5 0.0017 0.0012 116.9082 1.4237 0.8865
family3 - family6 0.0042 0.0024 116.9082 1.7447 0.7177
family3 - family7 0.0007 0.0014 116.9082 0.5126 0.9999
family3 - family8 0.0015 0.0015 116.9082 0.9851 0.9866
family3 - family9 0.0023 0.0012 117.1840 1.8839 0.6261
family5 - family6 0.0024 0.0024 116.9082 1.0175 0.9835
family5 - family7 -0.0010 0.0014 116.9082 -0.7512 0.9979
family5 - family8 -0.0003 0.0015 116.9082 -0.1774 1.0000
family5 - family9 0.0006 0.0012 117.1840 0.4611 0.9999
family6 - family7 -0.0035 0.0025 116.9082 -1.4014 0.8952
family6 - family8 -0.0027 0.0025 116.9082 -1.0614 0.9785
family6 - family9 -0.0019 0.0024 116.9801 -0.7818 0.9972
family7 - family8 0.0008 0.0016 116.9082 0.4729 0.9999
family7 - family9 0.0016 0.0014 117.1255 1.1601 0.9631
family8 - family9 0.0008 0.0015 117.0920 0.5538 0.9998

Pairwise comparisons by timepoint (Tukey):

contrast time_f estimate SE df t.ratio p.value
family1 - family10 1 0.0000 0.0020 323.8818 -0.0038 1.0000
family1 - family2 1 -0.0001 0.0020 323.8818 -0.0383 1.0000
family1 - family3 1 0.0000 0.0016 323.8818 -0.0060 1.0000
family1 - family5 1 0.0000 0.0016 323.8818 0.0149 1.0000
family1 - family6 1 -0.0001 0.0032 323.8818 -0.0272 1.0000
family1 - family7 1 0.0000 0.0018 323.8818 0.0162 1.0000
family1 - family8 1 -0.0001 0.0020 323.8818 -0.0358 1.0000
family1 - family9 1 0.0000 0.0016 323.8818 0.0194 1.0000
family10 - family2 1 -0.0001 0.0023 323.8818 -0.0298 1.0000
family10 - family3 1 0.0000 0.0020 323.8818 -0.0011 1.0000
family10 - family5 1 0.0000 0.0020 323.8818 0.0160 1.0000
family10 - family6 1 -0.0001 0.0034 323.8818 -0.0234 1.0000
family10 - family7 1 0.0000 0.0022 323.8818 0.0173 1.0000
family10 - family8 1 -0.0001 0.0023 323.8818 -0.0277 1.0000
family10 - family9 1 0.0000 0.0020 323.8818 0.0196 1.0000
family2 - family3 1 0.0001 0.0020 323.8818 0.0334 1.0000
family2 - family5 1 0.0001 0.0020 323.8818 0.0505 1.0000
family2 - family6 1 0.0000 0.0034 323.8818 -0.0031 1.0000
family2 - family7 1 0.0001 0.0022 323.8818 0.0490 1.0000
family2 - family8 1 0.0000 0.0023 323.8818 0.0022 1.0000
family2 - family9 1 0.0001 0.0020 323.8818 0.0541 1.0000
family3 - family5 1 0.0000 0.0016 323.8818 0.0210 1.0000
family3 - family6 1 -0.0001 0.0032 323.8818 -0.0242 1.0000
family3 - family7 1 0.0000 0.0018 323.8818 0.0216 1.0000
family3 - family8 1 -0.0001 0.0020 323.8818 -0.0309 1.0000
family3 - family9 1 0.0000 0.0016 323.8818 0.0254 1.0000
family5 - family6 1 -0.0001 0.0032 323.8818 -0.0349 1.0000
family5 - family7 1 0.0000 0.0018 323.8818 0.0030 1.0000
family5 - family8 1 -0.0001 0.0020 323.8818 -0.0480 1.0000
family5 - family9 1 0.0000 0.0016 323.8818 0.0044 1.0000
family6 - family7 1 0.0001 0.0033 323.8818 0.0353 1.0000
family6 - family8 1 0.0000 0.0034 323.8818 0.0046 1.0000
family6 - family9 1 0.0001 0.0032 323.8818 0.0371 1.0000
family7 - family8 1 -0.0001 0.0022 323.8818 -0.0467 1.0000
family7 - family9 1 0.0000 0.0018 323.8818 0.0010 1.0000
family8 - family9 1 0.0001 0.0020 323.8818 0.0516 1.0000
family1 - family10 2 -0.0010 0.0020 323.8818 -0.5077 0.9999
family1 - family2 2 0.0003 0.0020 323.8818 0.1591 1.0000
family1 - family3 2 -0.0005 0.0016 323.8818 -0.2867 1.0000
family1 - family5 2 0.0001 0.0016 323.8818 0.0432 1.0000
family1 - family6 2 0.0006 0.0032 323.8818 0.2021 1.0000
family1 - family7 2 -0.0008 0.0018 323.8818 -0.4310 1.0000
family1 - family8 2 0.0001 0.0020 323.8818 0.0276 1.0000
family1 - family9 2 0.0000 0.0016 330.9462 -0.0062 1.0000
family10 - family2 2 0.0013 0.0023 323.8818 0.5774 0.9997
family10 - family3 2 0.0005 0.0020 323.8818 0.2735 1.0000
family10 - family5 2 0.0011 0.0020 323.8818 0.5429 0.9998
family10 - family6 2 0.0017 0.0034 323.8818 0.4888 0.9999
family10 - family7 2 0.0002 0.0022 323.8818 0.1024 1.0000
family10 - family8 2 0.0011 0.0023 323.8818 0.4635 0.9999
family10 - family9 2 0.0010 0.0020 328.6037 0.4996 0.9999
family2 - family3 2 -0.0008 0.0020 323.8818 -0.3932 1.0000
family2 - family5 2 -0.0002 0.0020 323.8818 -0.1238 1.0000
family2 - family6 2 0.0003 0.0034 323.8818 0.0965 1.0000
family2 - family7 2 -0.0011 0.0022 323.8818 -0.5116 0.9999
family2 - family8 2 -0.0003 0.0023 323.8818 -0.1139 1.0000
family2 - family9 2 -0.0003 0.0020 328.6037 -0.1633 1.0000
family3 - family5 2 0.0005 0.0016 323.8818 0.3299 1.0000
family3 - family6 2 0.0011 0.0032 323.8818 0.3486 1.0000
family3 - family7 2 -0.0003 0.0018 323.8818 -0.1765 1.0000
family3 - family8 2 0.0005 0.0020 323.8818 0.2617 1.0000
family3 - family9 2 0.0005 0.0016 330.9462 0.2781 1.0000
family5 - family6 2 0.0006 0.0032 323.8818 0.1801 1.0000
family5 - family7 2 -0.0009 0.0018 323.8818 -0.4693 0.9999
family5 - family8 2 0.0000 0.0020 323.8818 -0.0077 1.0000
family5 - family9 2 -0.0001 0.0016 330.9462 -0.0491 1.0000
family6 - family7 2 -0.0014 0.0033 323.8818 -0.4351 1.0000
family6 - family8 2 -0.0006 0.0034 323.8818 -0.1739 1.0000
family6 - family9 2 -0.0007 0.0032 325.7353 -0.2049 1.0000
family7 - family8 2 0.0008 0.0022 323.8818 0.3905 1.0000
family7 - family9 2 0.0008 0.0018 329.4570 0.4225 1.0000
family8 - family9 2 -0.0001 0.0020 328.6037 -0.0325 1.0000
family1 - family10 3 -0.0010 0.0020 323.8818 -0.5041 0.9999
family1 - family2 3 0.0003 0.0020 323.8818 0.1437 1.0000
family1 - family3 3 -0.0015 0.0016 323.8818 -0.8920 0.9933
family1 - family5 3 -0.0004 0.0016 323.8818 -0.2511 1.0000
family1 - family6 3 0.0022 0.0032 323.8818 0.6766 0.9990
family1 - family7 3 -0.0011 0.0018 323.8818 -0.5906 0.9996
family1 - family8 3 -0.0007 0.0020 323.8818 -0.3366 1.0000
family1 - family9 3 0.0006 0.0016 323.8818 0.3796 1.0000
family10 - family2 3 0.0013 0.0023 323.8818 0.5610 0.9998
family10 - family3 3 -0.0004 0.0020 323.8818 -0.2243 1.0000
family10 - family5 3 0.0006 0.0020 323.8818 0.2991 1.0000
family10 - family6 3 0.0032 0.0034 323.8818 0.9329 0.9909
family10 - family7 3 -0.0001 0.0022 323.8818 -0.0361 1.0000
family10 - family8 3 0.0003 0.0023 323.8818 0.1451 1.0000
family10 - family9 3 0.0016 0.0020 323.8818 0.8140 0.9964
family2 - family3 3 -0.0017 0.0020 323.8818 -0.8720 0.9942
family2 - family5 3 -0.0007 0.0020 323.8818 -0.3487 1.0000
family2 - family6 3 0.0019 0.0034 323.8818 0.5518 0.9998
family2 - family7 3 -0.0014 0.0022 323.8818 -0.6325 0.9994
family2 - family8 3 -0.0010 0.0023 323.8818 -0.4159 1.0000
family2 - family9 3 0.0003 0.0020 323.8818 0.1663 1.0000
family3 - family5 3 0.0010 0.0016 323.8818 0.6410 0.9994
family3 - family6 3 0.0036 0.0032 323.8818 1.1322 0.9690
family3 - family7 3 0.0004 0.0018 323.8818 0.2011 1.0000
family3 - family8 3 0.0008 0.0020 323.8818 0.3918 1.0000
family3 - family9 3 0.0021 0.0016 323.8818 1.2716 0.9389
family5 - family6 3 0.0026 0.0032 323.8818 0.8048 0.9967
family5 - family7 3 -0.0007 0.0018 323.8818 -0.3678 1.0000
family5 - family8 3 -0.0003 0.0020 323.8818 -0.1316 1.0000
family5 - family9 3 0.0010 0.0016 323.8818 0.6306 0.9994
family6 - family7 3 -0.0032 0.0033 323.8818 -0.9825 0.9873
family6 - family8 3 -0.0028 0.0034 323.8818 -0.8344 0.9957
family6 - family9 3 -0.0015 0.0032 323.8818 -0.4827 0.9999
family7 - family8 3 0.0004 0.0022 323.8818 0.1904 1.0000
family7 - family9 3 0.0017 0.0018 323.8818 0.9276 0.9913
family8 - family9 3 0.0013 0.0020 323.8818 0.6465 0.9993
family1 - family10 5 -0.0005 0.0020 323.8818 -0.2538 1.0000
family1 - family2 5 0.0004 0.0020 323.8818 0.1958 1.0000
family1 - family3 5 -0.0028 0.0016 323.8818 -1.7298 0.7277
family1 - family5 5 -0.0004 0.0016 323.8818 -0.2669 1.0000
family1 - family6 5 0.0039 0.0032 323.8818 1.2162 0.9526
family1 - family7 5 -0.0013 0.0018 323.8818 -0.7089 0.9987
family1 - family8 5 -0.0002 0.0020 323.8818 -0.1210 1.0000
family1 - family9 5 0.0011 0.0016 323.8818 0.6883 0.9989
family10 - family2 5 0.0009 0.0023 323.8818 0.3893 1.0000
family10 - family3 5 -0.0023 0.0020 323.8818 -1.1586 0.9643
family10 - family5 5 0.0001 0.0020 323.8818 0.0358 1.0000
family10 - family6 5 0.0044 0.0034 323.8818 1.2932 0.9329
family10 - family7 5 -0.0008 0.0022 323.8818 -0.3668 1.0000
family10 - family8 5 0.0003 0.0023 323.8818 0.1150 1.0000
family10 - family9 5 0.0016 0.0020 323.8818 0.8158 0.9964
family2 - family3 5 -0.0032 0.0020 323.8818 -1.6082 0.7999
family2 - family5 5 -0.0008 0.0020 323.8818 -0.4138 1.0000
family2 - family6 5 0.0035 0.0034 323.8818 1.0287 0.9829
family2 - family7 5 -0.0017 0.0022 323.8818 -0.7807 0.9973
family2 - family8 5 -0.0006 0.0023 323.8818 -0.2744 1.0000
family2 - family9 5 0.0007 0.0020 323.8818 0.3662 1.0000
family3 - family5 5 0.0024 0.0016 323.8818 1.4629 0.8718
family3 - family6 5 0.0067 0.0032 323.8818 2.0997 0.4754
family3 - family7 5 0.0015 0.0018 323.8818 0.8265 0.9960
family3 - family8 5 0.0026 0.0020 323.8818 1.2914 0.9334
family3 - family9 5 0.0039 0.0016 323.8818 2.4181 0.2776
family5 - family6 5 0.0043 0.0032 323.8818 1.3525 0.9143
family5 - family7 5 -0.0009 0.0018 323.8818 -0.4720 0.9999
family5 - family8 5 0.0002 0.0020 323.8818 0.0969 1.0000
family5 - family9 5 0.0016 0.0016 323.8818 0.9553 0.9894
family6 - family7 5 -0.0052 0.0033 323.8818 -1.5699 0.8205
family6 - family8 5 -0.0041 0.0034 323.8818 -1.2150 0.9528
family6 - family9 5 -0.0028 0.0032 323.8818 -0.8646 0.9946
family7 - family8 5 0.0011 0.0022 323.8818 0.4890 0.9999
family7 - family9 5 0.0024 0.0018 323.8818 1.3199 0.9249
family8 - family9 5 0.0014 0.0020 323.8818 0.6831 0.9990
family1 - family10 7 0.0027 0.0020 323.8818 1.3461 0.9165
family1 - family2 7 0.0012 0.0020 323.8818 0.5949 0.9996
family1 - family3 7 -0.0040 0.0016 323.8818 -2.4755 0.2479
family1 - family5 7 0.0006 0.0016 323.8818 0.3913 1.0000
family1 - family6 7 0.0054 0.0032 323.8818 1.7005 0.7459
family1 - family7 7 -0.0021 0.0018 323.8818 -1.1542 0.9651
family1 - family8 7 -0.0005 0.0020 323.8818 -0.2533 1.0000
family1 - family9 7 0.0009 0.0016 323.8818 0.5748 0.9997
family10 - family2 7 -0.0015 0.0023 323.8818 -0.6506 0.9993
family10 - family3 7 -0.0067 0.0020 323.8818 -3.3673 0.0237
family10 - family5 7 -0.0020 0.0020 323.8818 -1.0267 0.9831
family10 - family6 7 0.0027 0.0034 323.8818 0.8074 0.9966
family10 - family7 7 -0.0048 0.0022 323.8818 -2.2171 0.3971
family10 - family8 7 -0.0032 0.0023 323.8818 -1.3851 0.9029
family10 - family9 7 -0.0017 0.0020 323.8818 -0.8768 0.9940
family2 - family3 7 -0.0052 0.0020 323.8818 -2.6162 0.1839
family2 - family5 7 -0.0005 0.0020 323.8818 -0.2755 1.0000
family2 - family6 7 0.0042 0.0034 323.8818 1.2494 0.9447
family2 - family7 7 -0.0033 0.0022 323.8818 -1.5254 0.8430
family2 - family8 7 -0.0017 0.0023 323.8818 -0.7346 0.9983
family2 - family9 7 -0.0003 0.0020 323.8818 -0.1256 1.0000
family3 - family5 7 0.0047 0.0016 323.8818 2.8667 0.1008
family3 - family6 7 0.0094 0.0032 323.8818 2.9649 0.0778
family3 - family7 7 0.0019 0.0018 323.8818 1.0431 0.9813
family3 - family8 7 0.0035 0.0020 323.8818 1.7679 0.7033
family3 - family9 7 0.0050 0.0016 323.8818 3.0503 0.0615
family5 - family6 7 0.0048 0.0032 323.8818 1.5007 0.8547
family5 - family7 7 -0.0028 0.0018 323.8818 -1.5015 0.8544
family5 - family8 7 -0.0011 0.0020 323.8818 -0.5727 0.9997
family5 - family9 7 0.0003 0.0016 323.8818 0.1835 1.0000
family6 - family7 7 -0.0075 0.0033 323.8818 -2.2858 0.3538
family6 - family8 7 -0.0059 0.0034 323.8818 -1.7484 0.7159
family6 - family9 7 -0.0045 0.0032 323.8818 -1.4069 0.8947
family7 - family8 7 0.0016 0.0022 323.8818 0.7444 0.9981
family7 - family9 7 0.0031 0.0018 323.8818 1.6644 0.7677
family8 - family9 7 0.0014 0.0020 323.8818 0.7226 0.9985

11 Area Under the Curve (AUC)

AUC computed per individual via the trapezoid rule across all timepoints. metabolism_per_* is the primary metric matching the paper; corrected_fc and raw_fluorescence are retained for reference.

compute_auc <- function(df, value_col, group_vars) {
  df %>%
    filter(is.finite(time_hr), is.finite(.data[[value_col]])) %>%
    group_by(across(all_of(group_vars))) %>%
    summarise(
      AUC          = trapezoid_auc(time_hr, .data[[value_col]]),
      n_timepoints = n(),
      .groups      = "drop"
    ) %>%
    filter(is.finite(AUC))
}

# Only include grouping columns that are actually present in the data
individual_vars <- intersect(
  c("trace_id", "family_id_group", "treatment_group"),
  names(metabolism_df)
)

auc_metab_list <- list()
if (nrow(metabolism_df) > 0) {
  for (mc in active_meas_cols) {
    col <- paste0("metabolism_per_", mc)
    if (col %in% names(metabolism_df)) {
      auc_metab_list[[col]] <-
        compute_auc(metabolism_df, col, individual_vars) %>%
        mutate(metric = col)
    }
  }
}

auc_all <- bind_rows(auc_metab_list)

str(auc_all)
tibble [126 × 6] (S3: tbl_df/tbl/data.frame)
 $ trace_id       : chr [1:126] "1" "100" "101" "102" ...
 $ family_id_group: chr [1:126] "9" "7" "3" "1" ...
 $ treatment_group: chr [1:126] NA NA NA NA ...
 $ AUC            : num [1:126] 0.00868 0.05792 0.05097 0.0316 0.04869 ...
 $ n_timepoints   : int [1:126] 5 5 5 5 5 5 5 5 5 5 ...
 $ metric         : chr [1:126] "metabolism_per_area_mm2_measurement" "metabolism_per_area_mm2_measurement" "metabolism_per_area_mm2_measurement" "metabolism_per_area_mm2_measurement" ...

11.1 AUC summary tables

sum_vars <- intersect(
  c("metric", "family_id_group", "treatment_group"),
  names(auc_all)
)
auc_summary <- auc_all %>%
  group_by(across(all_of(sum_vars))) %>%
  summarise(
    n      = n(),
    mean   = mean(AUC, na.rm = TRUE),
    sd     = sd(AUC, na.rm = TRUE),
    se     = sd / sqrt(n),
    median = median(AUC, na.rm = TRUE),
    .groups = "drop"
  )

print(auc_summary)
# A tibble: 9 × 8
  metric     family_id_group treatment_group     n   mean      sd      se median
  <chr>      <chr>           <chr>           <int>  <dbl>   <dbl>   <dbl>  <dbl>
1 metabolis… 1               <NA>               20 0.0322 0.0243  0.00544 0.0340
2 metabolis… 10              <NA>               10 0.0330 0.0377  0.0119  0.0135
3 metabolis… 2               <NA>               10 0.0295 0.0204  0.00644 0.0281
4 metabolis… 3               <NA>               20 0.0445 0.0296  0.00663 0.0465
5 metabolis… 5               <NA>               20 0.0329 0.0214  0.00479 0.0319
6 metabolis… 6               <NA>                3 0.0152 0.00698 0.00403 0.0169
7 metabolis… 7               <NA>               13 0.0393 0.0300  0.00833 0.0453
8 metabolis… 8               <NA>               10 0.0341 0.0171  0.00540 0.0393
9 metabolis… 9               <NA>               20 0.0279 0.0297  0.00664 0.0105

12 Statistical Analysis

Each individual oyster (sample_id_group) is the observational unit. The model is built from whichever grouping factors are present: both family and treatment (with interaction) when both exist, or a one-way model when only one factor is available. Each plate maps to a unique family × treatment combination, so plate-level and group-level variance are confounded; interpret accordingly.

run_auc_stats <- function(auc_df) {
  empty <- tibble()

  has_family    <- "family_id_group" %in% names(auc_df) &&
    length(unique(na.omit(auc_df$family_id_group))) > 1
  has_treatment <- "treatment_group" %in% names(auc_df) &&
    length(unique(na.omit(auc_df$treatment_group))) > 1

  if (!has_family && !has_treatment) {
    return(list(model = NULL, anova = NULL,
                pairs_full = empty, pairs_family = empty,
                pairs_trt = empty,
                has_family = FALSE, has_treatment = FALSE))
  }

  if (has_family)    auc_df <- auc_df %>% mutate(family    = factor(family_id_group))
  if (has_treatment) auc_df <- auc_df %>% mutate(treatment = factor(treatment_group))

  formula_str <- if (has_family && has_treatment) {
    "AUC ~ family * treatment"
  } else if (has_family) {
    "AUC ~ family"
  } else {
    "AUC ~ treatment"
  }
  model     <- lm(as.formula(formula_str), data = auc_df)
  anova_res <- anova(model)

  if (has_family && has_treatment) {
    pairs_full   <- as.data.frame(pairs(emmeans(model, ~ family * treatment),
                                        adjust = "tukey"))
    pairs_family <- as.data.frame(pairs(emmeans(model, ~ family),
                                        adjust = "tukey"))
    pairs_trt    <- as.data.frame(pairs(emmeans(model, ~ treatment),
                                        adjust = "tukey"))
  } else if (has_family) {
    pairs_family <- as.data.frame(pairs(emmeans(model, ~ family),
                                        adjust = "tukey"))
    pairs_full   <- pairs_family
    pairs_trt    <- empty
  } else {
    pairs_trt    <- as.data.frame(pairs(emmeans(model, ~ treatment),
                                        adjust = "tukey"))
    pairs_full   <- pairs_trt
    pairs_family <- empty
  }

  list(
    model         = model,
    anova         = anova_res,
    pairs_full    = pairs_full,
    pairs_family  = pairs_family,
    pairs_trt     = pairs_trt,
    has_family    = has_family,
    has_treatment = has_treatment
  )
}

metrics_to_test <- unique(auc_all$metric)
stats_results   <- map(
  set_names(metrics_to_test),
  ~ run_auc_stats(auc_all %>% filter(metric == .x))
)

12.1 Results by metric

for (met in metrics_to_test) {
  stats <- stats_results[[met]]
  cat("\n\n### Metric:", met, "\n\n")
  cat("**ANOVA:**\n\n")
  cat(knitr::kable(as.data.frame(stats$anova), digits = 4, format = "pipe"), sep = "\n")
  cat("\n")
  if (stats$has_family && stats$has_treatment) {
    cat("**Pairwise: family × treatment (Tukey):**\n\n")
    cat(knitr::kable(as.data.frame(stats$pairs_full), digits = 4, format = "pipe"), sep = "\n")
  cat("\n")
    cat("**Pairwise: family main effect:**\n\n")
    cat(knitr::kable(as.data.frame(stats$pairs_family), digits = 4, format = "pipe"), sep = "\n")
  cat("\n")
    cat("**Pairwise: treatment main effect:**\n\n")
    cat(knitr::kable(as.data.frame(stats$pairs_trt), digits = 4, format = "pipe"), sep = "\n")
  cat("\n")
  } else if (stats$has_family) {
    cat("**Pairwise: family (Tukey):**\n\n")
    cat(knitr::kable(as.data.frame(stats$pairs_family), digits = 4, format = "pipe"), sep = "\n")
  cat("\n")
  } else if (stats$has_treatment) {
    cat("**Pairwise: treatment (Tukey):**\n\n")
    cat(knitr::kable(as.data.frame(stats$pairs_trt), digits = 4, format = "pipe"), sep = "\n")
  cat("\n")
  }
}

12.1.1 Metric: metabolism_per_area_mm2_measurement

ANOVA:

Df Sum Sq Mean Sq F value Pr(>F)
family 8 0.0047 6e-04 0.8172 0.5889
Residuals 117 0.0835 7e-04 NA NA

Pairwise: family (Tukey):

contrast estimate SE df t.ratio p.value
family1 - family10 -0.0009 0.0103 117 -0.0822 1.0000
family1 - family2 0.0027 0.0103 117 0.2584 1.0000
family1 - family3 -0.0123 0.0084 117 -1.4564 0.8729
family1 - family5 -0.0008 0.0084 117 -0.0902 1.0000
family1 - family6 0.0170 0.0165 117 1.0276 0.9824
family1 - family7 -0.0071 0.0095 117 -0.7474 0.9979
family1 - family8 -0.0020 0.0103 117 -0.1907 1.0000
family1 - family9 0.0043 0.0084 117 0.5054 0.9999
family10 - family2 0.0035 0.0119 117 0.2950 1.0000
family10 - family3 -0.0115 0.0103 117 -1.1069 0.9721
family10 - family5 0.0001 0.0103 117 0.0086 1.0000
family10 - family6 0.0178 0.0176 117 1.0149 0.9838
family10 - family7 -0.0063 0.0112 117 -0.5573 0.9998
family10 - family8 -0.0011 0.0119 117 -0.0940 1.0000
family10 - family9 0.0051 0.0103 117 0.4949 0.9999
family2 - family3 -0.0150 0.0103 117 -1.4476 0.8767
family2 - family5 -0.0034 0.0103 117 -0.3321 1.0000
family2 - family6 0.0143 0.0176 117 0.8145 0.9962
family2 - family7 -0.0098 0.0112 117 -0.8710 0.9941
family2 - family8 -0.0046 0.0119 117 -0.3890 1.0000
family2 - family9 0.0016 0.0103 117 0.1542 1.0000
family3 - family5 0.0115 0.0084 117 1.3662 0.9081
family3 - family6 0.0293 0.0165 117 1.7715 0.7006
family3 - family7 0.0052 0.0095 117 0.5454 0.9998
family3 - family8 0.0103 0.0103 117 0.9984 0.9854
family3 - family9 0.0166 0.0084 117 1.9618 0.5728
family5 - family6 0.0178 0.0165 117 1.0737 0.9769
family5 - family7 -0.0064 0.0095 117 -0.6673 0.9991
family5 - family8 -0.0012 0.0103 117 -0.1171 1.0000
family5 - family9 0.0050 0.0084 117 0.5956 0.9996
family6 - family7 -0.0241 0.0171 117 -1.4090 0.8923
family6 - family8 -0.0190 0.0176 117 -1.0787 0.9762
family6 - family9 -0.0127 0.0165 117 -0.7695 0.9975
family7 - family8 0.0051 0.0112 117 0.4574 0.9999
family7 - family9 0.0114 0.0095 117 1.1959 0.9559
family8 - family9 0.0062 0.0103 117 0.6034 0.9996

13 AUC Box Plots with Statistical Annotations

Significance labels: *** p < 0.001, ** p < 0.01, * p < 0.05. Brackets are drawn only for significant pairs (p < 0.05). Plots are generated for whichever grouping factors are present: treatment-only, family-only, all-groups, within-family, and within-treatment.

sig_label <- function(p) {
  case_when(p < 0.001 ~ "***", p < 0.01 ~ "**", p < 0.05 ~ "*",
            TRUE ~ "ns")
}

# Add significance brackets to an existing ggplot.
# pairs_df   : data frame with $contrast and $p.value columns
# group_levels: ordered character vector matching x-axis factor levels
# y_vals     : numeric vector of AUC values used to set bracket heights
add_sig_brackets <- function(p, pairs_df, group_levels, y_vals) {
  sig_pairs <- pairs_df %>%
    mutate(label = sig_label(p.value)) %>%
    filter(label != "ns")
  if (nrow(sig_pairs) == 0) return(p)

  y_max   <- max(y_vals, na.rm = TRUE)
  y_range <- diff(range(y_vals, na.rm = TRUE))
  step    <- y_range * 0.12

  for (i in seq_len(nrow(sig_pairs))) {
    parts <- str_split(as.character(sig_pairs$contrast[i]), " - ", 2)[[1]]
    g1 <- trimws(parts[1])
    g2 <- trimws(parts[2])
    x1 <- match(g1, group_levels)
    x2 <- match(g2, group_levels)
    if (is.na(x1) || is.na(x2)) next
    bar_y <- y_max + i * step
    p <- p +
      annotate("segment", x = x1, xend = x2,
               y = bar_y, yend = bar_y,
               colour = "black", linewidth = 0.6) +
      annotate("segment", x = x1, xend = x1,
               y = bar_y, yend = bar_y - step * 0.3,
               colour = "black", linewidth = 0.6) +
      annotate("segment", x = x2, xend = x2,
               y = bar_y, yend = bar_y - step * 0.3,
               colour = "black", linewidth = 0.6) +
      annotate("text", x = (x1 + x2) / 2,
               y = bar_y + step * 0.15,
               label = sig_pairs$label[i], size = 4.5)
  }
  p
}
for (met in metrics_to_test) {
  df      <- auc_all %>% filter(metric == met)
  stats   <- stats_results[[met]]
  y_lab   <- auc_y_label(met)
  has_fam <- stats$has_family
  has_trt <- stats$has_treatment

  # ── Treatment main effect (x = treatment, tick = treatment name) ───────
  if (has_trt) {
    df_p <- df %>%
      mutate(x = factor(treatment_group, levels = sort(unique(treatment_group))))
    grps <- levels(df_p$x)
    p <- ggplot(df_p, aes(x = x, y = AUC, fill = x)) +
      geom_boxplot(alpha = 0.6, outlier.shape = NA) +
      geom_jitter(width = 0.15, alpha = 0.4, size = 1.5) +
      scale_fill_manual(values = trt_colours[grps], guide = "none") +
      labs(x = "Treatment", y = y_lab) +
      theme_classic(base_size = 13)
    p <- add_sig_brackets(p, stats$pairs_trt, grps, df_p$AUC)
    print(p)
    ggsave(file.path(fig_dir, paste0("auc_treatment_", met, ".png")),
           p, width = 5, height = 5)
  }

  # ── Family main effect (x = family, tick = family name) ───────────────
  if (has_fam) {
    df_p <- df %>%
      mutate(x = factor(family_id_group, levels = sort(unique(family_id_group))))
    grps <- levels(df_p$x)
    p <- ggplot(df_p, aes(x = x, y = AUC, fill = x)) +
      geom_boxplot(alpha = 0.6, outlier.shape = NA) +
      geom_jitter(width = 0.15, alpha = 0.4, size = 1.5) +
      scale_fill_manual(values = fam_colours[grps], guide = "none") +
      labs(x = "Family", y = y_lab) +
      theme_classic(base_size = 13)
    p <- add_sig_brackets(p, stats$pairs_family, grps, df_p$AUC)
    print(p)
    ggsave(file.path(fig_dir, paste0("auc_family_", met, ".png")),
           p, width = 5, height = 5)
  }

  # Remaining plots require both factors
  if (!has_fam || !has_trt) next

  # ── All family:treatment groups (x = family:treatment) ─────────────────
  # emmeans contrasts use spaces; convert to colon to match tick labels
  pairs_fc <- stats$pairs_full %>%
    mutate(contrast = str_replace_all(
      contrast,
      "([a-z]+) ([a-z]+)( - )([a-z]+) ([a-z]+)",
      "\\1:\\2\\3\\4:\\5"
    ))
  df_p <- df %>%
    mutate(x = factor(
      paste(family_id_group, treatment_group, sep = ":"),
      levels = sort(unique(paste(family_id_group, treatment_group, sep = ":")))
    ))
  grps     <- levels(df_p$x)
  fill_map <- setNames(make_palette(length(grps)), grps)
  p <- ggplot(df_p, aes(x = x, y = AUC, fill = x)) +
    geom_boxplot(alpha = 0.6, outlier.shape = NA) +
    geom_jitter(width = 0.15, alpha = 0.4, size = 1.5) +
    scale_fill_manual(values = fill_map, guide = "none") +
    labs(x = "Family : Treatment", y = y_lab) +
    theme_classic(base_size = 13) +
    theme(axis.text.x = element_text(angle = 20, hjust = 1))
  p <- add_sig_brackets(p, pairs_fc, grps, df_p$AUC)
  print(p)
  ggsave(file.path(fig_dir, paste0("auc_all_groups_", met, ".png")),
         p, width = 6, height = 5)

  # ── Within each family: treatment comparison (x = family:treatment) ────
  # Tick labels are family:treatment so these plots are visually distinct
  # from the treatment main-effect plot above.
  for (fam in sort(unique(df$family_id_group))) {
    df_p <- df %>%
      filter(family_id_group == fam) %>%
      mutate(x = factor(
        paste(family_id_group, treatment_group, sep = ":"),
        levels = sort(unique(paste(family_id_group, treatment_group, sep = ":")))
      ))
    grps     <- levels(df_p$x)
    pairs_sub <- pairs_fc %>%
      filter(str_count(contrast, paste0(fam, ":")) == 2)
    p <- ggplot(df_p, aes(x = x, y = AUC, fill = x)) +
      geom_boxplot(alpha = 0.6, outlier.shape = NA) +
      geom_jitter(width = 0.15, alpha = 0.4, size = 1.5) +
      scale_fill_manual(values = fill_map[grps], guide = "none") +
      labs(x = "Family : Treatment", y = y_lab) +
      theme_classic(base_size = 13)
    p <- add_sig_brackets(p, pairs_sub, grps, df_p$AUC)
    print(p)
    ggsave(file.path(fig_dir, paste0("auc_", fam, "_trt_", met, ".png")),
           p, width = 5, height = 5)
  }

  # ── Within each treatment: family comparison (x = family:treatment) ────
  # Tick labels are family:treatment so these plots are visually distinct
  # from the family main-effect plot above.
  for (trt in sort(unique(df$treatment_group))) {
    df_p <- df %>%
      filter(treatment_group == trt) %>%
      mutate(x = factor(
        paste(family_id_group, treatment_group, sep = ":"),
        levels = sort(unique(paste(family_id_group, treatment_group, sep = ":")))
      ))
    grps      <- levels(df_p$x)
    pairs_sub <- pairs_fc %>%
      filter(str_count(contrast, paste0(":", trt)) == 2)
    p <- ggplot(df_p, aes(x = x, y = AUC, fill = x)) +
      geom_boxplot(alpha = 0.6, outlier.shape = NA) +
      geom_jitter(width = 0.15, alpha = 0.4, size = 1.5) +
      scale_fill_manual(values = fill_map[grps], guide = "none") +
      labs(x = "Family : Treatment", y = y_lab) +
      theme_classic(base_size = 13)
    p <- add_sig_brackets(p, pairs_sub, grps, df_p$AUC)
    print(p)
    ggsave(file.path(fig_dir, paste0("auc_", trt, "_fam_", met, ".png")),
           p, width = 5, height = 5)
  }
}

14 Save Output Data

write_csv(auc_all,      file.path(out_dir, "auc_all_metrics.csv"))
write_csv(auc_summary,  file.path(out_dir, "auc_summary.csv"))

if (nrow(metabolism_df) > 0)
  write_csv(metabolism_df,
            file.path(out_dir, "metabolism.csv"))

stats_compiled <- map_dfr(metrics_to_test, function(met) {
  bind_rows(
    stats_results[[met]]$pairs_full %>%
      mutate(comparison = "family:treatment"),
    stats_results[[met]]$pairs_family %>%
      mutate(comparison = "family"),
    stats_results[[met]]$pairs_trt %>%
      mutate(comparison = "treatment")
  ) %>% mutate(metric = met)
})

write_csv(stats_compiled,
          file.path(out_dir, "pairwise_stats.csv"))

message("Output files written to: ", out_dir)