Package 'bayesTLS'

Description

Survival of three cereal-aphid species across a broad range of stressful high and low temperatures, the model-ready frame for the multi-species case study. Aphids of three ages (2, 6, 12 days old) were exposed to a heat branch (34–40 degrees C) or a cold branch (-11 to -3 degrees C) for a range of durations and scored alive/dead after recovery. One row per assay group; branch flags the heat vs cold series. Subset to one branch (and typically one age) and fit CTmax and z as functions of species in one joint 4PL to compare species with full posterior uncertainty.

Usage

aphid_tdt

Format

A data frame with 3041 rows and 7 variables:

species

Species, a factor with levels M_dirhodum (Metopolophium dirhodum), S_avenae (Sitobion avenae), R_padi (Rhopalosiphum padi).

age

Age in days, a factor with levels 2, 6, 12.

branch

Stress branch, a factor with levels heat (34–40 degrees C), cold (-11 to -3 degrees C).

temp

Assay temperature (degrees C).

duration_min

Exposure duration (minutes).

n_total

Number of aphids treated.

n_surv

Number surviving after treatment and recovery.

Source

Li Y-J, Chen S-Y, Jørgensen LB, Overgaard J, Renault D, Colinet H, Ma C-S (2023). Data for: Interspecific differences in thermal tolerance landscape explain aphid community abundance under climate change. Dryad, doi:10.5061/dryad.mcvdnck4j (Dryad CC0). Associated article: Journal of Thermal Biology 114: 103583, doi:10.1016/j.jtherbio.2023.103583. Raw file: system.file("extdata", "data_lethal_TDT_aphid.csv", package = "bayesTLS").

Examples

## Not run: 
a <- subset(aphid_tdt, branch == "heat" & age == "6")
std <- standardize_data(a, temp = "temp", duration = "duration_min",
                        n_total = "n_total", n_surv = "n_surv",
                        duration_unit = "minutes")
wf <- fit_4pl(std, ctmax = ~ 0 + species, z = ~ 0 + species, t_ref = 60)
tls(wf, by = "species", lethal = TRUE)   # z, CTmax, T_crit per species

## End(Not run)

Description

A thin tidy wrapper around brms::bayes_R2() for a fitted bayes_tls workflow: it pulls the underlying brmsfit (erroring clearly on an unfitted workflow, via get_brmsfit()) and returns the posterior median $R^2$ with its standard error and 95% credible interval as a one-row tibble, rather than the bare matrix brms::bayes_R2() returns. Pass ... straight through to brms::bayes_R2() (e.g. re_formula = NA to exclude group-level effects, or ⁠ndraws =⁠ to subsample). Map over a named list of workflows to build a multi-fit table.

Usage

bayes_R2_tls(workflow, ...)

Arguments

Value

A one-row tibble with columns estimate, est_error, lower, and upper (the lower and upper credible bounds, 2.5% and 97.5% by default; controlled by probs passed via ...).

See Also

brms::bayes_R2(), diagnose_tdt_fit()

Examples

## Not run: 
wf <- fit_4pl(std)
bayes_R2_tls(wf)
# Multiple fits in one table:
purrr::imap_dfr(list(binom = wf1, beta = wf2),
                ~ cbind(model = .y, bayes_R2_tls(.x)))

## End(Not run)

Description

Accepts POSIXt, hms / difftime, numeric fractions of a day (Excel time) or bare numeric minutes, and character strings: "HH:MM:SS", "HH:MM", bare numeric strings (minutes), and durations beyond 24 h (e.g. "25:30:00"). Character strings are parsed element-wise; malformed entries become NA.

Usage

clock_to_minutes(x)

Arguments

Value

Numeric vector of minutes.

Examples

clock_to_minutes("08:30:00")
clock_to_minutes("25:30")   # 25 h 30 min = 1530 min
clock_to_minutes(0.5)       # half a day = 720 min

Description

Returns the per-draw duration at which population-level survival crosses the chosen threshold, at each temperature in temp_grid.

Usage

derive_tdt_curve(
  workflow,
  temp_grid,
  duration_grid = NULL,
  target_surv = "relative",
  ndraws = NULL,
  probs = c(0.025, 0.5, 0.975),
  time_multiplier = NULL,
  output_time_unit = "min",
  by = NULL
)

Arguments

Details

Two threshold modes are supported via target_surv:

• "relative" (default): the duration at which survival reaches the midpoint between the fitted lower and upper asymptotes, i.e. (low + up)/2. This is the 4PL mid parameter on the natural time axis, returned directly from posterior_linpred(nlpar = "mid") — no numerical inversion. When low is about 0 and up about 1 it coincides with the classical LT50.

• "absolute" (or a numeric p in ⁠(0, 1)⁠): the duration at which survival crosses the literal probability p (0.5 by default). The inversion is numerical — predict survival on a dense duration grid, then approx() through p.

This is the horizontal read of the survival surface: fix a survival threshold, read off the time required to reach it at each temperature.

Value

A list with draws (per-draw threshold durations; target_surv column is a character label), summary (quantile summary by temperature, plus moderator columns for grouped fits), target_surv (the label), time_multiplier, output_time_unit, and by.

Examples

## Not run: 
wf  <- fit_4pl(std)
crv <- derive_tdt_curve(wf, temp_grid = c(38, 40, 42))  # relative LT curve
crv$summary

## End(Not run)

Description

Returns the posterior median (and 95% credible band) of survival probability on a 2-D grid covering the experimental range. Suitable for plotting as a heatmap or contour with plot_tdt_landscape() (Phase 1d). For curves at a handful of temperatures, use predict_survival_curves() instead.

Usage

derive_tdt_landscape(
  workflow,
  temp_grid = NULL,
  duration_grid = NULL,
  ndraws = NULL,
  probs = c(0.025, 0.5, 0.975),
  by = NULL
)

Arguments

Value

A list with the same shape as predict_survival_curves() (a summary tibble with temp, duration, survival_median, survival_lower, survival_upper, plus any moderator columns for grouped fits, draws_matrix and grid).

Examples

## Not run: 
wf  <- fit_4pl(d, ...)
lsp <- derive_tdt_landscape(wf)
lsp$summary

## End(Not run)

Description

The vertical read of the survival surface: fix an exposure duration, find the temperature at which the posterior survival reaches the chosen threshold. Returns one temperature per posterior draw.

Usage

derive_temperature_for_duration(
  workflow,
  exposure_duration,
  temp_grid,
  target_surv = "relative",
  ndraws = NULL,
  probs = c(0.025, 0.5, 0.975),
  seed = NULL,
  by = NULL
)

Arguments

Details

Threshold modes (via target_surv) match derive_tdt_curve():

• "relative" (default) → temperature at which mid(T) = log10(exposure_duration) per draw. The inversion is done analytically per draw: extract posterior_linpred(nlpar = "mid") over temp_grid, then approx() to the target log10-time.

• "absolute" (= 0.5) or numeric p in ⁠(0, 1)⁠ → existing numerical inversion of the 4PL survival surface at the literal probability p.

This is the primitive used by tls() to derive CTmax at t_ref.

Value

