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

## ----setup--------------------------------------------------------------------
library(freqTLS)

## ----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")]

## ----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."
)

## ----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.")
}

## ----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)
  )

## ----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."
)