The Zanco et al. (2025) trade-off matrix: unpaired comparisons and clustering

Overview

This vignette replicates the second worked example from the Nutrigonometry validation paper: the pairwise angular trade-off matrix across egg-to-adult viability and wing size in Drosophila melanogaster, measured at three temperature regimes and in two locally-adapted populations (Zanco et al. 2025).

Unlike the Lee example (see vignette("getting-started")), viability is measured per vial and wing size per individual fly – the two traits do not share a sampling unit, so every comparison here is unpaired (paired = FALSE in nutrigon_compare()).

With 2 traits x 3 temperatures x 2 populations = 12 fitted surfaces, there are $\binom{12}{2} = 66$ pairwise comparisons. At B = 200 bootstrap resamples each, building the full matrix from scratch takes several minutes, so this package ships the result pre-built as zanco_tradeoff_matrix. This vignette loads it directly; the code used to build it (fitting the 12 surfaces, then looping over all pairs) is shown alongside so you can see exactly how it was constructed and adapt it to your own data.

library(nutrigon)
data(zanco_survival)
data(zanco_wingsize)
data(zanco_tradeoff_matrix)

TEMPS <- c("18fluc", "25fluc", "28fluc")
POPS  <- c("Melbourne", "Townsville")

cat(sprintf("Survival data: %d rows | Wing size data: %d rows\n",
            nrow(zanco_survival), nrow(zanco_wingsize)))
#> Survival data: 1495 rows | Wing size data: 2073 rows
cat(sprintf("Environments: %d temperatures x %d populations = %d surfaces per trait\n",
            length(TEMPS), length(POPS), length(TEMPS) * length(POPS)))
#> Environments: 3 temperatures x 2 populations = 6 surfaces per trait
cat(sprintf("Total surfaces: %d | Pairwise comparisons: %d\n",
            2 * length(TEMPS) * length(POPS),
            choose(2 * length(TEMPS) * length(POPS), 2)))
#> Total surfaces: 12 | Pairwise comparisons: 66

How the matrix was built

The 12 surfaces were fitted with nutrigon_fit(), one per trait x temperature x population combination:

surfaces <- list()
for (temp in TEMPS) {
  for (pop in POPS) {
    sub_s <- zanco_survival[zanco_survival$Temperature == temp &
                             zanco_survival$Population  == pop, ]
    key_s <- paste("viability", temp, pop, sep = "_")
    surfaces[[key_s]] <- nutrigon_fit(
      sub_s, "Prot", "Carb", "Prop_viable",
      design_id_col = "design_id", peak = TRUE
    )

    sub_w <- zanco_wingsize[zanco_wingsize$Temperature == temp &
                             zanco_wingsize$Population  == pop, ]
    key_w <- paste("wingsize", temp, pop, sep = "_")
    surfaces[[key_w]] <- nutrigon_fit(
      sub_w, "Prot", "Carb", "WingSize",
      design_id_col = "design_id", peak = TRUE
    )
  }
}

zanco_tradeoff_matrix$surfaces already contains these 12 fitted objects, so we use that directly:

surfaces <- zanco_tradeoff_matrix$surfaces
for (nm in names(surfaces)) {
  cat(sprintf("  [%s] theta = %.2f deg | elongation = %.2f | fallback = %s\n",
              nm,
              surfaces[[nm]]$optimum$theta,
              surfaces[[nm]]$elongation,
              surfaces[[nm]]$optimum$used_fallback))
}
#>   [viability_18fluc_Melbourne] theta = 58.60 deg | elongation = 62.36 | fallback = TRUE
#>   [wingsize_18fluc_Melbourne] theta = 36.98 deg | elongation = 1.41 | fallback = TRUE
#>   [viability_18fluc_Townsville] theta = 59.19 deg | elongation = 7.75 | fallback = TRUE
#>   [wingsize_18fluc_Townsville] theta = 37.03 deg | elongation = 5.91 | fallback = TRUE
#>   [viability_25fluc_Melbourne] theta = 57.19 deg | elongation = 6.99 | fallback = TRUE
#>   [wingsize_25fluc_Melbourne] theta = 54.24 deg | elongation = 1.81 | fallback = FALSE
#>   [viability_25fluc_Townsville] theta = 60.27 deg | elongation = 46.91 | fallback = TRUE
#>   [wingsize_25fluc_Townsville] theta = 82.86 deg | elongation = 2.05 | fallback = TRUE
#>   [viability_28fluc_Melbourne] theta = 82.72 deg | elongation = 109.72 | fallback = TRUE
#>   [wingsize_28fluc_Melbourne] theta = 82.22 deg | elongation = 26.12 | fallback = TRUE
#>   [viability_28fluc_Townsville] theta = 61.66 deg | elongation = 4.69 | fallback = TRUE
#>   [wingsize_28fluc_Townsville] theta = 76.34 deg | elongation = 9.16 | fallback = FALSE