A list with draws (per-draw threshold temperatures; target_surv column is a character label), summary (quantile summary), exposure_duration, target_surv (the label), target_mode, target_prob. A grouped fit adds the moderator column(s).

Examples

## Not run: 
wf <- fit_4pl(std)
# Temperature giving the relative midpoint threshold after a 60-unit exposure:
tt <- derive_temperature_for_duration(wf, exposure_duration = 60,
                                       temp_grid = seq(36, 44, by = 0.1))
tt$summary

## End(Not run)

Description

Derives $z = -1 / (\mathrm{d}/\mathrm{d}T\,\log_{10}\mathrm{LT}(T))$ per posterior draw, read straight from the fitted 4PL coefficients — no regression. There are two regimes:

Usage

derive_z(
  workflow,
  target_surv = "relative",
  temp_grid = NULL,
  ndraws = NULL,
  probs = c(0.025, 0.5, 0.975),
  h = 0.001,
  seed = NULL,
  by = NULL
)

Arguments

Details

• Relative threshold (default; the $(\ell+u)/2$ midpoint): $\log_{10}\mathrm{LT}_{\text{rel}}(T) = \mathrm{mid}(T) = \beta_0 + \beta_1 (T-\bar T)$ is exactly linear, so $z = -1/\beta_1$ where $\beta_1$ is the temperature slope on mid (b_mid_temp_c). The asymptotes $\ell, u$ and slope $k$ do not enter (the midpoint cancels the curve asymmetry). z is constant in temperature.

• Absolute threshold $p$: the LT curve gains the asymmetry-correction term, $\log_{10}\mathrm{LT}_p(T) = \mathrm{mid}(T) + \tfrac{1}{k(T)}\log\tfrac{u(T)-p}{p-\ell(T)}$. When $\ell$, $u$ or $k$ carry temperature effects this bends the curve, so z varies with temperature. A local $z(T) = -1/m(T)$ is computed at each assay temperature, where the local slope $m(T)$ is obtained by a central finite difference of the closed-form LT curve (step h). When the shape parameters are constant in T the correction is flat and this reduces to $-1/\beta_1$.

The returned pooled z (the default single summary) is the per-draw mean of the local $z(T)$ over temp_grid. The full per-temperature local $z(T)$ is also returned. z is invariant to the time unit (a constant time-multiplier shifts the LT intercept, not its slope), so no time_multiplier is needed here.

Value

A list with:

• draws: tibble ⁠(.draw, z)⁠ — pooled per-draw z.

• summary: tibble ⁠(z_median, z_lower, z_upper)⁠.

• local_draws: tibble ⁠(.draw, temp, z)⁠ — local z(T) per draw.

• local_summary: tibble ⁠(temp, z_median, z_lower, z_upper)⁠.

• target_surv, temp_grid. For a grouped fit each tibble gains the moderator column(s).

Examples

## Not run: 
wf <- fit_4pl(std)
z  <- derive_z(wf)             # relative: z = -1 / b_mid_temp_c per draw
z$summary
derive_z(wf, target_surv = "absolute")$local_summary  # local z(T)

## End(Not run)

Description

Returns a tibble with the per-fit summary numbers a reviewer will want at a glance: max Rhat, minimum bulk and tail ESS, divergent transitions, tree-depth saturations, minimum BFMI per chain, and pass flags for each criterion plus a combined all_pass. Healthy values:

• Rhat < 1.01

• ESS bulk and tail > 400

• zero divergent transitions

• zero tree-depth saturations

• BFMI > 0.3 per chain

Usage

diagnose_tdt_fit(workflow)

Arguments

Details

The treedepth_max field is the model's max_treedepth setting (the ceiling passed to brms::brm(), default 10); a saturation is a post-warmup transition that hit that ceiling. A fit that merely tops out below the ceiling has zero saturations. Per-chain BFMI is computed from the energy diagnostic in brms::nuts_params(.) following the standard Stan definition (Var(Delta E)/Var(E)).

Value

A tibble with one row of diagnostic statistics.

Examples

## Not run: 
diagnose_tdt_fit(wf)

## End(Not run)

Description

Per-individual thermal-tolerance assays for the spotted-wing fly (Drosophila suzukii), one row per fly, carrying three thermal-tolerance endpoints measured under static heat exposures at 34–38 degrees C: a lethal endpoint (dead), a sublethal knockdown time-to-event (t_coma), and a sublethal reproductive endpoint (prod). The model-ready frame for Case Study 4; aggregate dead to counts for the beta-binomial lethal fit, or use t_coma / prod directly for the sublethal endpoints. lvl indexes the exposure-duration grid as a percentage of the estimated median time-to-coma from the authors' initial TDT curves; time is the realised duration in minutes.

Usage

dsuzukii

Format

A data frame with 1407 rows and 9 variables:

id

Unique individual identifier (temp-lvl-sex-rep).

temp

Assay temperature (degrees C).

lvl

Exposure duration as a percentage of the estimated median time-to-coma from the authors' initial TDT curves.

time

Exposure duration (minutes).

sex

Sex, a factor with levels F, M.

rep

Replicate vial within a temperature x lvl x sex cell.

prod

Reproductive productivity (offspring per female per day).

dead

Mortality indicator: 1 = died, 0 = survived.

t_coma

Time to heat coma (minutes); NA where no coma was recorded for that individual.

Source

Ørsted M, Willot Q, Olsen AK, Kongsgaard V, Overgaard J (2024). Data for: Thermal limits of survival and reproduction depend on stress duration: a case study of Drosophila suzukii. Zenodo, doi:10.5281/zenodo.10602268 (distributed under CC BY 4.0). Associated article: doi:10.1111/ele.14421. Raw file: system.file("extdata", "data_multitrait_TDT_drosophila_suzukii.csv", package = "bayesTLS").

Examples

# Lethal endpoint: aggregate per-individual deaths to cell counts, then
# standardise for the beta-binomial 4PL.
cells <- dplyr::summarise(
  dplyr::group_by(dsuzukii, temp, time, sex),
  n_total = dplyr::n(), n_dead = sum(dead), .groups = "drop")
std <- standardize_data(cells, temp = "temp", duration = "time",
                        n_total = "n_total", n_dead = "n_dead",
                        random_effects = "sex", duration_unit = "minutes")

Description

Faithfully extracts ALL FOUR natural-scale 4PL sub-parameters (low, up, k, mid) per posterior draw, at one or more assay temperatures and per moderator group, retaining whatever temperature and/or group structure each sub-parameter was fit with (~ 1, ~ temp_c, ~ temp_c * group, random effects, …). It surfaces the shared TDT engine (tls_eval_subpars() over a tls_build_grid() temperature × moderator grid via brms::posterior_linpred(nlpar=)), so the result is parameterisation- and coding-agnostic — no coefficient-name parsing.

Usage

extract_4pl_pars(
  workflow,
  temps = NULL,
  by = NULL,
  re_formula = NA,
  ndraws = NULL
)

Arguments

Details

mid is the natural-scale midpoint sub-parameter (the log10 time at the per-draw relative threshold) evaluated at each temperature; on an absolute fit it is still the bare mid nlpar (the asymmetry correction is applied downstream by the threshold-specific readers, not folded in here).

