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

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

## ----profile-shrimp-----------------------------------------------------------
data(shrimp_lethal)
shrimp_std <- standardize_data(
  shrimp_lethal,
  temp = "Temperature_assay", duration = "Duration_exposure_hours",
  n_total = "N_individuals_after_trial", mortality = "Mortality_after_trial",
  duration_unit = "hours"
)
shrimp_fit <- fit_4pl(shrimp_std, t_ref = 1, family = "beta_binomial", quiet = TRUE)
tls(shrimp_fit, method = "profile")$summary

## ----profile-zebrafish--------------------------------------------------------
data(zebrafish_lethal)
zebra_std <- standardize_data(
  zebrafish_lethal,
  temp = "assay_temp", duration = "duration_h",
  n_total = "n_total", n_surv = "n_surv", duration_unit = "hours"
)
zebra_fit <- suppressWarnings(fit_4pl(zebra_std, by = "life_stage",
                                      t_ref = 1, family = "beta_binomial", quiet = TRUE))
# per-stage CTmax and z with profile intervals
tls(zebra_fit, method = "profile")$summary

## ----three-way, echo = FALSE--------------------------------------------------
cache_path <- system.file("extdata", "bayesTLS_benchmark_cache.rds", package = "freqTLS")
if (!nzchar(cache_path) || !file.exists(cache_path)) {
  stop("The shipped bayesTLS benchmark cache is missing; reinstall freqTLS from a complete source tarball.")
}
cache <- readRDS(cache_path)

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

pro <- confint(shrimp_fit, c("CTmax", "z"), method = "profile")  # confint() delegates to the engine fit

pick <- function(df, parm, cols) {
  r <- df[df$dataset == "shrimp" & df$parameter == parm, ]
  stats::setNames(as.numeric(r[cols]), c("est", "lo", "hi"))
}
ts_ct <- pick(cache$two_stage, "CTmax", c("estimate", "lower", "upper"))
ts_z  <- pick(cache$two_stage, "z",     c("estimate", "lower", "upper"))
by_ct <- pick(cache$bayesian,  "CTmax", c("median",   "lower", "upper"))
by_z  <- pick(cache$bayesian,  "z",     c("median",   "lower", "upper"))
pr_ct <- pro[pro$parameter == "CTmax", ]
pr_z  <- pro[pro$parameter == "z", ]

knitr::kable(data.frame(
  Quantity = c("CTmax (°C)", "z (°C / decade)"),
  `Two-stage (delta CI)`    = c(fmt(ts_ct["est"], ts_ct["lo"], ts_ct["hi"]),
                                fmt(ts_z["est"],  ts_z["lo"],  ts_z["hi"])),
  `bayesTLS (95% CrI)`      = c(fmt(by_ct["est"], by_ct["lo"], by_ct["hi"]),
                                fmt(by_z["est"],  by_z["lo"],  by_z["hi"])),
  `freqTLS (profile CI)` = c(fmt(pr_ct$estimate, pr_ct$conf.low, pr_ct$conf.high),
                                fmt(pr_z$estimate,  pr_z$conf.low,  pr_z$conf.high)),
  check.names = FALSE
), caption = "Shrimp CTmax and z shown side by side: the two model fits use the relative midpoint; the classical estimator uses absolute LT50.")

## ----benchmark-provenance, echo = FALSE---------------------------------------
provenance <- data.frame(
  Field = c("bayesTLS version", "bayesTLS source commit", "CmdStan version",
            "Cache build date", "Seed", "Model threshold", "Shape design"),
  Value = c(
    cache$meta$bayesTLS_version,
    cache$meta$git_sha,
    cache$meta$cmdstan_version,
    cache$meta$date_built,
    cache$meta$seed,
    cache$meta$config$target_surv,
    cache$meta$config$temp_effects
  )
)
knitr::kable(provenance, caption = "Version-stamped benchmark-cache provenance.")

## ----perf-load, echo = FALSE--------------------------------------------------
perf_path <- system.file("extdata", "performance_results.rds", package = "freqTLS")
has_perf  <- nzchar(perf_path) && file.exists(perf_path)
if (has_perf) perf <- readRDS(perf_path)

## ----perf-speed, echo = FALSE-------------------------------------------------
if (has_perf) {
  knitr::kable(
    perf$speed[, c("family", "design", "n_obs", "median_fit_ms", "median_profile_ms")],
    caption = "Median wall-clock (ms): one fit, and one CTmax profile CI."
  )
}

## ----speed-head-to-head, echo = FALSE-----------------------------------------
data("shrimp_lethal", package = "freqTLS", envir = environment())
fmt_t <- function(s) if (is.na(s)) "--" else if (s < 1) sprintf("%.0f ms", s * 1000) else sprintf("%.1f s", s)

sp_std <- standardize_data(shrimp_lethal, temp = "Temperature_assay",
                           duration = "Duration_exposure_hours",
                           n_total = "N_individuals_after_trial",
                           mortality = "Mortality_after_trial", duration_unit = "hours")
t_fit <- system.time(
  sp_fit <- fit_4pl(sp_std, t_ref = 1, family = "beta_binomial", quiet = TRUE)
)[["elapsed"]]
t_wald <- system.time(invisible(confint(sp_fit, c("CTmax", "z"), method = "wald")))[["elapsed"]]
t_prof <- system.time(invisible(confint(sp_fit, c("CTmax", "z"), method = "profile")))[["elapsed"]]

