---
title: "Cross-case-study summary: three taxa, one panel"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Cross-case-study summary: three taxa, one panel}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7.5,
  fig.height = 6,
  dpi = 96
)
```

This vignette mirrors Manuscript Figure 5 of the `bayesTLS` supplement -- a single
multi-taxon panel of the thermal-sensitivity parameter `z` and the critical
temperature `CTmax` across the three redistributable case studies -- but draws the `freqTLS`
**Confidence Eye** for each estimate instead of a Bayesian posterior ridge. Each
row is a likelihood **confidence interval** (a pale lens with a hollow point
estimate), never a posterior density.

**This vignette builds without Stan.** Its `freqTLS` maximum-likelihood fits and
intervals are read from a version-stamped maintainer cache generated by the
package's documented fitting code. The individual case-study articles and test
suite run the same fit/profile paths live. The Bayesian multi-taxon ridge plot
is the `bayesTLS` complement; see the closing note.

```{r setup}
library(freqTLS)
```

## The three taxa and the cross-study question

The three shared case-study datasets span three corners of thermal physiology:

* **brown shrimp** *Crangon crangon* -- a crustacean, lethal mortality, ungrouped;
* **zebrafish** *Danio rerio* -- a fish, lethal mortality, across three life
  stages (young embryos, old embryos, larvae); this panel uses the
  `zebrafish_lethal` life-stage dataset — a *separate* experiment from the
  oxygen-gradient `zebrafish_o2` data in `vignette("case-study-zebrafish")`;
* **vinegar fly** *Drosophila suzukii* -- an insect, lethal mortality, across the
  two sexes.

The cross-study question is descriptive: laid side by side, how do thermal limits
compare across a crustacean, a fish (resolved by developmental stage), and an
insect (resolved by sex)? `CTmax` (the critical
temperature at a fixed reference exposure) places each taxon on the temperature
axis; `z` (the change in temperature, in degrees Celsius, that multiplies tolerated
exposure time tenfold) measures how sharply tolerated exposure responds to
temperature: a *smaller* `z` is a steeper response (a small warming sharply cuts
tolerated time), a *larger* `z` a more gradual one.

**One caveat governs the whole panel.** The reference exposures differ by study,
because each follows the convention of its source assay:

| Taxon | Endpoint | Reference exposure (`tref`) | Threshold |
| --- | --- | --- | --- |
| Shrimp | lethal | 1 hour | relative midpoint |
| Zebrafish (per stage) | lethal | 1 hour | relative midpoint |
| *D. suzukii* (per sex) | lethal | 4 hours (240 min) | relative midpoint |

Because a `CTmax` is defined *at* its reference exposure, the four `CTmax` values
are **not** on a single common time scale: the fly number is the fitted relative
midpoint temperature at 4 hours, whereas the shrimp and zebrafish fits use a
1-hour reference. The panel is therefore an **illustrative
cross-taxon synthesis**, not a single common-scale comparison. `z`, by contrast,
is a slope -- degrees per tenfold change in time -- and is comparable across taxa
regardless of the reference exposure. Read the `CTmax` facet as four study-specific
anchors and the `z` facet as a like-with-like comparison of duration sensitivity.

## The three fits (version-stamped cache)

Each fit uses the configuration locked for its case study. The grouped fits
(zebrafish, fly) estimate a separate `CTmax` and `z` per level with a shared shape
(`low`, `up`, `k`); the ungrouped shrimp fit estimates one of each. All three use
**beta-binomial** survival counts. The cache records the package/source version,
R and TMB versions, input checksums, and the exact family, grouping, reference
exposure, threshold, shape, and interval configuration. The build script is
`data-raw/build_case_study_summary_cache.R`; the per-study vignettes and tests
exercise the live paths and their data-adequacy warnings.

The vinegar-fly data ships per-individual (`dsuzukii`); we aggregate it to
`(temp, time, sex)` counts in base R, exactly as the per-study vignette does.

```{r fits}
summary_cache_path <- system.file(
  "extdata", "case_study_summary_cache.rds", package = "freqTLS"
)
if (!nzchar(summary_cache_path) || !file.exists(summary_cache_path)) {
  stop("The shipped cross-case-study cache is missing; reinstall freqTLS from a complete source tarball.")
}
summary_cache <- readRDS(summary_cache_path)
summary_cache$meta[c("schema_version", "generated_on", "freqTLS_version",
                     "freqTLS_source_commit", "R_version", "TMB_version")]