This replaces the earlier constant-shape extractor that discarded the temperature slopes on low/up/k. The constant-in-T view the heat-injury integral needs (⁠low, up, k⁠ at the centring temperature plus a linear mid_int/mid_temp) now lives in the internal helper hi_pars(), which predict_heat_injury() calls for shape = "constant".

Value

A long tibble with ⁠(.draw, temp, low, up, k, mid)⁠ columns (plus the moderator column(s) for a grouped fit, immediately after .draw), one row per draw × temperature × group, filtered to rows producing valid parameter values (k > 0, up > low, all finite). temp is the uncentred assay temperature (°C).

See Also

hi_pars() for the constant-shape view used by the heat-injury integral; get_4pl_est() which wraps this for the "draws" view.

Examples

## Not run: 
pars <- extract_4pl_pars(wf)                 # observed temps, fit's groups
extract_4pl_pars(wf, temps = c(30, 33, 36))  # at chosen temperatures
head(pars)

## End(Not run)

Description

Wraps make_4pl_formula() + make_4pl_priors() + brms::brm(). Returns a workflow object containing the fit, data, formula, prior, and metadata that downstream helpers (e.g. tls()) read from.

Usage

fit_4pl(
  data,
  random_effects = NULL,
  temp_effects = c("low", "up", "k", "mid"),
  lower = 0,
  upper = 1,
  family = NULL,
  prior = NULL,
  chains = 4,
  iter = 4000,
  warmup = floor(iter/2),
  cores = chains,
  seed = 123,
  backend = "cmdstanr",
  control = list(adapt_delta = 0.95, max_treedepth = 12),
  init = 0,
  file = NULL,
  file_refit = "on_change",
  fit = TRUE,
  bounds = NULL,
  ctmax = NULL,
  z = NULL,
  up = NULL,
  low = NULL,
  k = NULL,
  mid = NULL,
  by = NULL,
  threshold = c("relative", "absolute"),
  p = 0.5,
  t_ref = 60,
  ...
)

Arguments

Details

By default all four 4PL sub-parameters get a temp_c slope, and random intercepts attach to mid only. Use temp_effects = "mid" for the classical constant-shape TDT model (recommended for sparse designs — see temp_effects). To fit separate categories (life stages, species, populations), filter data per category and call this function once per subset.

Value

A list of class "bayes_tls" with elements fit, data, formula, prior, meta. See print.bayes_tls(), summary.bayes_tls(), and plot.bayes_tls() for the available methods.

Examples

## Not run: 
d  <- standardize_data(raw, ...)
wf <- fit_4pl(d, fit = FALSE)         # inspect spec only
wf <- fit_4pl(d, file = "output/models/my_fit")
wf <- fit_4pl(d, lower = 0.85)        # sublethal data

## End(Not run)

Description

Format a posterior median plus credible interval as a single string

Usage

format_interval(median, lower, upper, digits = 2)

Arguments

Value

Character like "5.12 [4.87, 5.4]". A non-finite median yields NA_character_ (rather than "NA [...]"); a non-finite bound is shown as an en dash, so the strings stay table-ready.

Examples

format_interval(5.123, 4.872, 5.401)
format_interval(NA, 1, 2)   # -> NA

Description

The 4PL-shape companion to get_tls_est(). Returns the posterior of the natural-scale 4PL curve parameters from a fitted fit_4pl() workflow, as either per-draw values or a median + 95% credible-interval summary, per moderator group. Whereas get_tls_est() operates on a tls() object and returns the derived TLS quantities (z, CTmax, Tcrit), this operates on the fitted workflow and returns the raw curve parameters those quantities are derived from. A thin wrapper: "summary" calls tdt_parameter_table(), "draws" calls extract_4pl_pars().

Usage

get_4pl_est(workflow, what = c("summary", "draws"), by = NULL, temps = NULL)

Arguments

Value

A tibble (plus the moderator column(s) for a grouped fit): for "summary", one row per parameter with median and credible bounds at T_bar; for "draws", one row per posterior draw × temperature with temp, low, up, k, mid (and .draw).

See Also

get_tls_est() for the derived TLS quantities; tdt_parameter_table() and extract_4pl_pars(), which this wraps.

Examples

## Not run: 
wf <- fit_4pl(standardize_data(my_data, temp = "temp", duration = "time",
                               n_total = "n", n_dead = "dead"))
get_4pl_est(wf, "summary")                 # low/up/k/mid median + 95% CrI
get_4pl_est(wf, "draws", by = "species")   # per-draw params, per group

## End(Not run)

Description

Returns the underlying brms::brmsfit object held in a bayes_tls workflow, ready for any brms / posterior / bayesplot helper (e.g. brms::bayes_R2(), brms::mcmc_plot(), posterior::as_draws_df()). Prefer this over reaching into workflow$fit directly: it errors clearly on an unfitted workflow instead of silently returning NULL. Pairs with the predicate has_fit().

Usage

get_brmsfit(workflow)

Arguments

Value

The brmsfit object stored in the workflow.

Examples

## Not run: 
wf  <- fit_4pl(std)
fit <- get_brmsfit(wf)
brms::bayes_R2(fit)

## End(Not run)

Description

Long-format tibble of per-draw HI, dose, survival and mortality at every time step of the supplied temperature trace. Requires that predict_heat_injury() was called with save_draws = TRUE.

Usage

get_hi_draws(hi)

Arguments

Value

A tibble with columns .draw, time, temp, dose, hi, survival, mortality.

Examples

## Not run: 
hi <- predict_heat_injury(trace, wf, save_draws = TRUE)
get_hi_draws(hi)

## End(Not run)

Description

Accepts either a predict_survival_curves() result (static survival on a temperature × duration grid) or a predict_heat_injury() result (dynamic survival along a temperature trace). In the heat-injury case requires save_draws = TRUE.

Usage

get_surv_draws(x)

Arguments

Value

A long-format tibble. For predict_survival_curves(): .draw, temp, duration, survival. For predict_heat_injury(): .draw, time, temp, survival. For a grouped fit the moderator column(s) are also carried (so per-group draws are not collapsed).

Examples

## Not run: 
psc <- predict_survival_curves(wf, temps = c(32, 34),
                               durations = c(0.5, 1))
get_surv_draws(psc)

hi <- predict_heat_injury(trace, wf, save_draws = TRUE)
get_surv_draws(hi)

## End(Not run)

Description

Pulls the posterior draws or the summary (median + 95% credible interval) for any or all thermal-load-sensitivity quantities (z, CTmax, and, when the fit was extracted with lethal = TRUE, Tcrit) from a tls() result. Draws are returned in tidy-long form (quantity, .draw, value, plus any grouping column), so a posterior contrast is a join on .draw followed by a difference of value.

Usage

get_tls_est(x, what = c("summary", "draws"), params = NULL)

Arguments

Value

A tibble. For what = "summary": columns quantity, median, lower, upper (plus grouping columns for a grouped fit). For what = "draws": columns quantity, .draw, value (plus grouping columns).

See Also

tls(), which produces the object this reads.

Examples

