Cross-checking freqTLS with bayesTLS • freqTLS

The packages teach the same canonical empirical examples but answer them with different inference engines. bayesTLS uses posterior sampling; freqTLS uses maximum likelihood with profile, Wald, or bootstrap confidence intervals. Neither package’s uncertainty object is a drop-in replacement for the other.

Experimental-software warning: this comparison does not validate freqTLS. Users remain responsible for checking the data, design, model specification, convergence, identifiability, diagnostics, and interpretation. Independently refit important analyses with bayesTLS. Agreement is a cross-check, not proof of correctness; shared mistakes can make both packages agree.

Comparison scope

A numerical difference is interpretable only when all of the following are identical:

exact dataset bytes and filter;
response endpoint and family;
every ctmax, z, low, up, and k formula;
grouping and random-effect structure;
resolved reference-time unit;
relative or absolute threshold and reported estimand.

The comparison uses the pinned bayesTLS supplement rendered 2026-07-14 from commit 76510412.

Reproducibility record

The curated cache records the exact analysis hashes, formulas, filters, thresholds, reference times, source versions, and sampler diagnostics. The full rebuild and publication procedure is in docs/design/48-canonical-comparator-cache.md.

data.frame(
  item = c(
    "bayesTLS", "freqTLS cache-build source", "CmdStan", "R", "built (UTC)",
    "OpenBLAS threads", "bounded cores"
  ),
  recorded_value = c(
    paste0(bayes_cache$meta$bayesTLS_version, " @ ",
           substr(bayes_cache$meta$bayesTLS_git_sha, 1, 12)),
    paste0("pre-release source @ ",
           substr(bayes_cache$meta$freqTLS_git_sha, 1, 12)),
    bayes_cache$meta$cmdstan_version,
    bayes_cache$meta$R_version,
    bayes_cache$meta$date_built_utc,
    bayes_cache$meta$openblas_num_threads,
    bayes_cache$meta$bounded_cores
  ),
  check.names = FALSE
)
#>                         item                    recorded_value
#> 1                   bayesTLS              1.0.0 @ 76510412e06c
#> 2 freqTLS cache-build source pre-release source @ b32c86001a7e
#> 3                    CmdStan                            2.39.0
#> 4                          R                             4.5.3
#> 5                built (UTC)           2026-07-16 15:48:18 UTC
#> 6           OpenBLAS threads                                 1
#> 7              bounded cores                                 4

Every Bayesian fit passed the recorded R-hat, effective-sample-size, divergence, tree-depth, and BFMI gates.

bayes_cache$diagnostics[c(
  "case_id", "rhat_max", "ess_bulk_min", "ess_tail_min", "divergences",
  "treedepth_hits", "bfmi_min", "all_pass"
)]
#>                case_id rhat_max ess_bulk_min ess_tail_min divergences
#> 1     zebrafish_oxygen   1.0017         2507         1750           0
#> 2           aphid_age6   1.0016         3546         4500           0
#> 3        aphid_all_age   1.0019         2205         3027           0
#> 4         snowgum_psii   1.0009         3432         4687           0
#> 5 drosophila_mortality   1.0014         2822         4680           0
#> 6     drosophila_awake   1.0013         2964         4053           0
#>   treedepth_hits bfmi_min all_pass
#> 1              0   0.8577     TRUE
#> 2              0   0.9824     TRUE
#> 3              0   0.9556     TRUE
#> 4              0   0.8033     TRUE
#> 5              0   0.8929     TRUE
#> 6              0   0.9163     TRUE

Refit the same frequentist specifications

The code below refits the exact current freqTLS specifications from the installed datasets. The fitting output is suppressed because the diagnostic table follows, but the code remains visible and copyable. The five case articles explain each preparation and fit in smaller steps.

For the awake/coma count endpoint, duration-zero controls are excluded because the 4PL uses log10(duration) and cannot represent zero duration—not because they are evidence against up. Short positive exposures or mild treatments are needed to inform that upper asymptote.

data(zebrafish_o2)
zf <- droplevels(subset(
  zebrafish_o2,
  ploidy == "diploid" & oxygen %in% c("normoxia", "hyperoxia")
))
zf_fit <- fit_4pl(
  standardize_data(
    zf, temp = "temp", duration = "duration_min",
    n_total = "n_total", n_surv = "n_surv", duration_unit = "minutes"
  ),
  ctmax = ~ 0 + oxygen, z = ~ 0 + oxygen, low = ~ 0 + oxygen,
  up = ~ 1, k = ~ 1, family = "beta_binomial", t_ref = 60,
  method = "wald", quiet = TRUE
)