```

## Assembling the cross-study table

The panel is built from one combined data frame of
`(label, parameter, estimate, conf.low, conf.high)`, obtained by calling
`confint(..., method = "profile")` on each fit. For the grouped fits we ask for the
per-level parameter names (`CTmax:young_embryos`, `z:M`, and so on); for the
ungrouped fits we ask for the bare `CTmax` and `z`. Every interval here is a
profile-likelihood confidence interval, and -- as the `conf.status` column
confirms -- every profile closes.

```{r combine}
panel <- summary_cache$panel
# Six rows, two parameters: a tidy printout of the headline numbers.
panel_wide <- data.frame(
  Group           = panel$label[panel$parameter == "CTmax"],
  `CTmax estimate` = round(panel$estimate[panel$parameter == "CTmax"], 2),
  `CTmax 95% CI`  = sprintf("[%.2f, %.2f]",
                            panel$conf.low[panel$parameter == "CTmax"],
                            panel$conf.high[panel$parameter == "CTmax"]),
  `z estimate`     = round(panel$estimate[panel$parameter == "z"], 2),
  `z 95% CI`       = sprintf("[%.2f, %.2f]",
                            panel$conf.low[panel$parameter == "z"],
                            panel$conf.high[panel$parameter == "z"]),
  check.names = FALSE
)
knitr::kable(
  panel_wide,
  caption = "Six taxon/group rows: CTmax (at the study reference exposure) and z, each with its profile-likelihood 95% confidence interval."
)
```

## Cross-taxon validation against bayesTLS and the two-stage estimator

The same headline numbers are shown beside the classical two-stage estimator
and the `bayesTLS` posterior, read from the maintainer-built benchmark cache.
The `freqTLS` and `bayesTLS` fits use the matched relative-midpoint,
constant-shape configuration and the same per-study reference exposure. The
classical estimator uses absolute LT50 and is therefore an approximate
comparator for these lethal datasets, whose asymptotes lie near zero and one.

```{r three-way-summary, echo = FALSE}
cache_path <- system.file("extdata", "bayesTLS_benchmark_cache.rds", package = "freqTLS")
fmt3 <- function(v) ifelse(is.na(v[1]), "--", sprintf("%.2f [%.2f, %.2f]", v[1], v[2], v[3]))

# Map each summary row label to its benchmark-cache dataset key.
keys <- c(
  "Shrimp"                   = "shrimp",
  "Zebrafish: young embryos" = "zebrafish:young_embryos",
  "Zebrafish: old embryos"   = "zebrafish:old_embryos",
  "Zebrafish: larvae"        = "zebrafish:larvae",
  "D. suzukii: female"       = "dsuzukii:F",
  "D. suzukii: male"         = "dsuzukii:M"
)

