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

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

## ----load-data, include = FALSE-----------------------------------------------
data(dsuzukii)
# Per-individual records -> survival counts per (temp, time, sex) cell.
.nd  <- aggregate(list(n_dead  = dsuzukii$dead), dsuzukii[c("temp", "time", "sex")], sum)
.nt  <- aggregate(list(n_total = dsuzukii$dead), dsuzukii[c("temp", "time", "sex")], length)
mort <- merge(.nd, .nt)

## ----data-shape---------------------------------------------------------------
str(mort)
table(mort$sex)

## ----fit----------------------------------------------------------------------
std <- standardize_data(
  mort, temp = "temp", duration = "time",
  n_total = "n_total", n_dead = "n_dead", duration_unit = "minutes"
)
# fit_4pl(by = "sex") gives a per-sex CTmax and z (the "sex on the midpoint"
# model); $fit is the engine fit the extractors below read directly.
fit <- fit_4pl(std, by = "sex", t_ref = 240, family = "beta_binomial",
               quiet = TRUE)$fit
fit$convergence$code   # 0 = converged

## ----per-sex------------------------------------------------------------------
confint(fit, c("CTmax:F", "CTmax:M", "z:F", "z:M"), method = "profile")[
  , c("parameter", "estimate", "conf.low", "conf.high", "method")]

## ----eye, fig.height = 5, fig.alt = "Confidence Eye for Drosophila suzukii CTmax and z, one lens per sex (female and male). The lenses for the two sexes overlap almost completely for both CTmax (near 35.2 C) and z (near 3 C per decade), each drawn as a pale confidence lens with a hollow point estimate."----
plot_confidence_eye(fit, parm = c("CTmax", "z"), method = "profile")

## ----survival, fig.height = 4.5, fig.alt = "Fitted Drosophila suzukii survival curves by temperature (34 to 38 C) and sex, with observed cell proportions overlaid. Survival falls from near one to near zero as exposure time increases, and the curves shift to shorter survival times at higher temperatures; the female and male curves are nearly indistinguishable at each temperature."----
plot_survival_curves(fit)

## ----tcrit--------------------------------------------------------------------
c(
  T_crit_F = round(derive_tcrit(fit, rate = 1, group = "F"), 2),
  T_crit_M = round(derive_tcrit(fit, rate = 1, group = "M"), 2)
)

## ----contrast-----------------------------------------------------------------
sex_contrasts <- confint(
  fit, c("dCTmax:F-M", "dz:F-M"), method = "profile",
  fallback = TRUE, nboot = 1000L, boot_seed = 20260713L
)
sex_contrasts[, c("parameter", "estimate", "conf.low", "conf.high", "method")]

## ----validation---------------------------------------------------------------
ml <- confint(fit, c("CTmax:F", "CTmax:M", "z:F", "z:M"), method = "profile")
data.frame(
  parameter   = ml$parameter,
  ML_estimate = round(ml$estimate, 2),
  published   = c(35.2, 35.2, 3.03, 3.28)
)

## ----three-way-suzukii, echo = FALSE------------------------------------------
cache_path <- system.file("extdata", "bayesTLS_benchmark_cache.rds", package = "freqTLS")
has_cache  <- nzchar(cache_path) && file.exists(cache_path)

fmt <- function(est, lo, hi) {
  ifelse(is.na(est), "--", sprintf("%.2f [%.2f, %.2f]", est, lo, hi))
}

sex_levels <- fit$group_levels
sex_label  <- c(F = "Female", M = "Male")

# Live freqTLS per-sex CTmax and z as profile-likelihood confidence
# intervals (with explicit open-profile status; no Stan).
pro <- suppressMessages(confint(
  fit,
  c(paste0("CTmax:", sex_levels), paste0("z:", sex_levels)),
  method = "profile"
))
pr_of <- function(parm, g) {
  r <- pro[pro$parameter == paste0(parm, ":", g), , drop = FALSE]
  if (nrow(r) == 0L) return(stats::setNames(rep(NA_real_, 3), c("est", "lo", "hi")))
  stats::setNames(c(r$estimate[1L], r$conf.low[1L], r$conf.high[1L]),
                  c("est", "lo", "hi"))
}

# The cache keys each sex as "dsuzukii:<sex>".
cache_has_suzukii <- FALSE
if (has_cache) {
  cache <- readRDS(cache_path)
  cache_has_suzukii <- !is.null(cache$bayesian) &&
    "dataset" %in% names(cache$bayesian) &&
    all(paste0("dsuzukii:", sex_levels) %in% cache$bayesian$dataset)
}

if (cache_has_suzukii) {
  pick <- function(df, parm, g, cols) {
    r <- df[df$dataset == paste0("dsuzukii:", g) & df$parameter == parm, , drop = FALSE]
    if (nrow(r) == 0L) return(stats::setNames(rep(NA_real_, 3), c("est", "lo", "hi")))
    stats::setNames(as.numeric(r[1L, cols]), c("est", "lo", "hi"))
  }
  one_row <- function(g, parm, quantity) {
    ts <- pick(cache$two_stage, parm, g, c("estimate", "lower", "upper"))
    by <- pick(cache$bayesian,  parm, g, c("median",   "lower", "upper"))
    pr <- pr_of(parm, g)
    data.frame(
      Sex      = sex_label[[g]],
      Quantity = quantity,
      `Two-stage (delta CI)`    = fmt(ts["est"], ts["lo"], ts["hi"]),
      `bayesTLS (95% CrI)`      = fmt(by["est"], by["lo"], by["hi"]),
      `freqTLS (profile CI)` = fmt(pr["est"], pr["lo"], pr["hi"]),
      check.names = FALSE
    )
  }
  rows <- do.call(rbind, lapply(sex_levels, function(g) rbind(
    one_row(g, "CTmax", "CTmax (°C)"),
    one_row(g, "z",     "z (°C / decade)")
  )))
  knitr::kable(rows, row.names = FALSE,
               caption = "D. suzukii CTmax and z per sex from three estimators (matched constant-shape configuration, tref = 240 min).")
} else {
  prof_tbl <- do.call(rbind, lapply(sex_levels, function(g) {
    ct <- pr_of("CTmax", g); z <- pr_of("z", g)
    data.frame(
      Sex = sex_label[[g]],
      `CTmax (°C), profile CI`      = fmt(ct["est"], ct["lo"], ct["hi"]),
      `z (°C / decade), profile CI` = fmt(z["est"], z["lo"], z["hi"]),
      check.names = FALSE
    )
  }))
  knitr::kable(prof_tbl, row.names = FALSE,
               caption = "Live freqTLS per-sex estimates (no Stan). The bayesTLS and two-stage columns appear once the maintainer-built cache carries the D. suzukii sexes.")
}

## ----session------------------------------------------------------------------
sessionInfo()