## Not run: 
fit <- fit_4pl(standardize_data(my_data, temp = "temp", duration = "time",
                                n_total = "n", n_dead = "dead"))
est <- tls(fit, lethal = TRUE)
get_tls_est(est, "summary")             # all quantities: median + 95% CrI
zd  <- get_tls_est(est, "draws", "z")   # z draws, ready for a contrast

## End(Not run)

Description

Whether a workflow has been fitted

Usage

has_fit(workflow)

Arguments

Value

Logical scalar.

Examples

wf <- list(fit = NULL); class(wf) <- "bayes_tls"
has_fit(wf)

Description

Constructs the joint 4PL with the disjoint-bounds reparameterisation on the asymptotes (low, up) and a positivity reparam on k. All four 4PL sub-parameters (lowraw, upraw, logk, mid) get a temp_c slope — the model assumes temperature can affect every aspect of the dose-response curve. Random intercepts (if any) attach to mid.

Usage

make_4pl_formula(
  random_effects = NULL,
  bounds = NULL,
  family = brms::beta_binomial(link = "identity"),
  response_var = "survival",
  temp_effects = c("low", "up", "k", "mid"),
  ctmax = NULL,
  z = NULL,
  up = NULL,
  low = NULL,
  k = NULL,
  mid = NULL,
  by = NULL,
  threshold = c("relative", "absolute"),
  p = 0.5,
  log10_tref = 0,
  lower = 0,
  upper = 1
)

Arguments

Details

The response term depends on the family:

• Count families (binomial, beta_binomial): the response is ⁠n_surv | trials(n_total) ~ <4PL>⁠.

• Beta family (continuous proportion): the response is ⁠<response_var> ~ <4PL>⁠ with no trials() term. With link = "identity" the 4PL value is the mean directly (safe because the reparam keeps it in ⁠(0, 1)⁠); with link = "logit" the 4PL is wrapped in logit() so the mean is still the 4PL value.

If a parameter has no real temperature effect, its temp_c slope shrinks toward zero under the prior and the fit reduces cleanly to the simpler case.

Value

A brmsformula object.

Examples

make_4pl_formula()
make_4pl_formula(family = binomial(link = "identity"))
make_4pl_formula(temp_effects = "mid")            # constant-shape TDT
make_4pl_formula(family = brms::Beta(link = "identity"),
                 response_var = "survival")
make_4pl_formula(ctmax = ~ life_stage, z = ~ life_stage)  # direct CTmax/z

Description

Priors are weakly informative and do not adapt to the data (which would understate uncertainty). They are constructed to match make_4pl_formula()'s disjoint-bounds reparameterisation low = low_min + inv_logit(lowraw) * low_w and up = up_min + inv_logit(upraw) * up_w for the asymptote intervals derived from lower and upper by compute_4pl_bounds().

Usage

make_4pl_priors(
  data,
  lower = 0,
  upper = 1,
  random_effects = NULL,
  prior_random_sd = "exponential(2)",
  prior_phi = "gamma(2, 0.1)",
  bounds = NULL,
  ctmax = NULL,
  z = NULL,
  up = NULL,
  low = NULL,
  k = NULL
)

Arguments

Details

Each of the four 4PL sub-parameters (lowraw, upraw, logk, mid) gets an Intercept-specific prior plus a general class = "b" prior covering the temp_c slope. Random-effect SDs on mid get one prior_random_sd row per grouping variable.

Value

A brmsprior object ready to pass to fit_4pl().

Examples

raw <- data.frame(
  temperature_C = rep(c(30, 32, 34), each = 4),
  exposure_h    = rep(c(1, 2, 4, 8), times = 3),
  n             = 30L,
  alive         = c(29, 28, 25, 5, 30, 27, 18, 2, 28, 22, 10, 1)
)
d <- standardize_data(raw, temp = "temperature_C", duration = "exposure_h",
                      n_total = "n", n_surv = "alive")
make_4pl_priors(d)

Description

Returns a list of four tibbles, each with columns time and temp, representing the canonical validation scenarios:

Usage

make_temperature_scenarios(
  baseline = 20,
  spike_temp = 30,
  n_hours = 96,
  dt_hours = 1,
  spike_times_single = 24,
  spike_times_multi = c(24, 48, 72),
  diurnal_n_days = 7,
  diurnal_night_temp = baseline,
  diurnal_day_peaks = NULL,
  diurnal_peak_fwhm = 6,
  diurnal_noise_sd = 0.3,
  diurnal_seed = 1L
)

Arguments

Details

• flat — constant baseline temperature for n_hours hours.

• single_spike — flat baseline with one time step (dt_hours wide) at spike_temp placed at each of spike_times_single (default one spike).

• multi_spike — flat baseline with one time step (dt_hours wide) at spike_temp placed at each of spike_times_multi (default three spikes).

• diurnal — multi-day diurnal cycle with day-to-day variation in peak temperature and small hourly fluctuations on top. Useful as a stand-in for a natural thermal regime where some days exceed T_crit and accrue injury while others stay below.

Each spike is one time step (dt_hours) wide, so the analytical dose it delivers is the instantaneous rate times the step width, $100\cdot 10^{(T_{spike}-CT_{max,1hr})/z}\cdot dt_{hours}$ % LT50-dose. With the default dt_hours = 1 this reduces to $100\cdot 10^{(T_{spike}-CT_{max,1hr})/z}$, so setting spike_temp = CTmax_1hr delivers ~100% LT50-dose per spike — a calibration target for predict_heat_injury() validation. For dt_hours != 1 the dose scales by dt_hours.

Value

A named list of four tibbles (flat, single_spike, multi_spike, diurnal), each with columns time (numeric, hours from start) and temp (°C). Note the traces differ in length: flat, single_spike, and multi_spike span n_hours (default 96) while diurnal spans diurnal_n_days * 24 (default 168) hours.

Examples

scens <- make_temperature_scenarios()
lapply(scens, head)
# Diurnal calibrated for a shrimp-like organism (T_crit ~ 25 °C):
make_temperature_scenarios(
  baseline = 18,
  diurnal_n_days = 7,
  diurnal_night_temp = 19,
  diurnal_day_peaks  = c(24, 27, 24.5, 27.5, 23, 26.5, 25.5)
)$diurnal

Description

Implements the classical HI integral exactly (Equation 8 of the manuscript), given known TDT parameters. Use this as the truth that predict_heat_injury() should recover when fed the same trace and a model fit whose posterior is consistent with ⁠(z, CTmax_1hr, T_c)⁠.

Usage

planted_dose_from_trace(trace, z, CTmax_1hr, T_c)

Arguments

Details

Integrated by forward Euler from zero: the interval ending at point $i$ accrues the rate at its start (point $i-1$) times that interval's width,

$\Delta HI_i = 100 \cdot 10^{(T_{i-1} - CT_{max,1hr}) / z} \cdot (t_i - t_{i-1})$

when ⁠T_{i-1} > T_c⁠, and zero otherwise (and hi_inc[1] = 0). Cumulative HI is the running sum. This matches predict_heat_injury()'s integrator.

Value