if (nzchar(cache_path) && file.exists(cache_path)) {
  cache <- readRDS(cache_path)
  cache_pick <- function(df, key, parm, cols) {
    if (is.null(df)) return(c(NA_real_, NA_real_, NA_real_))
    r <- df[df$dataset == key & df$parameter == parm, , drop = FALSE]
    if (nrow(r) == 0L) return(c(NA_real_, NA_real_, NA_real_))
    as.numeric(r[1L, cols])
  }
  pro_pick <- function(label, parm) {
    r <- panel[panel$label == label & panel$parameter == parm, , drop = FALSE]
    if (nrow(r) == 0L) return(c(NA_real_, NA_real_, NA_real_))
    c(r$estimate[1L], r$conf.low[1L], r$conf.high[1L])
  }
  quantities <- list(c("CTmax", "CTmax (°C)"), c("z", "z (°C / decade)"))
  rows <- do.call(rbind, lapply(names(keys), function(label) {
    key <- keys[[label]]
    do.call(rbind, lapply(quantities, function(q) {
      ts <- cache_pick(cache$two_stage, key, q[1], c("estimate", "lower", "upper"))
      by <- cache_pick(cache$bayesian,  key, q[1], c("median", "lower", "upper"))
      pr <- pro_pick(label, q[1])
      data.frame(
        Group = label, Quantity = q[2],
        `Two-stage (delta CI)`    = fmt3(ts),
        `bayesTLS (95% CrI)`      = fmt3(by),
        `freqTLS (profile CI)` = fmt3(pr),
        check.names = FALSE
      )
    }))
  }))
  knitr::kable(rows, row.names = FALSE,
               caption = "Cross-taxon three-way comparison: CTmax and z per taxon/group from the classical two-stage estimator, bayesTLS (posterior median + 95% CrI), and freqTLS (profile-likelihood estimate + 95% CI). The two model fits share each study's relative-threshold, constant-shape configuration; the two-stage estimate is an absolute-LT50 approximation for these near-0/near-1 lethal curves.")
} else {
  stop("The shipped bayesTLS benchmark cache is missing; reinstall freqTLS from a complete source tarball.")
}
```

## The panel: a Confidence Eye per taxon and group

Each row is drawn as an honest Confidence Eye -- a pale, shallow horizontal lens
whose width is exactly the 95% confidence interval and whose cosine taper is
tallest at the estimate, with a hollow point marking the estimate itself. The two
facets carry independent x-axes (`CTmax` in degrees Celsius; `z` in degrees Celsius
per tenfold change in exposure time), because the two parameters live on different
scales and, for `CTmax`, on different reference exposures.

The geometry follows the `freqTLS` Confidence-Eye contract: the shallow,
wide lens reads as a confidence *interval*, not a probability density, and a
profile that did not close would render as a hollow point with **no** lens. All
12 headline profiles here close (six groups for each of two parameters), so
every row carries a lens.

```{r panel, fig.width = 9, fig.height = 6, fig.alt = "Two-facet Confidence-Eye panel. Left facet, CTmax in degrees Celsius: six pale horizontal confidence-interval lenses with hollow point estimates, ordered top to bottom as shrimp near 31.8, the three zebrafish stages near 39.8 to 41.4, and the two D. suzukii sexes near 35.2. Right facet, z in degrees Celsius per tenfold time change: matching lenses with shrimp near 2.2, zebrafish stages near 1.8 to 2.0, and the two fly sexes near 3.0 to 3.2. Each CTmax row is annotated with its reference exposure (1 hour or 4 hours), underscoring that the CTmax values are not on a common time scale."}
# Build the honest Confidence-Eye geometry by hand so all three fits share one
# cross-taxon panel: a cosine-tapered pale lens (width = confidence interval)
# plus a hollow point, faceted by parameter with free x-axes. This reuses the
# package lens shape (see ?plot_confidence_eye) but spans every fit in one figure.
stopifnot(requireNamespace("ggplot2", quietly = TRUE))

# Fixed top-to-bottom row order (shrimp, zebrafish stages, fly sexes).
row_order <- c(
  "Shrimp",
  "Zebrafish: young embryos", "Zebrafish: old embryos", "Zebrafish: larvae",
  "D. suzukii: female", "D. suzukii: male"
)
panel$row <- match(panel$label, row_order)
panel$parameter <- factor(panel$parameter, levels = c("CTmax", "z"),
                          labels = c("CTmax (°C)", "z (°C / decade)"))

# One cosine-tapered lens polygon per row (tallest at the estimate, zero at each
# bound), built per facet so the free x-axes do not distort the taper.
lens_df <- do.call(rbind, lapply(seq_len(nrow(panel)), function(i) {
  x <- seq(panel$conf.low[i], panel$conf.high[i], length.out = 80)
  d_lo <- max(panel$estimate[i] - panel$conf.low[i], .Machine$double.eps)
  d_hi <- max(panel$conf.high[i] - panel$estimate[i], .Machine$double.eps)
  frac <- ifelse(x <= panel$estimate[i], (panel$estimate[i] - x) / d_lo,
                 (x - panel$estimate[i]) / d_hi)
  w <- 0.32 * cos((pi / 2) * pmin(pmax(frac, 0), 1))
  data.frame(id = i, parameter = panel$parameter[i], x = x,
             ymin = panel$row[i] - w, ymax = panel$row[i] + w)
}))

# Reference-exposure annotation for the CTmax facet only.
tref_lab <- data.frame(
  parameter = factor("CTmax (°C)",
                     levels = levels(panel$parameter)),
  row = panel$row[panel$parameter == "CTmax (°C)"],
  x   = panel$conf.high[panel$parameter == "CTmax (°C)"],
  lab = c("1 h", "1 h", "1 h", "1 h", "4 h", "4 h")
)

