bayesTLS R Package

bayesTLS is an R package for fitting joint Bayesian four-parameter logistic (4PL) models to thermal death time (TDT) and thermal load sensitivity (TLS) data. It works with raw counts or continuous proportions, extracts classical quantities such as z, CTmax, Tcrit, and LT curves as posterior summaries, and provides tools for uncertainty propagation, group contrasts, heat-injury prediction, and bundled case-study data.

Citation

Please cite the companion paper when using bayesTLS:

Noble, D.W.A., Arnold, P.A., Nakagawa, S. & Pottier, P. (2026). A flexible modelling framework for estimating thermal tolerance and sensitivity. bioRxiv. https://doi.org/10.64898/2026.07.16.738378

This preprint is the primary citation for both the statistical framework and the bayesTLS package. You can also run citation("bayesTLS") in R for the formatted reference and a BibTeX entry.

Installing the bayesTLS package

The analytical workflow ships as an installable R package, available on CRAN. Install the released version with:

install.packages("bayesTLS")
library(bayesTLS)

To install the development version from GitHub instead:

# install.packages("remotes")  # if needed
remotes::install_github("daniel1noble/bayesTLS")

The package depends on brms (which needs a Stan backend; cmdstanr is recommended). All other imported dependencies (dplyr, ggplot2, MASS, patchwork, posterior, tibble) are CRAN packages and resolve automatically. The classical two-stage beta-binomial Stage 1 additionally uses the suggested glmmTMB package.

Key functions in bayesTLS

The coding workflow is decoupled — each step is a standalone function so users can stop, swap, or skip at any stage:

Function	Role
standardize_data()	Rename user columns to standard names; supports count responses (survived/dead/total) and continuous proportions; records response type + metadata.
make_4pl_priors() / make_4pl_formula()	Disjoint-bounds reparameterised 4PL; asymptote bounds adjustable for sublethal/PSII ranges.
fit_4pl()	Joint Bayesian 4PL via brms. Family resolved from the response — beta_binomial(identity) for counts, Beta(identity) for proportions. Centred temperature; temp_effects selects which 4PL parameters depend on temperature (all four by default; "mid" for sparse designs). Two parameterisations: midpoint (default) and direct CTmax/z (supply ctmax =/z = formulas). To compare moderator groups (species, life stage, …) in one joint fit, use by = "<moderator>" (midpoint — fits ~ temp_c * moderator on all four sub-parameters) or ctmax = ~ 0 + <moderator> (direct); the moderator is recorded so downstream readers auto-group. mid = gives explicit control of the midpoint sub-model. Random effects with a known covariance matrix — phylogenetic models or SNP-relatedness animal models — are supported via gr(group, cov = A) terms plus data2 = (validated at fit time); prior = replaces the default priors wholesale.
tls() / tls_z() / tls_ctmax() / tls_tcrit()	One call to derive z, CTmax, and (lethal data only, lethal = TRUE) Tcrit per moderator group (by =; sex, life stage, clone, …) from any fitted 4PL — a bayes_tls workflow or a hand-written brms 4PL — all from one consistent posterior. z and CTmax share one posterior subsample; seed = makes the subsample + Tcrit draws reproducible. target_surv picks the relative (default) or absolute threshold.
derive_z(), derive_tdt_curve(), derive_temperature_for_duration(), derive_tdt_landscape()	Exported TDT primitives (the building blocks tls() composes, and available directly for advanced use such as the per-temperature local z(T) profile from derive_z()); all support per-group results via by =.
extract_4pl_pars(), tdt_parameter_table(), tdt_quantile(), summarise_observed_survival()	Per-draw 4PL parameters, natural-scale posterior parameter tables (group-aware via by =), TDT-friendly quantile summaries, and observed-survival cell summaries (mean ± SE per temp × duration).
predict_survival_curves(), predict_heat_injury(), repair_rate_schoolfield()	Posterior survival curves and heat-injury accumulation under a temperature trace (both per-group via by =), with an optional Sharpe–Schoolfield repair kernel supplied by repair_rate_schoolfield().
make_temperature_scenarios(), planted_dose_from_trace()	Reference temperature traces and the analytical HI integral for validation.
ts_stage1(), ts_stage2(), ts_ci(), ts_curve()	The classical two-stage TDT pipeline for side-by-side comparison: per-temperature dose-response (stats::glm binomial or glmmTMB beta-binomial) → OLS of log10(LT50) on temperature → z, CTmax, Tcrit, with delta-method (Normal/t) or MVN intervals. One implementation, shared by the manuscript tutorial, case studies, and the bias simulation.
diagnose_tdt_fit(), bayes_R2_tls()	Sampling diagnostics (R-hat, ESS, divergences) and Bayesian R^2 summaries for a fitted workflow.
plot_survival_curves(), plot_tdt_curve(), plot_tdt_landscape(), plot_heat_injury(), plot_temperature_scenarios(), plot_temperature_density(), plot_repair_rate(), theme_tdt()	Plotting helpers for survival curves, TDT curves, tolerance landscapes, heat injury, temperature traces/densities, and repair TPCs, with a shared project theme. Where a fitted curve/surface is overlaid on data (plot_survival_curves(), plot_tdt_landscape(), plot_tdt_curve()), it is by default capped at each group's longest observed exposure (shorter durations still shown) and, for grouped fits, faceted with free per-group axes (clip_to_observed, facet_scales); pass clip_to_observed = FALSE for the full prediction grid.
get_brmsfit(), get_tls_est(), get_4pl_est(), get_surv_draws(), get_hi_draws(), has_fit(), …	Accessors for the underlying fit, posterior draws, and posterior summaries. get_tls_est() pulls draws or summary for any/all TLS quantities from a tls() object, merged on .draw so quantities keep their joint pairing (for contrasts and bivariate plots); get_4pl_est() is its 4PL-shape counterpart (draws or summary of low/up/k/mid from a fitted workflow).
clock_to_minutes(), format_interval()	Utilities: parse clock-time strings to elapsed minutes; format a posterior median with its interval as median [lower, upper].

