Frequentist and Bayesian thermal-load sensitivity: strengths, weaknesses, and the blurry line

library(freqTLS)

This vignette is for a thermal biologist deciding how to fit a thermal-load- sensitivity (TLS) model, not just which package. freqTLS (maximum likelihood, profile-likelihood intervals) and bayesTLS (Bayesian, posterior intervals) fit the same four-parameter logistic thermal-death-time model. They are complementary lenses on one model, and the choice between them is really a choice between two ways of handling the same hard realities: weak data, weak identifiability, and non-convergence. The aim here is an honest account of the trade-offs — including the uncomfortable fact, developed at the end, that the “prior-free” frequentist path quietly uses prior-like machinery of its own.

It builds without Stan: the freqTLS fits run live, while the Bayesian comparison and the coverage evidence are read from cached simulations rather than recomputed.

One model, two inferential philosophies

Both packages fit survival as a 4PL curve in log10(duration) whose midpoint moves with temperature through CTmax (the critical thermal maximum) and z (thermal sensitivity). They differ only in how they turn data into a statement about the parameters:

On strong data these agree closely — on the brown-shrimp benchmark the two packages’ CTmax estimates differ by only about 0.07 °C (cached in benchmark_vs_bayes.rds; the full three-way comparison, which adds the classical two-stage method, is in vignette("comparing-to-bayesTLS")). They diverge exactly where the data are weak and the prior starts to do the work — which is precisely where you most need to know what is data and what is assumption.

Priors: regularisation versus sensitivity

A prior is a double-edged tool.

freqTLS takes the opposite trade. Being prior-free, it needs nothing to defend and reports only what the data carry — but it offers no free regularisation: when the data are thin, the likelihood is flat and the interval is wide or fails to close. That is not a defect to paper over; it is information (see the next section). When you genuinely have prior knowledge and want to use it, that is a positive reason to reach for bayesTLS.

Identifiability: revealing versus masking

This is the sharpest practical difference.

When a parameter is weakly identified — the data barely constrain it — the likelihood is nearly flat along that direction. The two paths respond in opposite ways:

# Sparse design: few temperatures, so the CTmax/z slope is weakly identified.
sparse <- simulate_tls(
  temps = c(35, 36), times = c(1, 4), reps = 2, n = 20,
  CTmax = 36, z = 4, family = "binomial", seed = 11
)
fit_sparse <- suppressWarnings(fit_tls(
  sparse, y = survived, n = total, time = duration, temp = temp,
  family = "binomial", tref = 1
))
# freqTLS says so: a profile that does not close returns NA on the open side
# (here with the prior-free bootstrap fallback turned off to show the raw signal).
suppressWarnings(confint(fit_sparse, "z", method = "profile", fallback = FALSE))
#> # A tibble: 1 × 8
#>   parameter conf.low conf.high estimate level method  scale conf.status
#>   <chr>        <dbl>     <dbl>    <dbl> <dbl> <chr>   <chr> <chr>      
#> 1 z               NA        NA     2.95  0.95 profile log   open_both

Neither behaviour is “right” in the abstract. The point is that they answer different questions: what do the data alone support? versus what should I believe, given the data and my prior? freqTLS makes weak identifiability visible; the Bayesian path lets you act through it, at the cost of leaning on the prior. Reporting both — or at least being explicit about which you used — is the honest stance (Amrhein, Greenland & McShane, 2019).

Non-convergence and diagnostics

Both paths can fail, in different ways and with different tells:

freqTLS (ML) bayesTLS (MCMC)
Failure modes optimiser non-convergence (code != 0); non-positive-definite Hessian (pdHess = FALSE); a profile that does not close divergent transitions; low effective sample size; R-hat > 1.01; poor mixing
Diagnostics convergence code, pdHess, profile geometry, the identifiability warnings divergences, ESS, R-hat, trace plots
Typical fixes better starts, the built-in BFGS retry, a bootstrap fallback, or simplifying the model; unstable refits can still leave NA reparameterising, more iterations/adapt-delta, or stronger priors