data(aphid_tdt)
aphid6 <- droplevels(subset(aphid_tdt, branch == "heat" & age == "6"))
aphid6_fit <- fit_4pl(
  standardize_data(
    aphid6, temp = "temp", duration = "duration_min",
    n_total = "n_total", n_surv = "n_surv", duration_unit = "minutes"
  ),
  ctmax = ~ 0 + species, z = ~ 0 + species, low = ~ 1, up = ~ 1,
  k = ~ temp_c, family = "beta_binomial", t_ref = 60,
  method = "wald", quiet = TRUE
)

aphid_all <- droplevels(subset(aphid_tdt, branch == "heat"))
aphid_all_std <- standardize_data(
  aphid_all, temp = "temp", duration = "duration_min",
  n_total = "n_total", n_surv = "n_surv", duration_unit = "minutes"
)
aphid_all_fit <- fit_4pl(
  aphid_all_std,
  ctmax = ~ 1 + species * age, z = ~ 1 + species * age,
  low = ~ 1, up = ~ 1, k = ~ temp_c, family = "beta_binomial",
  t_ref = 60, method = "wald", quiet = TRUE
)

data(snowgum_psii)
snowgum_fit <- fit_4pl(
  standardize_data(
    snowgum_psii, temp = "Temp", duration = "Time",
    proportion = "fvfm_prop", duration_unit = "minutes"
  ),
  ctmax = ~ 0 + recovery + (1 | plant), z = ~ 0 + recovery,
  low = ~ 1, up = ~ 1, k = ~ 1, family = "beta", t_ref = 60,
  method = "wald", quiet = TRUE
)

data(dsuzukii)
mort <- stats::aggregate(
  cbind(
    n_dead = as.integer(dsuzukii$dead),
    n_total = rep.int(1L, nrow(dsuzukii))
  ) ~ temp + time + sex,
  data = dsuzukii,
  FUN = sum
)
mort$n_surv <- mort$n_total - mort$n_dead
mort_fit <- fit_4pl(
  standardize_data(
    mort, temp = "temp", duration = "time",
    n_total = "n_total", n_surv = "n_surv", duration_unit = "minutes"
  ),
  ctmax = ~ 0 + sex, z = ~ 0 + sex,
  low = ~ temp_c, up = ~ temp_c, k = ~ temp_c,
  family = "beta_binomial", t_ref = 240, method = "wald", quiet = TRUE
)

coma_cell <- interaction(
  dsuzukii$temp, dsuzukii$lvl, dsuzukii$sex, drop = TRUE
)
coma <- do.call(rbind, lapply(split(dsuzukii, coma_cell), function(d) {
  data.frame(
    temp = d$temp[1], lvl = d$lvl[1], sex = d$sex[1],
    duration = d$time[1], n_total = nrow(d),
    n_awake = sum(is.na(d$t_coma))
  )
}))
coma <- droplevels(subset(coma, duration > 0))
coma_fit <- fit_4pl(
  standardize_data(
    coma, temp = "temp", duration = "duration",
    n_total = "n_total", n_surv = "n_awake", duration_unit = "minutes"
  ),
  ctmax = ~ sex, z = ~ sex,
  low = ~ temp_c, up = ~ temp_c, k = ~ temp_c,
  family = "beta_binomial", t_ref = 60, method = "wald", quiet = TRUE
)

All displayed ML fits converged with a positive-definite Hessian and passed the package’s raw-gradient gate. The interval method for the direct-coordinate rows below is Wald; the all-age aphid cell values and mortality absolute-LT50 values are point predictions and are labelled accordingly. The organism-specific articles show the shape-identification warnings: temperature-varying low, up, or k can be weakly identified, so readers should inspect check_tls() before interpreting a shape effect.