Full reference: ?fit_4pl, ?tls, ?predict_heat_injury, etc.

bayesTLS Data

Four publicly available datasets are included with the bayesTLS package to make it easier to reproduce analyses and results in the paper and for testing purposes. These datasets are:

Dataset	Help	Endpoint
snowgum_psii	?snowgum_psii	Snowgum leaf PSII (continuous proportion)
dsuzukii	?dsuzukii	Drosophila suzukii multi-trait TDT (lethal, knockdown, fertility; per individual)
zebrafish_o2	?zebrafish_o2	Zebrafish lethal TDT across an oxygen gradient (survival counts)
aphid_tdt	?aphid_tdt	Cereal-aphid lethal TDT — three species × ages, heat + cold branches (survival counts)

Each dataset can be loaded easily using data(aphid_tdt) (as an example of loading the aphid_tdt dataset). If you want other datasets loaded then simply replace aphid_tdt. If you want to learn more about a dataset you can explore it's helpful ?aphid_tdt.

Planned features

Feature requests we are actively tracking (see the linked issues/plans for detail):

• Time-to-event / repeated-measures designs (target 1.2.0): individuals monitored through time (interval checks or continuous biosensors) rather than independent temperature × duration batches. Needs: standardize_data() to accept id/last_alive/first_dead (or exact death times) and build at-risk interval rows, plus a conditional-interval Bernoulli response branch in fit_4pl() — same 4PL parameters, so tls() (z, CTmax, relative and absolute) already works. Plan: notes/time_to_event/2026-09-26-tte-native-support-plan.qmd. (Thanks to Drew Villeneuve for bringing up this option)

Shipped:

• Covariance-matrix random effects (#12, shipped in 1.1.0 — development version; on CRAN from the next release): ctmax = ~ 1 + (1 | gr(species, cov = A)) with data2 = list(A = A) fits phylogenetic models — and the same machinery serves quantitative-genetic animal models with a SNP relatedness matrix in place of the phylogenetic correlation (gr(id, cov = A)). When by names a grouping variable that enters only through such a random effect, tls(by = ) includes its group effects automatically (override with re_formula). Design-feasibility simulations for the animal-model use (how many families, which assay designs can estimate genetic variance in CTmax and z) live in notes/quantgen_sim/.

Reproducing the associated paper

This repository also contains the code needed to reproduce the companion paper, its supplement, and the simulation results. To do this, install the package from this repository, render the Quarto documents through the Makefile, and use the tests to verify the package functions used by the manuscript and supplement. All code chunks are provided in ms/ms.qmd and ms/supplement.qmd; use make supp to render the supplement and scripts/simulations/run_simulations.R to rerun the full two-stage-bias simulation suite. The large cached artifacts the renders read (simulation outputs and brms model fits) are not stored in the repository; make data downloads them from the project's public OSF deposit (see Fetching cached outputs below) so you need not re-fit or re-simulate. Run all scenarios with Rscript scripts/simulations/run_simulations.R, or a subset by passing scenario labels, e.g. Rscript scripts/simulations/run_simulations.R scen9_tmax_060. Each scenario and the per-simulation pipeline are visible in that one script; the reusable building blocks live in scripts/simulations/sim_functions.R.

Repository layout

Use this map to find the files and outputs involved in reproduction:

• Makefile — main entry point for rendering the manuscript and supplement.
• ms/ms.qmd — manuscript source file; renders to _output/ms.html, _output/ms.docx, and _output/ms.pdf.
• ms/supplement.qmd — supplement source file; renders to _output/supp.html, _output/supp.docx, and _output/supp.pdf.
• R/ — package source for the analysis functions used by the paper; see Key functions.
• data/ — bundled package datasets as .rda files.
• inst/extdata/ — source CSV files used to build the bundled datasets.
• data-raw/make_datasets.R — script that rebuilds the .rda datasets from the source CSV files.
• output/models/ — cached brms fits (.rds) created and reused during renders.
• scripts/simulations/ — the two-stage-bias simulation: run_simulations.R (scenarios + runner) and sim_functions.R (reusable building blocks).
• scripts/ — other support scripts (presentation figures, document cleanup).
• tests/testthat/ — fast unit tests plus optional brms integration tests.
• bib/ — bibliography and the Ecology Letters CSL file.
• notes/ — dated development notes for derivations, literature checks, and simulations, organised into topic subfolders (literature/, theory/, twostage_sim/, package_design/, case_studies/, manuscript/, time_to_event/, quantgen_sim/); simulation threads keep their R scripts in an R/ subfolder and their outputs in output/.
• man/ — roxygen-generated package help files.
• _output/ — final rendered manuscript and supplement files; created by the render commands below.

Reproducing the analysis

1. Install the package from this repository so the renders use the same code you are reproducing:

   # install.packages("devtools")  # if needed
   devtools::install()

2. Fetch the cached artifacts the renders read (simulation outputs and brms model fits), which are not stored in the repository:

   make data        # results + model fits (~85 MB) -- enough to render the paper
   make data-full   # ALSO the 2.6 GB raw/draws tier (only to re-aggregate sims)

   make data is safe to re-run: a tier already present locally is skipped. No credentials are needed (the OSF deposit is public).

3. Render the manuscript, supplement, or both from the repository root:

   make all      # manuscript and supplement, all formats
   make ms       # manuscript only: _output/ms.{html,docx,pdf}
   make supp     # supplement only: _output/supp.{html,docx,pdf}

   You can also render one format at a time, for example make ms-pdf or make supp-html.

   The first supplement render may take several minutes because it fits or checks cached models. Subsequent renders reuse the .rds files in output/models/ and should be much faster.

4. Inspect the rendered files in _output/. Use make clean to remove only _output/, or make build-clean to also remove the out-of-tree Quarto build cache at ~/Library/Caches/tls-render/.

Fetching cached outputs from OSF

The simulation outputs and brms model fits the manuscript reads are large and are not committed to the repository. They live in the project's public OSF deposit (node c6dxy) as tiered zip archives, fetched with make data:

Command	Downloads	Size	Use
make data	results.zip + models.zip	~85 MB	render the paper (default)
make data-full	the above plus sim_full.zip	+2.6 GB	re-aggregate the raw simulation

make data is idempotent — a tier already present under output/ is skipped (pass TIER=--force to re-download), each download is md5-verified against data-raw/osf_manifest.json, and the archives unpack straight into output/. No credentials are required. Maintainers refresh the deposit with make osf-push (needs an OSF_PAT token); see scripts/osf_publish.R.

4. Run the package tests if you want to verify the analysis functions separately from the manuscript render:

   devtools::test()
   
   Sys.setenv(RUN_BRMS_TESTS = "true")
   devtools::test()

   The first command runs the fast unit tests. The second also runs the gated brms integration tests, which fit small cached models and check parameter recovery.

Rendering details

The manuscript and supplement are independent Quarto documents. There is no project-level _quarto.yml; each .qmd file carries its own YAML so the files can be rendered separately. The Makefile renders outside the Dropbox-synced project directory, then copies only the final HTML, DOCX, and PDF files back to _output/. This avoids conflicted intermediate files while preserving cached model fits in output/models/.

Within a document, use Quarto's @ references as normal (@fig-X, @tbl-X, @eq-X, @sec-X). Between the manuscript and supplement, use plain text, such as "Equation 7 of the manuscript" in the supplement or "Figure S2" in the manuscript. The supplement render automatically labels figures, tables, and equations with an "S" prefix.

See CLAUDE.md §8a for more detail on the render architecture, cross-reference handling, and the brms cache.

Case studies in the supplement

The supplement first walks through the full workflow with a simulated dataset (including the classical two-stage pipeline for comparison), then repeats it on four empirical case studies — one per bundled dataset (see bayesTLS Data) — each showcasing a different modelling need:

1. Zebrafish × oxygen (zebrafish_o2) — lethal TDT with oxygen treatment as a covariate on CTmax and z; per-treatment thermal limits.
2. Cereal aphids (aphid_tdt) — three species compared in one joint model (species on CTmax and z); per-species limits and heat injury under a field temperature trace.
3. Snow gum leaf PSII (snowgum_psii) — continuous-proportion response fit with a Beta likelihood; light vs dark recovery comparison.
4. Drosophila suzukii (dsuzukii) — per-individual multi-trait TDT: lethal endpoint compared across sexes, then sublethal knockdown and fertility endpoints.

The supplement also contains the extended two-stage-bias simulation results and sensitivity tests, the code behind the manuscript's cross-case-study figures, a frequentist companion analysis (freqTLS), and the mathematical derivations for the threshold-summary correction factors.