Comparing freqTLS to bayesTLS

Credit and origins

The thermal-load-sensitivity (TLS) modelling framework implemented here was introduced by Daniel W. A. Noble, Pieter A. Arnold, and Patrice Pottier in the bayesTLS package (manuscript in preparation). The 4PL thermal death-time model and the mapping from the midpoint slope to z and CTmax are theirs. freqTLS is an independent likelihood implementation of that framework; please cite bayesTLS when you use it. The bundled shrimp_lethal and zebrafish_lethal datasets are vendored from bayesTLS under CC BY 4.0 (see ?shrimp_lethal, ?zebrafish_lethal, and citation("freqTLS")).

freqTLS contributes a TMB maximum-likelihood likelihood, the direct CTmax/z reparameterisation that makes both quantities profile-able, and profile-likelihood confidence intervals.

The three-way design

The benchmark compares three estimators of the same constant-shape 4PL with a matched time unit and reference time. The two model fits use the relative midpoint threshold. The classical two-stage estimator uses absolute LT50; it is shown beside them only for lethal datasets whose fitted asymptotes are near zero and one, where the two thresholds are close. This is the fairness boundary used throughout this article:

Under this matched configuration the Bayesian and profile-likelihood fits target the same likelihood and the same fitted curve (see vignette("model-math")); they differ in how they summarise uncertainty — a posterior (with priors) versus a prior-free likelihood interval.

The reproducible recipe (not run here)

This is the exact bayesTLS recipe used to build the cache on a machine with Stan installed. It is a plain fenced code block, not an executable vignette chunk.

library(bayesTLS)
data(shrimp_lethal, package = "freqTLS")

# 1. Standardise: name the temperature / duration / count columns and the time
#    unit. fit_4pl() and the two-stage path both consume the standardised frame.
std <- bayesTLS::standardize_data(
  data          = shrimp_lethal,
  temp          = "Temperature_assay",
  duration      = "Duration_exposure_hours",
  n_total       = "N_individuals_after_trial",
  mortality     = "Mortality_after_trial",
  duration_unit = "hours"
)

# 2. Bayesian fit, matched: constant shape via temp_effects = "mid", beta-binomial.
bfit <- bayesTLS::fit_4pl(
  data         = std,
  temp_effects = "mid",
  family       = brms::brmsfamily("beta_binomial", link = "identity"),
  chains       = 4, iter = 4000, seed = 123, backend = "cmdstanr"
)

# 3. Relative-threshold CTmax + z at tref = 1 hour. NOTE: t_ref / time_multiplier
#    live on extract_tdt() (and ts_stage2 / ts_ci), NOT on fit_4pl().
btdt  <- bayesTLS::extract_tdt(
  bfit, target_surv = "relative",
  t_ref = 1, time_multiplier = 1, output_time_unit = "hours"
)
ctmax <- bayesTLS::get_ctmax_summary(btdt)  # temp_lower / temp_median / temp_upper
zsumm <- bayesTLS::get_z_summary(btdt)      # z_median / z_lower / z_upper

# 4. Classical two-stage, same standardised data and reference time.
s1  <- bayesTLS::ts_stage1(std, family = "betabinomial")
s2  <- bayesTLS::ts_stage2(s1, t_ref = 1, time_multiplier = 1)
tci <- bayesTLS::ts_ci(s2, method = "delta", level = 0.95,
                       t_ref = 1, time_multiplier = 1)  # $CTmax_1hr and $z blocks

To populate the cache, a maintainer works from a clone of the source repository and runs the repository-only script data-raw/build_benchmark_cache.R on a Stan machine, with bayesTLS installed from a pinned checkout and BAYESTLS_GIT_SHA set to its verified 40-character commit when the installed package lacks RemoteSha metadata. The builder refuses to write a cache if it cannot record that exact source commit. The script is not installed with the package; it writes inst/extdata/bayesTLS_benchmark_cache.rds (the Bayesian and two-stage summaries plus a meta provenance block: bayesTLS_version, git_sha, source_url, cmdstan_version, date_built, seed, the configuration, and the data-reconstruction note).