Each pairwise comparison was then run with nutrigon_compare(..., paired = FALSE, resample = "stratified", B = 200) (resampling within design points, since each diet has several replicate vials or flies):

res <- nutrigon_compare(
  raw_i, "Prot", "Carb", z_i,
  raw_j, "Prot", "Carb", z_j,
  design_id_col1 = "design_id", design_id_col2 = "design_id",
  paired   = FALSE,
  resample = "stratified",
  B        = 200,
  seed     = pair_id
)

The result of running this over all 66 pairs is stored in zanco_tradeoff_matrix:

Theta_matrix <- zanco_tradeoff_matrix$Theta
round(Theta_matrix[1:4, 1:4], 2)
#>                             viability_18fluc_Melbourne
#> viability_18fluc_Melbourne                        0.00
#> wingsize_18fluc_Melbourne                        21.85
#> viability_18fluc_Townsville                       0.15
#> wingsize_18fluc_Townsville                       21.52
#>                             wingsize_18fluc_Melbourne
#> viability_18fluc_Melbourne                      21.85
#> wingsize_18fluc_Melbourne                        0.00
#> viability_18fluc_Townsville                     21.46
#> wingsize_18fluc_Townsville                       0.06
#>                             viability_18fluc_Townsville
#> viability_18fluc_Melbourne                         0.15
#> wingsize_18fluc_Melbourne                         21.46
#> viability_18fluc_Townsville                        0.00
#> wingsize_18fluc_Townsville                        21.85
#>                             wingsize_18fluc_Townsville
#> viability_18fluc_Melbourne                       21.52
#> wingsize_18fluc_Melbourne                         0.06
#> viability_18fluc_Townsville                      21.85
#> wingsize_18fluc_Townsville                        0.00

Headline result: viability-wing size trade-off by temperature and population

pairwise_df <- zanco_tradeoff_matrix$pairwise

headline <- pairwise_df[
  pairwise_df$trait_1 != pairwise_df$trait_2 &
  pairwise_df$temperature_1 == pairwise_df$temperature_2 &
  pairwise_df$population_1  == pairwise_df$population_2, ]
headline$temp_label <- sub("fluc", "\u00b15\u00b0C", headline$temperature_1)

print(headline[order(headline$temperature_1, headline$population_1),
               c("temp_label", "population_1", "theta_deg", "ci_lower", "ci_upper")],
      row.names = FALSE)
#>  temp_label population_1 theta_deg ci_lower ci_upper
#>      18±5°C    Melbourne     21.85    16.67    26.80
#>      18±5°C   Townsville     21.85    17.34    27.95
#>      25±5°C    Melbourne      4.14     0.25    20.24
#>      25±5°C   Townsville     21.67     0.03    26.32
#>      28±5°C    Melbourne      1.06     0.07    28.09
#>      28±5°C   Townsville     15.87     0.77    24.91

The viability-wing size trade-off collapses towards zero for Melbourne at warmer temperatures, while Townsville’s angular separation stays fairly stable across temperatures – consistent with population-specific thermal adaptation of the nutritional optima.

Hierarchical clustering and multidimensional scaling

Treating Theta as a distance matrix lets us cluster the 12 surfaces and visualise them in two dimensions:

dist_obj <- as.dist(Theta_matrix)
hc       <- hclust(dist_obj, method = "average")
mds      <- cmdscale(dist_obj, k = 2, eig = TRUE)

cat(sprintf("MDS: %.1f%% of variance explained by first 2 axes.\n",
            100 * sum(mds$eig[1:2]) / sum(abs(mds$eig))))
#> MDS: 92.8% of variance explained by first 2 axes.

library(ggplot2)

if (requireNamespace("ggdendro", quietly = TRUE)) {
  ddata <- ggdendro::dendro_data(hc)
  ggplot() +
    geom_segment(data = ggdendro::segment(ddata),
                 aes(x = x, y = y, xend = xend, yend = yend)) +
    geom_text(data = ggdendro::label(ddata),
              aes(x = x, y = y, label = label),
              hjust = 1, size = 3, nudge_y = -1) +
    coord_flip() +
    scale_y_reverse(expand = expansion(mult = c(0.3, 0.05))) +
    labs(x = NULL, y = "Angular distance (degrees)",
         title = "Hierarchical clustering of the 12 fitted surfaces") +
    theme_minimal(base_size = 11) +
    theme(axis.text.y = element_blank(), panel.grid = element_blank())
} else {
  plot(hc, main = "Hierarchical clustering of the 12 fitted surfaces",
       xlab = "", sub = "")
}