speed <- data.frame(
  Estimator = c("freqTLS", "freqTLS", "freqTLS", "classical two-stage", "bayesTLS"),
  Task = c("fit (ML)", "fit + Wald CTmax & z", "fit + profile CTmax & z",
           "fit + delta CI", "fit (4 chains x 4000 MCMC)"),
  `Wall-clock` = c(fmt_t(t_fit), fmt_t(t_fit + t_wald), fmt_t(t_fit + t_prof), "--", "--"),
  Source = c("live", "live", "live", "cached", "cached"),
  check.names = FALSE, stringsAsFactors = FALSE
)
tpath <- system.file("extdata", "timing_results.rds", package = "freqTLS")
if (nzchar(tpath) && file.exists(tpath)) {
  tm <- readRDS(tpath)$timing
  speed$`Wall-clock`[speed$Estimator == "classical two-stage"] <- fmt_t(tm$seconds[tm$method == "two_stage"])
  speed$`Wall-clock`[speed$Estimator == "bayesTLS"] <- fmt_t(tm$seconds[tm$method == "bayesTLS"])
}
knitr::kable(speed, caption = "Wall-clock on the shrimp benchmark. freqTLS is timed live as this page renders; the two-stage and bayesTLS times are cached from the repository-only maintainer script data-raw/timing-study.R (not installed with the package). The bayesTLS time is post-compile sampling plus overhead, and a first fit also pays a one-time Stan compilation.")

## ----perf-acc, echo = FALSE---------------------------------------------------
if (has_perf) {
  knitr::kable(
    perf$accuracy[, c("family", "truth_setting", "parameter", "bias", "rmse", "n_converged")],
    caption = sprintf("Near-unbiased recovery (nsim = %d).", perf$meta$nsim_accuracy)
  )
}

## ----perf-cov, echo = FALSE---------------------------------------------------
if (has_perf) {
  knitr::kable(
    perf$coverage[, c("family", "method", "parameter", "coverage", "median_width", "nominal")],
    caption = sprintf("Profile and Wald 95%% interval coverage (nsim = %d).",
                      perf$meta$nsim_coverage)
  )
}

## ----bb-phi, echo = FALSE-----------------------------------------------------
bbp_path <- system.file("extdata", "beta_binomial_phi_results.rds", package = "freqTLS")
if (nzchar(bbp_path) && file.exists(bbp_path)) {
  bbp <- readRDS(bbp_path)
  z <- bbp$coverage[bbp$coverage$parameter == "z", ]
  phis <- sort(unique(z$phi_true))
  lab <- c("5" = "strong", "50" = "mild", "200" = "very mild")[as.character(phis)]
  pick <- function(m) z$coverage[match(paste(phis, m), paste(z$phi_true, z$method))]
  tab <- data.frame(
    `phi (overdispersion)` = sprintf("%g (%s)", phis, ifelse(is.na(lab), "", lab)),
    profile = pick("profile"), Wald = pick("wald"), bootstrap = pick("bootstrap"),
    check.names = FALSE)
  knitr::kable(tab, digits = 3, caption = sprintf(
    paste("Empirical 95%% coverage of z by method as overdispersion weakens (phi",
          "up = nearer binomial, phi weakly identified; %d sims/cell). The profile",
          "collapses; Wald holds; bootstrap is intermediate."),
    bbp$meta$nsim))
}

## ----bootstrap-demo, eval = FALSE---------------------------------------------
# # A deliberately sparse design where the CTmax profile does not close.
# sparse <- simulate_tls(family = "binomial", temps = c(35, 36), times = c(1, 2),
#                        reps = 2, n = 10, CTmax = 36, z = 4, seed = 9)
# sfit_sparse <- suppressWarnings(fit_tls(sparse, y = survived, n = total,
#                                         time = duration, temp = temp,
#                                         family = "binomial", tref = 1))
# # Strict profile: NA on an open side. Default: a parametric bootstrap interval.
# strict <- suppressWarnings(
#   confint(sfit_sparse, "CTmax", method = "profile", fallback = FALSE))
# boot <- suppressWarnings(
#   confint(sfit_sparse, "CTmax", method = "profile", nboot = 1000, boot_seed = 1))
# data.frame(
#   setting   = c("strict profile (fallback = FALSE)", "default (bootstrap fallback)"),
#   conf.low  = c(strict$conf.low,  boot$conf.low),
#   conf.high = c(strict$conf.high, boot$conf.high),
#   method    = c(strict$method,    boot$method)
# )

## ----eye, fig.alt = "Confidence Eye for the shrimp CTmax and z: pale confidence lenses with hollow point estimates, the freqTLS uncertainty display."----
plot_confidence_eye(shrimp_fit, parm = c("CTmax", "z"), method = "profile")

## ----v02-stage-shape----------------------------------------------------------
stage_shape <- suppressWarnings(fit_4pl(
  zebra_std, by = "life_stage",
  low = ~ life_stage, up = ~ life_stage, k = ~ life_stage,
  t_ref = 1, family = "beta_binomial", quiet = TRUE
))
# zebra_fit (above) is the shared-shape fit; stage_shape lets low / up / k vary.
c(shared_shape_AIC = round(AIC(zebra_fit), 1),
  stage_shape_AIC  = round(AIC(stage_shape), 1))

## ----v02-up-------------------------------------------------------------------
stage_est <- stage_shape$fit$estimates
stage_est[grepl("^up:", stage_est$parameter), c("parameter", "estimate", "std.error")]

## ----v02-derive---------------------------------------------------------------
c(
  CTmax_50pct_1h = round(derive_ctmax(shrimp_fit$fit, surv = 0.5, duration = 1), 2),
  surv_32C_4h    = round(tail(predict_heat_injury(shrimp_fit$fit,
                     data.frame(time = seq(0, 4, by = 0.1), temp = 32))$survival, 1), 3)
)