The input trace augmented with two new columns: - hi_inc — instantaneous HI contribution at each hour (%). - hi_cumulative — running sum of hi_inc (%).

Examples

scens <- make_temperature_scenarios(spike_temp = 30)
planted_dose_from_trace(scens$single_spike, z = 5, CTmax_1hr = 32, T_c = 25)

Description

Two stacked panels: cumulative heat injury (%) on top, predicted survival fraction on the bottom, both with 95% credible bands.

Usage

plot_heat_injury(hi, lt50_threshold = 100)

Arguments

Value

A ggplot object (combined via patchwork).

Examples

## Not run: 
hi <- predict_heat_injury(scens$single_spike, wf)
plot_heat_injury(hi)

## End(Not run)

Description

Plot a Sharpe-Schoolfield repair TPC

Usage

plot_repair_rate(temp_grid, repair_pars)

Arguments

Value

A ggplot object.

Examples

rp <- list(TA = 14065, TAL = 50000, TAH = 120000,
           TL = 10.5 + 273.15, TH = 22.5 + 273.15,
           TREF = 17 + 273.15, r_ref = 0.01)
plot_repair_rate(seq(5, 30, length.out = 200), rp)

Description

Plots posterior median ± 95% CrI of survival as a function of duration, one coloured line per assay temperature. Optionally overlays observed survival proportions per replicate.

Usage

plot_survival_curves(
  pred,
  observed = NULL,
  log_time = FALSE,
  palette = "viridis",
  response_label = "Survival probability",
  clip_to_observed = TRUE,
  facet_scales = c("auto", "fixed", "free", "free_x", "free_y")
)

Arguments

Details

Colour mapping uses the viridis palette by default (yellow at the highest temperature, dark purple at the lowest); the exposure-duration axis is linear by default, with the classical log-time TDT view available via log_time = TRUE.

Value

A ggplot object.

Examples

## Not run: 
pred <- predict_survival_curves(wf, temps = c(30, 32, 34, 36))
plot_survival_curves(pred, observed = wf$data)
plot_survival_curves(pred, log_time = TRUE)  # classical TDT view

## End(Not run)

Description

Plots posterior median +/- 95% CrI of the duration to reach target_surv across temperature. By default, returns a two-panel figure: a linear-time panel (the absolute time to reach the target) alongside a log10-time panel (the classical TDT view where the relationship is near-linear with slope -1/z). Both panels share posterior draws.

Usage

plot_tdt_curve(
  ltx,
  panels = c("both", "linear", "log"),
  colour = "#146C7C",
  observed = NULL,
  clip_to_observed = TRUE
)

Arguments

Details

Grouped output from derive_tdt_curve() with by = ... is faceted automatically so each moderator level gets its own TDT curve panel.

Internal y-axis behaviour: if the 95% credible interval of the LT_x curve extends above 48 hours (converted to ltx$output_time_unit), both panels' y-axes are clamped to 48 hours so the readable part of the curve isn't compressed by an extreme upper tail at cool temperatures. If the data already stay below 48 hours, no limits are applied and ggplot2 chooses a data-driven range. Callers can override either way by appending their own + ggplot2::coord_cartesian(ylim = ...) to the returned plot.

Value

A ggplot object (single panel) or patchwork composition (two panels).

Examples

## Not run: 
ltx <- derive_tdt_curve(wf, temp_grid = seq(29, 37, by = 0.5))
plot_tdt_curve(ltx)                  # two-panel default
plot_tdt_curve(ltx, panels = "log")  # classical TDT view only

# Override the internal 48-h cap by appending a coord_cartesian:
plot_tdt_curve(ltx) +
  ggplot2::coord_cartesian(ylim = c(0, 240))

## End(Not run)

Description

Exposure-duration axis is linear by default; pass log_time = TRUE for the classical log-time TDT presentation.

Usage

plot_tdt_landscape(
  landscape,
  observed = NULL,
  contours = c(0.25, 0.5, 0.75),
  log_time = FALSE,
  clip_to_observed = TRUE,
  facet_scales = c("auto", "fixed", "free", "free_x", "free_y")
)

Arguments

Details

A grouped fit (a moderator on CTmax/z, e.g. oxygen or species) produces one heatmap panel per group, mirroring plot_survival_curves(); a single-condition fit returns the one ungrouped heatmap.

Value

A ggplot object.

Examples

## Not run: 
lsp <- derive_tdt_landscape(wf)
plot_tdt_landscape(lsp, observed = wf$data)

## End(Not run)

Description

Density plot of a per-draw temperature posterior, with a horizontal segment marking the 95% CrI and a point at the median. Useful for visualising derive_temperature_for_duration() output.

Usage

plot_temperature_density(temp_post, truth = NULL, x_label = "Temperature (°C)")

Arguments

Value

A ggplot object.

Examples

## Not run: 
tp <- derive_temperature_for_duration(wf, exposure_duration = 60,
                                       temp_grid = seq(36, 44, by = 0.1))
plot_temperature_density(tp)

## End(Not run)

Description

One panel per scenario stacked vertically, with optional horizontal dashed line at the damage threshold T_c.

Usage

plot_temperature_scenarios(scens, T_c = NULL)

Arguments

Value

A ggplot object.

Examples

scens <- make_temperature_scenarios()
plot_temperature_scenarios(scens, T_c = 24)

Description

Delegates to plot() on the underlying brmsfit (via get_brmsfit()), which produces brms's default mcmc_combo layout (per-parameter density on the left, post-warmup trace on the right, chains coloured). Use this to eyeball chain mixing alongside the numeric Rhat / ESS in summary.bayes_tls().

Usage

## S3 method for class 'bayes_tls'
plot(x, ...)

Arguments

Value

The brms plot output (invisibly): a list of bayesplot_grid grid objects, one per plotted page.

Examples

## Not run: 
plot(wf)
plot(wf, variable = "^b_mid", regex = TRUE)
plot(wf, combo = c("dens_overlay", "trace"))

## End(Not run)

Description

Propagates the model's full posterior through an Eulerian damage- accumulation integration along the supplied temperature trace, returning the posterior median and 95% credible band of:

Usage

predict_heat_injury(
  trace,
  workflow,
  target_surv = "relative",
  T_c = NULL,
  trace_unit = "hours",
  ndraws = NULL,
  repair = FALSE,
  repair_pars = NULL,
  repair_scales_with_survival = TRUE,
  irreversible_mortality = TRUE,
  save_draws = FALSE,
  seed = NULL,
  by = NULL,
  shape = c("constant", "varying")
)

Arguments

Details

• HI(t) — cumulative heat injury, in percent of an LT_target_surv dose. When HI(t) = 100, the population has accumulated one full dose at the chosen survival threshold (default 50% mortality).

• S(t) — predicted survival fraction, mapped from the cumulative dose through the fitted 4PL.

Optionally adds a temperature-dependent repair rate via repair_rate_schoolfield().

Value

A list with elements:

• summary: tibble with time, temp, and posterior median + 95% CrI for hi, survival, and mortality at each time step (plus the moderator column(s) for a grouped fit).

• draws: optional per-draw trajectories (when save_draws = TRUE).

• meta: list of inputs used.