fit_list <- list(
  zebrafish_oxygen = zf_fit,
  aphid_age6 = aphid6_fit,
  aphid_all_age = aphid_all_fit,
  snowgum_psii = snowgum_fit,
  drosophila_mortality = mort_fit,
  drosophila_awake = coma_fit
)
freq_diagnostics <- do.call(rbind, lapply(names(fit_list), function(case_id) {
  ans <- diagnose_tdt_fit(fit_list[[case_id]])
  data.frame(
    case_id = case_id,
    converged = ans$converged,
    pd_hessian = ans$pd_hessian,
    max_abs_gradient = ans$max_abs_gradient,
    gradient_pass = ans$gradient_pass,
    interval_method = if (case_id == "aphid_all_age") {
      "ML point prediction"
    } else if (case_id == "drosophila_mortality") {
      "Wald direct coordinates; absolute LT50 point only"
    } else {
      "Wald 95% confidence interval"
    }
  )
}))
rownames(freq_diagnostics) <- NULL
freq_diagnostics
#>                case_id converged pd_hessian max_abs_gradient gradient_pass
#> 1     zebrafish_oxygen      TRUE       TRUE     1.948692e-04          TRUE
#> 2           aphid_age6      TRUE       TRUE     2.091929e-06          TRUE
#> 3        aphid_all_age      TRUE       TRUE     9.832409e-06          TRUE
#> 4         snowgum_psii      TRUE       TRUE     2.898515e-04          TRUE
#> 5 drosophila_mortality      TRUE       TRUE     7.996796e-05          TRUE
#> 6     drosophila_awake      TRUE       TRUE     1.204766e-05          TRUE
#>                                     interval_method
#> 1                      Wald 95% confidence interval
#> 2                      Wald 95% confidence interval
#> 3                               ML point prediction
#> 4                      Wald 95% confidence interval
#> 5 Wald direct coordinates; absolute LT50 point only
#> 6                      Wald 95% confidence interval

Actual paired differences

The primary table puts the maximum-likelihood estimate and its Wald confidence interval beside the Bayesian posterior median and credible interval. Interval widths are descriptive, method-specific quantities: a confidence-interval width and credible-interval width are not interchangeable measures of uncertainty.

as_freq_rows <- function(fit, case_id, by, endpoint, t_ref) {
  out <- tls(fit, by = by, lethal = FALSE, method = "wald")$summary
  data.frame(
    case_id = case_id,
    endpoint = endpoint,
    threshold = "relative",
    t_ref = t_ref,
    group = as.character(out[[by]]),
    parameter = out$quantity,
    freq_estimate = out$median,
    freq_lower = out$lower,
    freq_upper = out$upper,
    freq_interval = "Wald 95% confidence interval"
  )
}

freq_rows <- rbind(
  as_freq_rows(
    zf_fit, "zebrafish_oxygen", "oxygen", "survival counts", 60
  ),
  as_freq_rows(
    aphid6_fit, "aphid_age6", "species", "survival counts", 60
  ),
  as_freq_rows(
    snowgum_fit, "snowgum_psii", "recovery",
    "retained PSII function proportion", 60
  ),
  as_freq_rows(
    coma_fit, "drosophila_awake", "sex",
    "awake counts (no recorded t_coma)", 60
  )
)

bayes_rows <- bayes_cache$summaries
# Use the current reader-facing wording before joining the pinned cache to the
# freshly refitted rows. The cache's older "missing t_coma" label denotes the
# same awake/coma endpoint.
bayes_rows$endpoint[bayes_rows$case_id == "drosophila_awake"] <-
  "awake counts (no recorded t_coma)"
bayes_rows$group <- ifelse(
  !is.na(bayes_rows$oxygen), as.character(bayes_rows$oxygen),
  ifelse(
    !is.na(bayes_rows$recovery), as.character(bayes_rows$recovery),
    ifelse(
      !is.na(bayes_rows$sex), as.character(bayes_rows$sex),
      as.character(bayes_rows$species)
    )
  )
)

paired <- merge(
  freq_rows,
  bayes_rows[c(
    "case_id", "endpoint", "threshold", "t_ref", "group", "parameter",
    "median", "lower", "upper", "interval_method"
  )],
  by = c("case_id", "endpoint", "threshold", "t_ref", "group", "parameter"),
  all = FALSE,
  sort = FALSE
)
names(paired)[names(paired) == "median"] <- "bayes_median"
names(paired)[names(paired) == "lower"] <- "bayes_lower"
names(paired)[names(paired) == "upper"] <- "bayes_upper"
names(paired)[names(paired) == "interval_method"] <- "bayes_interval"
paired$difference_freq_minus_bayes <- paired$freq_estimate - paired$bayes_median
paired$freq_width <- paired$freq_upper - paired$freq_lower
paired$bayes_width <- paired$bayes_upper - paired$bayes_lower
stopifnot(nrow(paired) == 18L)