library(ggplot2)
mds_df <- data.frame(
  surface = rownames(Theta_matrix),
  mds1    = mds$points[, 1],
  mds2    = mds$points[, 2]
)
mds_df$trait <- ifelse(startsWith(mds_df$surface, "viability"),
                        "Viability", "Wing size")

ggplot(mds_df, aes(x = mds1, y = mds2, colour = trait)) +
  geom_point(size = 3) +
  ggrepel::geom_text_repel(aes(label = surface), size = 2.8, show.legend = FALSE) +
  scale_colour_manual(values = c(Viability = "#E69F00", `Wing size` = "#009E73"),
                      name = "Trait") +
  labs(x = "MDS axis 1", y = "MDS axis 2",
       title = "Multidimensional scaling of angular separations") +
  theme_bw(base_size = 11) +
  theme(panel.grid = element_blank())

Surface overlay plots

The plot helpers nutrigon_grid(), nutrigon_hull(), and nutrigon_region_polygon() let us overlay both traits’ peak regions on the same design hull, faceted by temperature and population.

build_overlay_data <- function(surfaces, temps, pops) {
  region_rows <- list(); hull_rows <- list()
  for (temp in temps) {
    for (pop in pops) {
      v_key <- paste("viability", temp, pop, sep = "_")
      w_key <- paste("wingsize",  temp, pop, sep = "_")
      temp_label <- sub("fluc", "\u00b15\u00b0C", temp)

      for (kk in c(v_key, w_key)) {
        trait_label <- if (startsWith(kk, "viability")) "Viability" else "Wing size"
        s <- surfaces[[kk]]
        reg <- nutrigon_region_polygon(s)
        names(reg)[1:2] <- c("Prot", "Carb")
        region_rows[[length(region_rows) + 1L]] <- cbind(
          reg, trait = trait_label, temperature = temp_label, population = pop
        )
      }
      hull_df <- nutrigon_hull(surfaces[[v_key]])
      names(hull_df)[1:2] <- c("Prot", "Carb")
      hull_rows[[length(hull_rows) + 1L]] <- cbind(
        hull_df, temperature = temp_label, population = pop
      )
    }
  }
  list(regions = do.call(rbind, region_rows),
       hulls   = do.call(rbind, hull_rows))
}

ov <- build_overlay_data(surfaces, TEMPS, POPS)

temp_order <- paste0(c("18", "25", "28"), "\u00b15\u00b0C")
ov$regions$temperature <- factor(ov$regions$temperature, levels = temp_order)
ov$hulls$temperature    <- factor(ov$hulls$temperature,  levels = temp_order)
ov$regions$population   <- factor(ov$regions$population, levels = POPS)
ov$hulls$population     <- factor(ov$hulls$population,   levels = POPS)

library(ggplot2)
TRAIT_COLS <- c("Viability" = "#E69F00", "Wing size" = "#009E73")

ggplot() +
  geom_polygon(data = ov$hulls,
               aes(x = Prot, y = Carb),
               fill = NA, colour = "black", linewidth = 0.5) +
  geom_polygon(data = ov$regions,
               aes(x = Prot, y = Carb,
                   fill   = trait,
                   colour = trait,
                   group  = interaction(trait, temperature, population)),
               alpha = 0.35, linewidth = 0.6, inherit.aes = FALSE) +
  scale_fill_manual(values = TRAIT_COLS,   name = "Trait") +
  scale_colour_manual(values = TRAIT_COLS, name = "Trait") +
  facet_grid(population ~ temperature) +
  labs(x = "Protein (g/L)", y = "Carbohydrate (g/L)",
       title = "Peak regions: viability vs wing size") +
  theme_bw(base_size = 11) +
  theme(panel.grid      = element_blank(),
        strip.text      = element_text(face = "bold", size = 10),
        legend.position = "right")

Refitting viability with a quasibinomial GLM

The matrix above fits every surface, including viability, with ordinary least squares on the raw proportion Prop_viable – this matches the original published analysis and is what zanco_tradeoff_matrix was built with. But viability is naturally a bounded proportion derived from binomial counts (Total_Adults out of Total_Eggs per vial), with visibly more spread than a binomial count alone would predict – exactly the situation a quasibinomial GLM is built for: it models the proportion on the logit scale while inflating standard errors to account for overdispersion, without requiring individual-egg-level data.