Examples

## Not run: 
scens <- make_temperature_scenarios()
hi    <- predict_heat_injury(scens$single_spike, wf)
hi$summary

## End(Not run)

Description

Predicts survival probability at every combination of temps × durations under a fitted 4PL, then returns the posterior median and 95% credible interval at each grid point. Random effects are marginalised out.

Usage

predict_survival_curves(
  workflow,
  temps = NULL,
  durations = NULL,
  ndraws = NULL,
  probs = c(0.025, 0.5, 0.975),
  by = NULL
)

Arguments

Details

Use this for curves (a handful of temperatures, dense in duration). For a 2-D survival heatmap call derive_tdt_landscape() instead.

Value

A list with elements summary (tibble of temp, duration, survival_median, survival_lower, survival_upper; plus the moderator column(s) for a grouped fit), draws_matrix (the raw posterior matrix ⁠[ndraws x grid points]⁠ of survival probabilities), and grid (the prediction grid tibble from new_tdt_grid()). Use get_surv_draws() for the long-tibble form of the draws.

Examples

## Not run: 
wf <- fit_4pl(d, ...)
pred <- predict_survival_curves(wf, temps = c(30, 32, 34))

## End(Not run)

Description

One-screen header: data shape, asymptote/response bounds, parameterisation (midpoint, or direct CTmax/z with its threshold and t_ref; for grouped direct fits, the moderator(s)), random-effect grouping (if any), and fit status (draws if fitted, "spec only" otherwise). Call summary() for parameter posteriors and HMC diagnostics, plot() for MCMC trace plots.

Usage

## S3 method for class 'bayes_tls'
print(x, ...)

Arguments

Value

The object, invisibly.

Examples

## Not run: 
wf <- fit_4pl(std)
print(wf)

## End(Not run)

Description

Computes a temperature-dependent repair rate using the Sharpe-Schoolfield formulation (Sharpe & Schoolfield, 1981). All Arrhenius temperatures inside this function are in Kelvin; the user-facing temp_celsius argument is in degrees Celsius and converted internally.

Usage

repair_rate_schoolfield(temp_celsius, TA, TAL, TAH, TL, TH, TREF, r_ref)

Arguments

Details

The Sharpe-Schoolfield rate at temperature T_K (in Kelvin) is

$r(T_K) = \frac{r_{ref} \cdot \exp\bigl(T_A (T_{ref}^{-1} - T_K^{-1})\bigr)} {1 + \exp\bigl(T_{AL} (T_K^{-1} - T_L^{-1})\bigr) + \exp\bigl(T_{AH} (T_H^{-1} - T_K^{-1})\bigr)}$

where the numerator is the Arrhenius enzymatic rate at the reference temperature, and the two terms in the denominator suppress the rate at temperatures below TL and above TH respectively (the low- and high-temperature inactivation arms).

Value

Numeric vector of repair rates at the supplied temperatures, in the same units as r_ref. Negative or non-finite values are coerced to zero.

Examples

# Default shrimp-style parameters from the prototype: optimum ~17 °C
repair_rate_schoolfield(
  temp_celsius = seq(5, 30, by = 5),
  TA  = 14065, TAL = 50000, TAH = 120000,
  TL  = 10.5 + 273.15, TH = 22.5 + 273.15, TREF = 17 + 273.15,
  r_ref = 0.01
)

Description

Chlorophyll-fluorescence ($F_v/F_m$) measurements on excised snow gum (Eucalyptus pauciflora) leaf sections before and 16–24 h after heat exposure, from Experiment 1 (post-heat light vs dark recovery) of Arnold et al. (2026). Short branches were cut from six mature trees grown outdoors in Canberra, ACT; $\sim$1 cm$^2$ leaf sections were dark-adapted, given an initial $F_v/F_m$, then submerged in a temperature-controlled water bath under sub-saturating light across a grid of assay temperatures (30–56 degrees C) and exposure durations (5–120 min). After heat, paired arrays were held for 90 min in moderate light (recovery = "Light") or in darkness (recovery = "Dark"); a final $F_v/F_m$ was taken 16–24 h later. The response is the continuous proportion fvfm_prop (post/pre ratio), modelled with a Beta likelihood. The model-ready frame for the leaf PSII case study (sublethal, continuous-proportion endpoint). The light/dark contrast is a two-group categorical moderator; in the source experiment post-heat light lowered apparent heat tolerance.

Usage

snowgum_psii

Format

A data frame with 394 rows and 8 variables:

Temp

Assay temperature (degrees C); 30, 35, 40, 44, 48, 52, 56.

Time

Exposure duration (minutes); 5, 15, 30, 60, 120.

recovery

Post-heat recovery light condition: "Dark" (darkness immediately after heat) or "Light" (90 min moderate light post-heat). A two-level moderator.

plant

Replicate mature tree (factor, 6 levels); the natural random-effect grouping.

meas_day

Assay day (factor, 2 levels). Two levels only, so a poor random-effect grouping; better treated as fixed or omitted.

initial_fvfm

$F_v/F_m$ measured before heat exposure.

final_fvfm

$F_v/F_m$ measured 16–24 h after heat exposure.

fvfm_prop

Retained PSII function, final_fvfm / initial_fvfm (a proportion in the unit interval; 0 indicates complete loss of measurable PSII function).

Details

Two of the 396 raw rows have post/pre $F_v/F_m$ marginally above 1 (both Dark, low dose) where a leaf measured slightly higher after heat than before; retained function cannot exceed 1, so these are treated as measurement noise and excluded, leaving 394 rows.

Source

Arnold PA, Harris RJ, Aitken SM, Hoek MM, Cook AM, Leigh A, Nicotra AB (2026) Towards a standard approach to investigating the thermal load sensitivity of photosystem II via chlorophyll fluorescence. bioRxiv doi:10.64898/2026.04.09.717599 (CC BY-NC 4.0), Experiment 1, snow gum slice. Raw file: system.file("extdata", "data_function_PSII_TDT_snowgum.csv", package = "bayesTLS").

Examples

std <- standardize_data(snowgum_psii, temp = "Temp", duration = "Time",
                        proportion = "fvfm_prop",
                        random_effects = "plant",
                        duration_unit = "minutes")

Description

Rewrites user column names into a single project-standard schema and attaches metadata used by every downstream fitting and prediction helper. This is the single entry point for raw data — everything else in the library assumes the output of this function.

Usage

standardize_data(
  data,
  temp,
  duration,
  n_total = NULL,
  n_surv = NULL,
  n_dead = NULL,
  survival = NULL,
  mortality = NULL,
  proportion = NULL,
  proportion_eps = 0.001,
  random_effects = NULL,
  duration_unit = "hours",
  temp_mean = NULL
)

Arguments

Details

Two response types are supported:

• Count data (binomial / beta-binomial): supply n_total plus exactly one of n_surv, n_dead, survival, or mortality. The other counts are derived and the standardised columns include n_total, n_surv, n_dead, survival. response_type is recorded as "count".