paired_display <- paired[c(
  "case_id", "endpoint", "threshold", "t_ref", "group", "parameter",
  "freq_estimate", "freq_lower", "freq_upper", "bayes_median",
  "bayes_lower", "bayes_upper", "freq_width", "bayes_width",
  "difference_freq_minus_bayes"
)]
names(paired_display) <- c(
  "Case", "Endpoint", "Threshold", "tref", "Group", "Quantity",
  "ML estimate", "ML CI lower", "ML CI upper", "Bayes median",
  "Bayes CrI lower", "Bayes CrI upper", "ML CI width", "Bayes CrI width",
  "ML minus Bayes"
)
knitr::kable(paired_display, digits = 3,
             caption = "Matched quantities: ML Wald confidence intervals and Bayesian credible intervals remain labelled separately.")

Matched quantities: ML Wald confidence intervals and Bayesian credible intervals remain labelled separately.
Case	Endpoint	Threshold	tref	Group	Quantity	ML estimate	ML CI lower	ML CI upper	Bayes median	Bayes CrI lower	Bayes CrI upper	ML CI width	Bayes CrI width	ML minus Bayes
zebrafish_oxygen	survival counts	relative	60	normoxia	CTmax	38.740	38.432	39.048	38.889	38.558	39.375	0.616	0.817	-0.150
zebrafish_oxygen	survival counts	relative	60	hyperoxia	CTmax	39.301	39.060	39.543	39.348	39.170	39.499	0.483	0.330	-0.047
zebrafish_oxygen	survival counts	relative	60	normoxia	z	6.006	4.371	8.251	5.587	4.056	7.754	3.880	3.698	0.419
zebrafish_oxygen	survival counts	relative	60	hyperoxia	z	2.501	2.139	2.923	2.488	2.048	2.933	0.784	0.885	0.013
aphid_age6	survival counts	relative	60	M_dirhodum	CTmax	35.215	35.025	35.405	35.173	34.972	35.353	0.380	0.381	0.043
aphid_age6	survival counts	relative	60	S_avenae	CTmax	36.527	36.431	36.623	36.509	36.410	36.600	0.192	0.189	0.019
aphid_age6	survival counts	relative	60	R_padi	CTmax	37.170	37.063	37.277	37.145	37.044	37.251	0.213	0.207	0.025
aphid_age6	survival counts	relative	60	M_dirhodum	z	4.747	4.505	5.002	4.775	4.530	5.043	0.497	0.514	-0.028
aphid_age6	survival counts	relative	60	S_avenae	z	3.610	3.461	3.766	3.624	3.471	3.781	0.305	0.310	-0.014
aphid_age6	survival counts	relative	60	R_padi	z	3.965	3.698	4.252	3.959	3.686	4.247	0.554	0.562	0.007
snowgum_psii	retained PSII function proportion	relative	60	Dark	CTmax	45.721	45.178	46.264	45.695	45.046	46.291	1.085	1.245	0.026
snowgum_psii	retained PSII function proportion	relative	60	Light	CTmax	44.075	43.558	44.593	44.054	43.433	44.698	1.035	1.265	0.021
snowgum_psii	retained PSII function proportion	relative	60	Dark	z	4.707	4.290	5.165	4.674	4.230	5.122	0.875	0.893	0.033
snowgum_psii	retained PSII function proportion	relative	60	Light	z	3.641	3.202	4.140	3.624	3.178	4.106	0.938	0.927	0.016
drosophila_awake	awake counts (no recorded t_coma)	relative	60	F	CTmax	36.497	36.393	36.601	36.480	36.363	36.587	0.208	0.224	0.017
drosophila_awake	awake counts (no recorded t_coma)	relative	60	M	CTmax	36.296	36.203	36.389	36.270	36.159	36.372	0.186	0.214	0.026
drosophila_awake	awake counts (no recorded t_coma)	relative	60	F	z	2.433	2.288	2.587	2.413	2.254	2.581	0.299	0.328	0.019
drosophila_awake	awake counts (no recorded t_coma)	relative	60	M	z	2.393	2.262	2.532	2.404	2.259	2.569	0.270	0.310	-0.011

The age-six aphid row is also one cell of the all-age extension, but the two fits are intentionally distinct. For all nine species-by-age cells, the largest absolute freqTLS-minus-bayesTLS point difference is reported separately rather than printing another 18-row table.