Two things are worth noting. First, the Bayesian fix of last resort — strengthening the prior — is the very lever that trades identifiability for prior dependence (above). Second, when a profile does not close or the Hessian is not positive definite, confint() can fall back to a prior-free parametric bootstrap. That route often closes the interval without a prior, but too few stable refits still produce NA; the failure remains visible. The likelihood path needs no Stan toolchain, whereas an MCMC fit needs a working sampler.

Confidence versus credible intervals

The interval labels are not interchangeable. A freqTLS interval is a confidence interval: parameter values not rejected by the data at the stated level, with a coverage interpretation and no prior. A bayesTLS interval is a credible interval: a posterior-probability statement that depends on the prior. freqTLS therefore never calls its intervals “posterior” or “credible”, and its default visual (the Confidence Eye) deliberately avoids posterior-density iconography. Some statisticians argue these intervals should be renamed compatibility intervals — to discourage misreading the coverage as a probability about the parameter — and that “significance” language should be retired (Rafi & Greenland 2020; Amrhein, Greenland & McShane 2019); freqTLS keeps the familiar term “confidence interval” but shares that caution.

Does the interval actually cover? Small-sample calibration

A confidence interval is only worth the name if it contains the truth at the advertised rate. That is a checkable claim — and, unlike the Bayesian path whose calibration is entangled with the prior, the frequentist path lets us check it cleanly: simulate data from a known CTmax and z, refit, and count how often the nominal 95% interval covers the generating value.

The asymptotic theory behind a likelihood interval is exact only as the sample grows. In small samples the curvature-based standard error is a little too small and the normal / chi-squared reference a little too generous, so the naive interval comes out too narrow and under-covers. freqTLS corrects this with the Bates–Watts profile-t calibration: it refers the signed-root likelihood statistic to a t distribution with the residual degrees of freedom of the fit (df = n_obs − n_parameters) rather than to a standard normal. This widens the Wald and profile intervals precisely when the sample is small, and not at all when it is large (Bates & Watts, 1988).

The pay-off, measured over ~500 simulated datasets per cell (cached in calibration_results.rds):

95% confidence-interval coverage for CTmax (nominal 0.95), over ~500 simulated datasets per cell. Monte-Carlo SE ≈ 0.01.
Sample size Residual df Coverage (Wald, z-ref) Coverage (Wald, t-ref) Extra width
small 10 0.927 0.964 +14%
medium 35 0.946 0.964 +4%
large 100 0.970 0.970 +1%

At df ≈ 10 the asymptotic 95% interval covers only about 93% of the time; the t-calibration restores it to ~96%, at the cost of ~14% extra width. By df ≈ 100 the two references coincide (t ≈ z) and the correction costs essentially nothing. The widening is therefore self-cancelling: it pays for the small-sample optimism of the asymptotics, then steps out of the way as data accumulate.

The coverages above are for the Wald interval, where the z-versus-t reference is a clean one-line change. The default profile interval inherits the same small-sample correction — its cutoff is the squared t quantile qt(df)^2 rather than qchisq(level, 1) — so it tracks the Wald t-reference column, not the asymptotic-z one.

This bears on the Bayesian comparison in two ways. First, it is the honest basis for the “confidence” label of the previous section — here the coverage is measured, not assumed. Second, it is itself a mild small-sample adjustment, a frequentist cousin of the regularisation discussed next: a correction that leans on the t distribution much as a weakly informative prior leans on its scale. The difference is that this one is explicit, validated by simulation, and disappears as the sample grows.

The blurry line: clamping, shrinkage, and penalties as implicit priors

It is tempting to frame this as “prior versus no prior”. That is too clean. Frequentist fits routinely embed prior-like regularisation, and freqTLS is no exception. Three forms appear in this very package.