• Continuous proportion (Beta), e.g. a chlorophyll-fluorescence $F_v/F_m$ ratio with no denominator: supply proportion and omit the count arguments. The value is stored in survival (clamped into the open interval ⁠(proportion_eps, 1 - proportion_eps)⁠ so the Beta likelihood is finite); no n_total/n_surv columns are created. response_type is recorded as "proportion".

If the dataset spans multiple categories (life stages, species, populations, etc.), filter to one category before calling this function and fit a separate model per subset — the fitter does not estimate category-level effects.

Value

A tibble with the standardised columns plus a "tdt_meta" attribute storing temp_mean, duration_unit, random_effects, response_type ("count" or "proportion"), response_var (the response column name for a proportion fit, else NULL), and proportion_eps (the clamp half-width used for a proportion fit, else NULL).

Examples

# Count data
raw <- data.frame(
  temperature_C = rep(c(30, 32, 34), each = 4),
  exposure_h    = rep(c(1, 2, 4, 8), times = 3),
  n             = 30L,
  alive         = c(29, 28, 25, 5, 30, 27, 18, 2, 28, 22, 10, 1)
)
standardize_data(raw,
                 temp     = "temperature_C",
                 duration = "exposure_h",
                 n_total  = "n",
                 n_surv   = "alive")

# Continuous proportion (Beta) data
raw_p <- data.frame(
  temperature_C = rep(c(30, 32, 34), each = 4),
  exposure_h    = rep(c(1, 2, 4, 8), times = 3),
  fvfm_ratio    = c(0.95, 0.9, 0.7, 0.2, 0.92, 0.6, 0.3, 0, 0.8, 0.4, 0.1, 0)
)
standardize_data(raw_p,
                 temp       = "temperature_C",
                 duration   = "exposure_h",
                 proportion = "fvfm_ratio")

Description

Useful as an overlay on plotted posterior curves. Returns one row per temperature × duration combination, with the mean observed survival proportion across replicates and a standard error.

Usage

summarise_observed_survival(observed)

Arguments

Value

A tibble with columns temp, duration, survival_mean, survival_se, survival_lower, survival_upper, n_units, n_total_sum.

Examples

raw <- data.frame(T = rep(30, 6), hrs = rep(c(1, 5), each = 3),
                  n = 30, alive = c(28, 27, 29, 10, 12, 8))
d <- standardize_data(raw, "T", "hrs", n_total = "n", n_surv = "alive")
summarise_observed_survival(d)

Description

Delegates to summary() on the underlying brmsfit (via get_brmsfit()), so the returned object is a brms brmssummary with the population-level coefficient table, group-level standard deviations, family parameters, and HMC diagnostics already laid out by brms.

Usage

## S3 method for class 'bayes_tls'
summary(object, ...)

Arguments

Details

The high-level workflow context (data shape, centred-temperature mean, asymptote bounds, random-effect grouping, draw count) is available via print.bayes_tls(). For natural-scale 4PL parameters (low, up, k, z), use tdt_parameter_table(). For the TDT quantities (z, CTmax_1hr, optionally T_crit), use tls().

Value

A brms brmssummary object (brms handles printing).

Examples

## Not run: 
summary(wf)
summary(wf, prob = 0.9, robust = TRUE)

## End(Not run)

Description

Pulls the population-level posterior of the reparameterised 4PL parameters, transformed back to the natural scale via the model's bounds, as a one-row-per-parameter tibble with median + 95% CrI.

Usage

tdt_parameter_table(workflow, by = NULL)

Arguments

Details