aphid_cells <- expand.grid(
  species = levels(aphid_all$species),
  age = levels(aphid_all$age),
  temp = mean(aphid_all_std$temp),
  KEEP.OUT.ATTRS = FALSE
)
aphid_pred <- predict(aphid_all_fit, aphid_cells, type = "parameters")
aphid_freq <- rbind(
  data.frame(
    species = aphid_cells$species, age = aphid_cells$age,
    parameter = "CTmax", freq_estimate = aphid_pred$CTmax
  ),
  data.frame(
    species = aphid_cells$species, age = aphid_cells$age,
    parameter = "z", freq_estimate = aphid_pred$z
  )
)
aphid_bayes <- subset(bayes_rows, case_id == "aphid_all_age")
aphid_join <- merge(
  aphid_freq,
  aphid_bayes[c("species", "age", "parameter", "median")],
  by = c("species", "age", "parameter")
)
aphid_join$difference_freq_minus_bayes <-
  aphid_join$freq_estimate - aphid_join$median
stopifnot(nrow(aphid_join) == 18L)
aggregate(
  abs(difference_freq_minus_bayes) ~ parameter,
  data = aphid_join,
  FUN = max
)
#>   parameter abs(difference_freq_minus_bayes)
#> 1     CTmax                       0.02725888
#> 2         z                       0.02246847

Snow-gum is a paired refit of the locked shared-shape analogue. It is not a claim that freqTLS reproduces the richer recovery-by-temperature shape model displayed in the pinned supplement.

Drosophila mortality: match the reported estimand

The pinned mortality summary uses the absolute 50% threshold. freqTLS fits the relative direct-coordinate backbone, then can derive the same absolute 240-minute LT50 by solving the fitted survival curve. This article compares that shared point estimand; it does not relabel relative CTmax or z as absolute quantities.

mort_lt50 <- data.frame(
  group = mort_fit$fit$group_levels,
  freq_estimate = vapply(mort_fit$fit$group_levels, function(sex_level) {
    temp_grid <- seq(min(mort$temp), max(mort$temp), length.out = 1001)
    grid_difference <- predict(
      mort_fit,
      data.frame(temp = temp_grid, duration = 240, group = sex_level),
      type = "survival"
    ) - 0.5
    crossing <- which(
      grid_difference[-length(grid_difference)] *
        grid_difference[-1L] <= 0
    )
    stopifnot(length(crossing) == 1L)
    objective <- function(temp) {
      predict(
        mort_fit,
        data.frame(temp = temp, duration = 240, group = sex_level),
        type = "survival"
      ) - 0.5
    }
    uniroot(
      objective,
      temp_grid[c(crossing, crossing + 1L)],
      tol = 1e-10
    )$root
  }, numeric(1))
)
mort_bayes <- subset(
  bayes_rows,
  case_id == "drosophila_mortality" & parameter == "CTmax"
)
mort_compare <- merge(
  mort_lt50,
  mort_bayes[c("group", "median", "lower", "upper")],
  by = "group"
)
names(mort_compare)[names(mort_compare) == "median"] <- "bayes_median"
mort_compare$difference_freq_minus_bayes <-
  mort_compare$freq_estimate - mort_compare$bayes_median
mort_compare$freq_interval <- "none: exact-model bootstrap unstable"
stopifnot(nrow(mort_compare) == 2L)
mort_compare
#>   group freq_estimate bayes_median    lower    upper
#> 1     F      35.15584     35.13553 35.01181 35.24284
#> 2     M      35.17396     35.16390 35.04825 35.26617
#>   difference_freq_minus_bayes                        freq_interval
#> 1                  0.02031098 none: exact-model bootstrap unstable
#> 2                  0.01005650 none: exact-model bootstrap unstable

The Bayesian cache also contains absolute-threshold z. In this temperature- varying-shape model, freqTLS does not return a temperature-specific local absolute z: its direct z is the relative-midpoint backbone sensitivity. The mortality case article reports that relative z explicitly and limits the absolute comparison to the four-hour LT50 point.

What agreement means

Maximum-likelihood estimates and posterior medians need not be numerically identical. Priors, finite-sample skew, boundary behaviour, random effects, and weak identification can produce legitimate differences. The audit confirms the data and model specification first, reports the actual difference, and never averages estimates or hides a discrepancy.

Frequentist extensions and unsupported analyses

Confidence Eyes, open-profile status, likelihood contrasts, bootstrap failure provenance, and explicit convergence/Hessian/gradient diagnostics are freqTLS distinctions. Simulations, profile theory, survival surfaces, limited random intercepts, and deterministic heat-injury scenarios remain separate extension pages with synthetic examples.

Censored-time, hurdle-productivity, posterior inference, and fitted repair dynamics remain bayesTLS-only.