Clamping. To keep the likelihood well behaved, freqTLS clamps the fitted probability off the boundary (to [1e-12, 1 - 1e-12]), clamps a beta- family response off {0, 1}, and floors the Beta shape parameters. Each is a hard interior bound — behaviourally a degenerate prior that assigns no mass to the edge. This is the same disease that separation cures with a penalty or a weak prior: with perfectly separated binary data the unpenalised maximum- likelihood estimate runs to infinity, and a mild penalty (equivalently, a weakly informative prior) brings it back to a finite value (Firth, 1993; Gelman, Jakulin, Pittau & Su, 2008). The advantage is numerical stability and finite estimates; the cost is a small, usually hidden, boundary bias and an arbitrary clamp constant.

Shrinkage / partial pooling. freqTLS’s random intercepts on CTmax, log_z, low, and log_k (e.g. CTmax ~ <fixed> + (1 | group)) assume the group deviations are Gaussian, b_g ~ N(0, sigma). That Gaussian is a prior on the deviations, and the predicted group effects (the BLUPs returned by ranef()) are shrunk toward zero exactly as a posterior mean would be — the empirical-Bayes view of mixed models (James & Stein, 1961; Efron & Morris, 1975; Robinson, 1991). Shrinkage borrows strength across groups and lowers variance; the costs are that the shrinkage scale is itself a modelling assumption, and the maximum-likelihood variance is biased low with few groups, a boundary problem usually addressed with — again — a penalty or prior on the variance component (Chung, Rabe-Hesketh, Dorie, Gelman & Liu, 2013).

# Data-poor colonies (few assays each) so partial pooling has work to do.
d_re <- simulate_tls(
  family = "binomial", temps = c(34, 36, 38), times = c(1, 4), reps = 1, n = 8,
  CTmax = 36, z = 4, re_sd = 1.5, n_re_groups = 12, seed = 42
)
# Partial pooling: a random intercept shrinks the colony effects toward 0.
fit_re <- suppressWarnings(fit_tls(
  tls_bf(survived | trials(total) ~ time(duration) + temp(temp),
         CTmax ~ 1 + (1 | colony)),
  data = d_re, family = "binomial", tref = 1
))
blup <- ranef(fit_re)$estimate
# No pooling: a fixed CTmax per colony (the unshrunk per-group estimates).
fit_fix <- suppressWarnings(fit_tls(
  d_re, y = survived, n = total, time = duration, temp = temp,
  group = colony, family = "binomial", tref = 1
))
ct_fix <- fit_fix$estimates$estimate[startsWith(fit_fix$estimates$parameter, "CTmax:")]
# The Gaussian prior on b_g (scale sigma_CTmax) pulls the spread in.
c(no_pooling_sd = round(sd(ct_fix - mean(ct_fix)), 2),
  partial_pooling_sd = round(sd(blup), 2),
  sigma_CTmax = round(fit_re$estimates$estimate[fit_re$estimates$parameter == "sigma_CTmax"], 2))
#>      no_pooling_sd partial_pooling_sd        sigma_CTmax 
#>               2.32               1.48               1.48

Here the spread of the colony effects is pulled in from the no-pooling fit (a separate CTmax per colony) to the shrunk random-effect BLUPs — the Gaussian prior b_g ~ N(0, sigma) doing exactly what a prior does. The pull is strongest for data-poor colonies; with many assays per colony there is little to shrink.

Penalised likelihood in general. More broadly, adding a penalty to the log- likelihood is maximum a posteriori (MAP) estimation under a prior: ridge is a Gaussian prior (Hoerl & Kennard, 1970), the lasso is a Laplace prior (Tibshirani, 1996; Park & Casella, 2008), and Firth’s bias-reduction penalty is the Jeffreys prior (Firth, 1993). The continuum from maximum likelihood through profile likelihood to penalised likelihood is laid out by Cole, Chu & Greenland (2014), and the penalty-as-prior equivalence — penalised likelihood equals maximum a posteriori (MAP) estimation under a prior — is standard (Bishop, 2006; Hastie, Tibshirani & Friedman, 2009).

So the real distinction is not prior versus no prior. It is how the regularisation is expressed and propagated:

freqTLS keeps its regularisation light and transparent — a documented clamp, an optional and clearly labelled random effect — and reports confidence intervals rather than a posterior. But it is honest to say it is lightly regularised, not assumption-free.

Practical guidance: when to reach for which

References