The freqTLS side (live)

The freqTLS fits run live and need nothing beyond this package. We use the shrimp_lethal data (ungrouped) and zebrafish_lethal (grouped by life stage). The shrimp survival counts are reconstructed from the source CSV proportions (see ?shrimp_lethal); because all three estimators share the vendored data, they inherit the same reconstructed counts.

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
#> # A tibble: 2 × 4
#>   quantity median lower upper
#>   <chr>     <dbl> <dbl> <dbl>
#> 1 CTmax     31.8  31.6  31.9 
#> 2 z          2.19  1.96  2.46

fit_4pl() returns a freq_tls workflow object (the twin of bayesTLS’s bayes_tls): it bundles the engine fit, the standardised data, and the formula. Its S3 methods — tls(), confint(), summary(), plot_confidence_eye() — delegate to the engine fit, so you call them on the workflow object directly. The underlying engine fit is also available as shrimp_fit$fit if you want to reach it.

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
#> # A tibble: 6 × 5
#>   life_stage    quantity median lower upper
#>   <chr>         <chr>     <dbl> <dbl> <dbl>
#> 1 young_embryos CTmax     39.9  39.8  40.0 
#> 2 old_embryos   CTmax     41.4  41.2  41.6 
#> 3 larvae        CTmax     39.8  39.7  39.9 
#> 4 young_embryos z          2.00  1.82  2.19
#> 5 old_embryos   z          1.80  1.53  2.16
#> 6 larvae        z          1.98  1.76  2.22

(The zebrafish fit emits a data-adequacy warning about temperatures with fewer than three durations; this is the kind of identifiability signal freqTLS surfaces explicitly. It is suppressed here only to keep the vignette output tidy.)

The three-way comparison, with real numbers

The cache holds the maintainer-built bayesTLS (posterior) and classical two-stage summaries; the freqTLS column is computed live as this page renders. The two model fits use the matched beta-binomial, relative-threshold, constant-shape configuration at tref = 1 hour, so they target the same fitted curve. The classical column uses absolute LT50 and is the approximate comparator described above.

For the shrimp data the three estimators land on essentially the same CTmax and z, with comparable interval widths: the profile-likelihood confidence interval and the Bayesian credible interval nearly coincide. That is the point of the complementary framing — under the matched configuration the likelihood and the posterior summarise the same fitted curve, one with a prior and MCMC, the other prior-free and by optimisation.

How fast, accurate, and calibrated

Accuracy and calibration below are descriptive characteristics of the likelihood path, measured by simulation (the repository-only maintainer script data-raw/performance-study.R, which is not installed with the package) — they ask whether freqTLS’s own intervals are trustworthy, not whether they beat bayesTLS (which buys priors, full posteriors, and the heat-injury sub-models). Speed, though, is one place a head-to-head is both fair and stark.

How fast — a full fit, and one CTmax profile interval, by design size:

Fits are milliseconds and a profile interval is well under a second. Putting the three estimators side by side on the shrimp benchmark makes the speed gap concrete:

freqTLS’s profile path runs in about a second — comparable to the classical two-stage and roughly an order of magnitude faster than the Bayesian MCMC fit — while its Wald path is near-instant. That speed is the likelihood path’s concrete advantage, and it is why freqTLS keeps profile the default while offering Wald as a fast opt-in for well-identified fits. The whole three-way table above is itself computed live as this page renders, with no Stan.

How accurate — bias and RMSE for CTmax and z:

How calibrated — empirical coverage of the 95% intervals:

Profile coverage tracks the nominal 95% closely; where it dips (small samples, beta-binomial z), freqTLS reports it honestly rather than hiding it.

Bootstrap fallback when a profile stays open