ggplot2::ggplot() +
  ggplot2::geom_ribbon(
    data = lens_df,
    ggplot2::aes(x = x, ymin = ymin, ymax = ymax, group = id),
    fill = "#1b7837", colour = NA, alpha = 0.30
  ) +
  ggplot2::geom_point(
    data = panel,
    ggplot2::aes(x = estimate, y = row),
    shape = 21, fill = "white", colour = "#1b7837", size = 3, stroke = 1
  ) +
  ggplot2::geom_text(
    data = tref_lab,
    ggplot2::aes(x = x, y = row, label = lab),
    hjust = -0.25, size = 2.7, colour = "grey35"
  ) +
  ggplot2::scale_y_continuous(
    breaks = seq_along(row_order), labels = row_order,
    trans = "reverse", expand = ggplot2::expansion(add = 0.7)
  ) +
  ggplot2::scale_x_continuous(
    expand = ggplot2::expansion(mult = c(0.05, 0.20))
  ) +
  ggplot2::facet_wrap(~ parameter, scales = "free_x") +
  ggplot2::labs(
    x = NULL, y = NULL,
    title = "Thermal limits across three taxa: CTmax and z",
    subtitle = "Confidence Eyes: pale lens = 95% confidence interval; hollow point = estimate.",
    caption = paste(
      "Profile-likelihood confidence intervals (freqTLS, no Stan).",
      "CTmax reference exposure differs by study (annotated): not a common time scale."
    )
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(
    panel.grid.major.y = ggplot2::element_blank(),
    panel.grid.minor = ggplot2::element_blank(),
    plot.caption = ggplot2::element_text(hjust = 0)
  )
```

## Within-taxon contrasts

The two grouped studies invite a within-taxon comparison: do the zebrafish life
stages differ, and do the fly sexes differ? `freqTLS` requests a profile interval
**on the difference** (`dCTmax:A-B`, `dz:A-B`) and uses its documented bootstrap
fallback when a contrast profile does not close. Both are prior-free frequentist
counterparts to the `bayesTLS` pMCMC bracket. A difference interval that excludes
zero is a clear separation; one that spans zero is not.

```{r contrasts}
contrast_tbl <- summary_cache$contrasts
knitr::kable(
  data.frame(
    Contrast   = contrast_tbl$parameter,
    Difference = round(contrast_tbl$estimate, 3),
    `95% CI`   = sprintf("[%.3f, %.3f]",
                         contrast_tbl$conf.low, contrast_tbl$conf.high),
    `Excludes 0` = ifelse(
      contrast_tbl$conf.low > 0 | contrast_tbl$conf.high < 0, "yes", "no"
    ),
    Method = contrast_tbl$method,
    check.names = FALSE
  ),
  caption = "Within-taxon contrasts: 95% confidence intervals from the requested profile or its documented bootstrap fallback."
)
```

For zebrafish, the one clear separation in `CTmax` is **old embryos versus young
embryos** and **larvae versus old embryos**: old embryos sit about 1.46 degrees
Celsius above young embryos, with a confidence interval well clear of zero,
while larvae are essentially indistinguishable from young embryos in `CTmax`. None
of the `z` contrasts excludes zero, so the three stages share a common duration
sensitivity even where their critical temperatures differ. For *D. suzukii*, both
the `CTmax` and the `z` sex contrasts span zero: there is no clear sex difference
in either thermal limit, the same conclusion Ørsted et al. (2024) reached for these
data.

## What this shows

Two things stand out from the panel:

* **Resolving a taxon by an internal axis can matter or not.** Zebrafish life stage
  shifts `CTmax` by over a degree (old embryos are the most heat-tolerant stage),
  whereas *D. suzukii* sex shifts neither `CTmax` nor `z` detectably. The
  Confidence Eyes make this visible at a glance: the zebrafish lenses separate on
  the `CTmax` axis, while the two fly lenses overlap almost completely.

* **Every estimate carries an honest, prior-free interval.** All 12 headline
  profiles close (six groups for each of two parameters), so each row is a closed
  lens rather than a hollow point. Where a profile
  did **not** close (a weakly identified design or a boundary asymptote), the same
  display would show a hollow point with no lens -- never a fabricated closed eye.

These are confidence intervals throughout: they summarise the
likelihood's support for each parameter and make no probability statement about the
parameter itself.

## See bayesTLS for the Bayesian multi-taxon ridge

This panel is the `freqTLS` (likelihood) counterpart to Manuscript Figure 5 of
the `bayesTLS` supplement. Its full figure also includes the snow-gum PSII
dataset, which is not redistributed here because its source is CC BY-NC 4.0.
The Bayesian supplement draws the taxa as posterior **ridge
densities** with median points and 95% credible bars. The two displays answer the
same cross-study question from complementary inferential engines: the `bayesTLS`
ridge is a posterior density shaped by priors and the data, whereas the
`freqTLS` Confidence Eye is a prior-free likelihood interval that cannot be
read as a posterior. For the Bayesian multi-taxon ridge plot, the within-taxon
pMCMC contrasts, and the full posterior workflow, see `bayesTLS`
(<https://github.com/daniel1noble/bayesTLS>); for the side-by-side
posterior-versus-Confidence-Eye contrast on a single dataset, see
`vignette("comparing-to-bayesTLS")`.