nutrigon_fit() and nutrigon_compare() support this directly via model_type = "glm" and family = quasibinomial(). As one example, refit the Melbourne 180b1C viability surface both ways and compare:

mel18 <- zanco_survival[zanco_survival$Temperature == "18fluc" &
                         zanco_survival$Population  == "Melbourne", ]

fit_ols  <- zanco_tradeoff_matrix$surfaces[["viability_18fluc_Melbourne"]]
fit_quasi <- nutrigon_fit(mel18, "Prot", "Carb", "Prop_viable",
                           design_id_col = "design_id",
                           model_type = "glm", family = quasibinomial())

data.frame(
  fit         = c("lm (shipped matrix)", "glm, quasibinomial"),
  theta_deg   = c(fit_ols$optimum$theta, fit_quasi$optimum$theta),
  elongation  = c(fit_ols$elongation,    fit_quasi$elongation)
)
#>                   fit theta_deg elongation
#> 1 lm (shipped matrix)  58.59938   62.36493
#> 2  glm, quasibinomial  57.82001  299.35175

The optimum’s location can shift somewhat between the two fits: the OLS surface is quadratic in the raw proportion, while the quasibinomial surface is quadratic on the logit scale (see ?nutrigon_fit, Details) – these are genuinely different statistical models of the same data, not two numerical routes to the same answer. nutrigon_curvature() and nutrigon_grid() both work identically on the quasibinomial fit; for nutrigon_grid(), remember that predicted is then on the logit scale (see ?nutrigon_grid, GLM fits section).

A full quasibinomial refit of all 6 viability surfaces (paired with the existing lm-fitted wing-size surfaces via nutrigon_compare()’s model_type1/model_type2) would give a more principled version of the matrix above, at the cost of re-running the ~66-comparison bootstrap. That full refit is not included here to keep the vignette fast to build; the single-surface comparison above illustrates the mechanics.

Joint uncertainty across the matrix

Pairwise intervals from separate nutrigon_compare() calls are marginal: they do not carry the dependence among entries that share a surface. nutrigon_matrix() bootstraps each surface once and uses the same replicate index in every pair, so each replicate is a complete matrix and any summary (e.g. a dendrogram) can be given bootstrap uncertainty. It also returns the ordering of the optimum angles; with two nutrients the matrix is one-dimensional, so this ordering carries all of its information. The analysis in the paper (survival as a quasi-binomial GLM) corresponds to:

specs <- list()
for (temp in TEMPS) for (pop in POPS) {
  s <- zanco_survival[zanco_survival$Temperature == temp &
                        zanco_survival$Population == pop, ]
  w <- zanco_wingsize[zanco_wingsize$Temperature == temp &
                        zanco_wingsize$Population == pop, ]
  specs[[paste("survival", temp, pop, sep = "_")]] <- list(
    data = s, x_col = "Prot", y_col = "Carb", z_col = "Prop_viable",
    model_type = "glm", family = quasibinomial(), design_id_col = "design_id")
  specs[[paste("wingsize", temp, pop, sep = "_")]] <- list(
    data = w, x_col = "Prot", y_col = "Carb", z_col = "WingSize",
    design_id_col = "design_id")
}
m <- nutrigon_matrix(specs, B = 200, resample = "stratified", seed = 2025)
m$clustering$stability

Headline strip chart

library(ggplot2)
headline$pop_f  <- factor(headline$population_1, levels = POPS)
headline$temp_f <- factor(headline$temp_label,   levels = temp_order)

ggplot(headline,
       aes(x = temp_f, y = theta_deg,
           colour = pop_f, fill = pop_f, shape = pop_f,
           ymin = ci_lower, ymax = ci_upper)) +
  geom_hline(yintercept = 0, linetype = "dashed", colour = "grey60") +
  geom_linerange(position = position_dodge(0.35), linewidth = 0.8) +
  geom_point(position = position_dodge(0.35), size = 3.5, stroke = 0.8) +
  scale_colour_manual(values = c(Melbourne = "#4477AA", Townsville = "#AA3377"),
                      name = "Population") +
  scale_fill_manual(values   = c(Melbourne = "#4477AA", Townsville = "#AA3377"),
                    name = "Population") +
  scale_shape_manual(values  = c(Melbourne = 16, Townsville = 17),
                     name = "Population") +
  labs(x = "Temperature regime", y = "Angular separation (\u00b0)",
       title = "Viability-wing size trade-off across environments") +
  theme_bw(base_size = 12) +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor   = element_blank())