When a profile does not close — for example under a weakly identified design or a boundary asymptote — confint() can fall back to a prior-free parametric bootstrap. The fallback often yields an asymmetry-respecting interval, but it does not manufacture certainty: too few valid or non-degenerate bootstrap fits can still produce an unavailable bound. Set fallback = FALSE to retain the open profile explicitly.

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

When enough bootstrap refits remain stable and non-degenerate, freqTLS can provide a prior-free interval after an open profile. Sparse or boundary designs can still leave a bound unavailable; the status column reports that outcome. This 1,000-refit recipe is displayed rather than executed during package checks; run it interactively for the full comparison. The live strict-profile example and fallback recipe are also explained in vignette("profile-likelihood").

The teaching device: posterior density versus Confidence Eye

The clearest way to see the Bayesian-versus-likelihood distinction is to draw, for the same CTmax (or z):

• the bayesTLS posterior density — a probability distribution over the parameter, shaped by the prior and the data; and
• the freqTLS Confidence Eye — a confidence lens with a hollow point estimate, carrying no prior and making no probability statement about the parameter.

plot_confidence_eye(shrimp_fit, parm = c("CTmax", "z"), method = "profile")

The Confidence Eye above is the freqTLS half of that contrast. When the cache (and bayesTLS) are available, the posterior density for the same quantity can be drawn beside it; the side-by-side makes explicit that one is a posterior and the other is a likelihood confidence interval. freqTLS deliberately never renders a posterior-style density for its own intervals, and its prose uses “confidence” language, never “posterior” or “credible”.

Beyond the matched shape: stage-specific curves

The three-way comparison holds the shape constant (low, up, k shared across temperatures and groups) so all three estimators target the same curve — that is what makes the benchmark fair. freqTLS can also relax that restriction: the shape parameters low, up, and log_k may vary by a grouping factor. Re-fitting zebrafish_lethal with stage-specific shapes and comparing by AIC asks whether the life stages differ in more than thermal location (CTmax, z):

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))
#> shared_shape_AIC  stage_shape_AIC 
#>           1222.5           1187.9

The stage-specific model has the substantially lower AIC, so the data support per-stage shapes. The difference is concentrated in the upper asymptote up (the maximum survival at benign exposures), not the steepness k:

stage_est <- stage_shape$fit$estimates
stage_est[grepl("^up:", stage_est$parameter), c("parameter", "estimate", "std.error")]
#>          parameter  estimate  std.error
#> 4 up:young_embryos 0.7177178 0.03159488
#> 5   up:old_embryos 0.9180763 0.01502271
#> 6        up:larvae 0.9393117 0.02523538

Young embryos have a markedly lower survival ceiling (up near 0.7) than older embryos and larvae (near 0.9) — a real biological difference that the matched constant-shape configuration used for the freqTLS and bayesTLS model-based benchmark cannot express. The classical two-stage workflow is a separate absolute-LT50 approximation and likewise does not model a stage-specific 4PL shape.

freqTLS also reads absolute critical temperatures and predicts heat injury off the same fitted curve. For the (ungrouped) shrimp fit, the absolute critical temperature for 50% survival at one hour, and the survival predicted under a four-hour exposure at 32 °C:

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)
)
#> CTmax_50pct_1h    surv_32C_4h 
#>         31.730          0.018

(derive_tcrit() similarly returns a rate-multiplier T_crit, with an explicit lethal-endpoint caveat.) These are deterministic transforms of the fitted CTmax / z, not new fits.

When to prefer which

• Use freqTLS when you want fast, prior-free, asymmetry-respecting intervals and an explicit identifiability check. When the profile does not close it falls back to a parametric bootstrap; unstable or degenerate refits can still leave the interval unavailable.
• Use bayesTLS when you want a full Bayesian workflow, prior information, or the heat-injury and repair sub-models.

The two are complementary lenses on the same model. freqTLS flags weak identifiability explicitly and never claims the profile is universally superior; see vignette("profile-likelihood") for the non-closing behaviour and the bootstrap fallback.