For the midpoint parameterisation the rows are low, up, k, the mid intercept and temperature slope, and z = -1 / mid_temp. For the direct CTmax/z parameterisation the rows are low, up, k, CTmax (at the model's reference dose) and z. All quantities are read by evaluating brms::posterior_linpred(nlpar=) at temp_c = 0 and temp_c = 1 (so the result is parameterisation- and coding-agnostic – no coefficient-name parsing). For a fit whose CTmax/z vary by a moderator the table is returned per group, with the moderator column(s) prepended.

Value

A tibble with columns parameter, median, lower, upper (plus the moderator column(s) for a grouped fit).

Examples

## Not run: 
tdt_parameter_table(wf)

## End(Not run)

Description

Quantile wrapper with TDT-friendly defaults

Usage

tdt_quantile(x, probs = c(0.025, 0.5, 0.975))

Arguments

Value

Numeric vector of length length(probs).

Examples

tdt_quantile(rnorm(100))

Description

A minimal classic theme with bold axis titles, used by all plotting helpers in this library. Caller can override with + theme(...).

Usage

theme_tdt(base_size = 13)

Arguments

Value

A ggplot2 theme object.

Examples

library(ggplot2)
ggplot(mtcars, aes(mpg, hp)) + geom_point() + theme_tdt()

Description

One call that derives the classical TDT quantities — thermal sensitivity z, CTmax, and (for lethal endpoints) T_crit — per moderator group from a fitted joint 4PL, including hand-written brms models with moderators (sex, life stage, clone, ...) on any sub-parameter. The four sub-parameters are evaluated at a moderator x temperature grid with brms::posterior_linpred(), then z, CTmax and T_crit are derived from the same posterior draws (so they are mutually consistent and T_crit pairs CTmax and z per draw).

Usage

tls(
  object,
  by = NULL,
  params = "all",
  target_surv = "relative",
  t_ref = NULL,
  time_multiplier = NULL,
  lethal = FALSE,
  TC_rate_range = c(0.1, 1),
  temp = "temp_c",
  temp_mean = NULL,
  temp_grid = NULL,
  re_formula = NA,
  lower = 0,
  upper = 1,
  nlpars = c("lowraw", "upraw", "logk", "mid"),
  ndraws = NULL,
  newdata = NULL,
  probs = c(0.025, 0.5, 0.975),
  seed = NULL,
  mode = NULL,
  p = NULL
)

tls_z(object, ...)

tls_ctmax(object, ...)

tls_tcrit(object, ...)

Arguments

Details

tls() derives z, CTmax and T_crit from the same posterior draws, so on a common temperature grid they are mutually consistent for every threshold. z is the central-difference local z(T) (negative reciprocal of the local slope of log10(LT) on temperature) pooled over temp_grid; CTmax is the per-draw temperature at which the fitted LT curve crosses t_ref, from the exact closed form when the relative midpoint is linear and from a true per-draw inversion of the (possibly bent) curve otherwise. Under target_surv = "relative" (the default) log10(LT) = mid(T) is linear, so z is constant in temperature and equals -1 / b_mid_temp_c; under an absolute threshold with temperature effects on the asymptotes or steepness the curve bends and both the local z(T) and the inversion account for it, with no linear approximation.

Value

A tls object: ⁠$summary⁠ (per-group, per-quantity median + interval), ⁠$draws⁠ (per-group, per-quantity posterior draws), and ⁠$meta⁠.

Examples

## Not run: 
tls(joint_sex_fit, by = "sex", lethal = TRUE, temp_mean = 36.1)  # z, CTmax, T_crit per sex
tls(wf_leaf, params = "z")                                       # z only, workflow

## End(Not run)

Description

Two propagation methods on the Stage-2 fit:

• "delta" — delta-method standard errors for z and CTmax, with both Normal and t quantiles (the t-quantile is the small-sample correction for the few Stage-2 residual degrees of freedom). This is the method the bias simulation reports.

• "mvn" — slope-CI inversion for z (defined only when the slope CI is wholly negative) plus MVN simulation of the Stage-2 coefficients for CTmax and T_crit, and a predict.lm confidence band for the LT-vs-T line over temp_grid. This is the method the case studies report.

Usage

ts_ci(
  stage2,
  method = c("delta", "mvn"),
  level = 0.95,
  t_ref = 60,
  time_multiplier = 1,
  TC_rate_range = c(0.1, 1),
  temp_grid = NULL,
  n_sim = 1000,
  seed = 123
)

Arguments

Value

For "delta", a list with z and CTmax_1hr, each list(point, lower, upper, lower_t, upper_t, se), plus df_resid. As in ts_stage2(), the CTmax_1hr block reports CTmax at t_ref despite its name (the two coincide at the default t_ref = 60). For "mvn", a list with summary_ci (z/CTmax/T_crit bounds) and curve_ci (per-temp_grid line band).

Examples

d <- data.frame(
  temp = rep(c(30, 32, 34, 36, 38), each = 12),
  dur  = rep(c(1, 5, 15, 45, 135, 405), times = 10),
  surv = rbinom(60, 20, 0.4), tot = 20)
s2 <- ts_stage2(ts_stage1(d, "temp", "dur", "surv", "tot"))
ts_ci(s2, method = "delta")$z

Description

Median LT-vs-temperature line from a two-stage fit

Usage

ts_curve(stage2, temp_grid, time_multiplier = 1)

Arguments

Value

A tibble with temp and duration_median (minutes).

Examples

d <- data.frame(
  temp = rep(c(30, 34, 38), each = 12),
  dur  = rep(c(1, 5, 15, 45), times = 9),
  surv = rbinom(36, 20, 0.4), tot = 20)
ts_curve(ts_stage2(ts_stage1(d, "temp", "dur", "surv", "tot")),
         temp_grid = seq(30, 38, 1))

Description

Fits a separate logistic dose-response curve at each assay temperature and reads off log10(LT50) (the duration at 50% survival). The binomial family uses stats::glm; the beta-binomial family uses glmmTMB::glmmTMB with a betabinomial family (overdispersion at Stage 1).

Usage

ts_stage1(
  data,
  temp = "temp",
  duration = "duration",
  n_surv = "n_surv",
  n_total = "n_total",
  family = c("binomial", "betabinomial")
)

Arguments

Details

Two validity flags are returned so callers can choose their own success rule: finite_ok (finite coefficients, negative non-trivial slope) and bracket_ok (the fitted LT50 lies within the observed duration range, padded by 0.5 on the log10 scale). stage1_ok is their conjunction.

Value

A tibble with one row per temperature: temp, log10_lt50, se_log10_lt50, slope, phi (beta-binomial precision, else NA), finite_ok, bracket_ok, stage1_ok.

Examples

d <- data.frame(
  temp = rep(c(30, 34, 38), each = 12),
  dur  = rep(rep(c(1, 5, 15, 45), 3), times = 3),
  surv = rbinom(36, 20, 0.5), tot = 20)
ts_stage1(d, "temp", "dur", "surv", "tot", family = "binomial")

Description

Regresses Stage-1 log10(LT50) on assay temperature by ordinary least squares and derives the classical quantities. z = -1/slope; CTmax(t_ref) = (log10(t_ref) - intercept) / slope; T_crit follows the rate-multiplier definition, ⁠CTmax(1h) + z * mean(log10(TC_rate_range/100))⁠.

Usage

ts_stage2(
  stage1,
  t_ref = 60,
  time_multiplier = 1,
  TC_rate_range = c(0.1, 1),
  rows = c("stage1_ok", "finite_ok")
)

Arguments

Details

The Stage-2 regression is unweighted ordinary least squares: it ignores the Stage-1 standard errors (se_log10_lt50) by design. This is the generated-regressor bias the manuscript demonstrates, so se_log10_lt50 is exposed for diagnostics only and is not propagated into the Stage-2 quantities.

Value

list(fit, summary); fit is NULL if fewer than 3 valid Stage-1 estimates remain. summary has intercept, slope_T, z, CTmax_1hr, T_crit, r_squared, n_stage1, n_excluded. Note that despite its name CTmax_1hr reports CTmax at t_ref (the two coincide at the default t_ref = 60); the 1 h anchor used for T_crit is computed separately and is invariant to t_ref.

Examples

d <- data.frame(
  temp = rep(c(30, 32, 34, 36, 38), each = 12),
  dur  = rep(c(1, 5, 15, 45, 135, 405), times = 10),
  surv = rbinom(60, 20, 0.4), tot = 20)
s1 <- ts_stage1(d, "temp", "dur", "surv", "tot")
ts_stage2(s1)$summary

Description

Survival of zebrafish (Danio rerio) larvae assayed for upper thermal tolerance under three oxygen treatments, the model-ready frame for the oxygen-gradient case study. Diploid and triploid larvae were held at assay temperatures of 26 (control), 38, 39 and 40 degrees C for 3.8–240 minutes under hypoxia, normoxia or hyperoxia, and scored alive/dead. One row per assay group; oxygen is the categorical moderator. Fit CTmax and z as functions of oxygen (optionally ploidy) in one joint 4PL to compare thermal tolerance across the gradient with full posterior uncertainty.

Usage

zebrafish_o2

Format

A data frame with 905 rows and 10 variables:

cohort

Larval cohort identifier.

ploidy

Ploidy, a factor with levels diploid, triploid.

oxygen

Oxygen treatment, a factor with levels hypoxia, normoxia, hyperoxia (the modelling moderator).

o2_nominal

Nominal oxygen target as percent air saturation (25 / 100 / 225).

o2_measured

Measured oxygen level (percent air saturation).

temp

Target assay temperature (degrees C).

temp_measured

Measured assay temperature (degrees C).

duration_min

Exposure duration (minutes).

n_total

Number of larvae in the assay group.

n_surv

Number surviving after exposure and recovery.

Source

Saruhashi S, Boerrigter JGJ, Hooymans MHL, Sinclair BJ, Verberk WCEP (2026). Data and code for: Oxygen availability and oxygen delivery but not oxidative stress shape heat tolerance in diploid and triploid zebrafish larvae. Zenodo, doi:10.5281/zenodo.20075355 (distributed under CC BY 4.0). Associated article: Journal of Experimental Biology 229(10): jeb251548, doi:10.1242/jeb.251548. Raw file: system.file("extdata", "data_lethal_TDT_zebrafish_oxygen.csv", package = "bayesTLS").

Examples

## Not run: 
std <- standardize_data(zebrafish_o2, temp = "temp", duration = "duration_min",
                        n_total = "n_total", n_surv = "n_surv",
                        duration_unit = "minutes")
wf <- fit_4pl(std, ctmax = ~ 0 + oxygen, z = ~ 0 + oxygen, t_ref = 60)
tls(wf, by = "oxygen", lethal = TRUE)   # z, CTmax, T_crit per oxygen treatment

## End(Not run)