| Title: | Frequentist Inference for Thermal Load Sensitivity Models |
|---|---|
| Description: | A maximum-likelihood implementation of the thermal-load-sensitivity framework for thermal death-time modelling introduced by Noble, Arnold and Pottier in the 'bayesTLS' package, providing the frequentist counterpart to that Bayesian workflow. The modelling idea and the four-parameter logistic parameterisation are theirs; 'freqTLS' contributes a 'Template Model Builder' ('TMB') likelihood whose midpoint is written directly in terms of critical thermal maximum ('CTmax') and thermal sensitivity (z) so both headline quantities are estimable, and reports uncertainty through a unified trio of frequentist intervals -- Wald (delta), profile-likelihood, and bootstrap -- for binomial and beta-binomial survival counts and beta-distributed proportions. Column and formula interfaces support fixed and grouped designs plus limited independent random intercepts; prediction includes survival curves and deterministic heat-injury scenarios. Equivalence claims are restricted to the matched relative-threshold, constant-shape 'bayesTLS' configuration. |
| Authors: | Shinichi Nakagawa [aut, cre, cph] (ORCID: <https://orcid.org/0000-0002-7765-5182>), Pieter A. Arnold [aut] (ORCID: <https://orcid.org/0000-0002-6158-7752>, co-author of the bayesTLS framework), Patrice Pottier [aut] (ORCID: <https://orcid.org/0000-0003-2106-6597>, co-author of the bayesTLS framework), Daniel W. A. Noble [aut] (ORCID: <https://orcid.org/0000-0001-9460-8743>, senior author of the bayesTLS thermal-load-sensitivity framework) |
| Maintainer: | Shinichi Nakagawa <[email protected]> |
| License: | GPL (>= 3) |
| Version: | 0.1.0 |
| Built: | 2026-07-21 14:49:33 UTC |
| Source: | https://github.com/cran/freqTLS |
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 profile-likelihood
confidence intervals on those direct parameters.
aphid_tdtaphid_tdt
A data frame with 3041 rows and 7 variables:
Species, a factor with levels M_dirhodum
(Metopolophium dirhodum), S_avenae
(Sitobion avenae), R_padi
(Rhopalosiphum padi).
Age in days, a factor with levels 2, 6, 12.
Stress branch, a factor with levels heat
(34–40 degrees C), cold (-11 to -3 degrees C).
Assay temperature (degrees C).
Exposure duration (minutes).
Number of aphids treated.
Number surviving after treatment and recovery.
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 = "freqTLS").
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 speciesa <- 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
check_tls() re-runs the data-adequacy diagnostics
on a fitted profile_tls object, including the two post-fit checks that
fit_tls() cannot run before the model exists: whether any fitted CTmax is
extrapolated beyond the assayed temperatures (item 7), and whether phi has
reached the binomial limit (item 8). Each concern is emitted as a
cli::cli_warn(). Use it to audit a fit, or after suppressWarnings() around
fit_tls().
check_tls(fit)check_tls(fit)
fit |
A |
The profile-geometry diagnostics (items 9-12) are emitted by confint() and
profile() when those are called, not here, because they require the profile
likelihood.
Invisibly, a character vector of the diagnostic codes that fired.
The warning code returned by check_tls() identifies the next action:
temps: assay at least three distinct temperatures spanning the survival
transition.
durations / durations_per_temp: assay at least three distinct durations
per temperature, with times on both sides of the transition.
no_mortality: extend to hotter or longer exposures until mortality occurs.
all_mortality: add cooler or shorter exposures that retain survivors.
threshold: extend the design until observed survival straddles 0.5.
up_not_approached / low_not_approached: add milder / harsher conditions
approaching survival 1 / 0, or do not interpret that asymptote.
ctmax_extrapolated: expand the assayed temperature range to bracket CTmax;
otherwise report it explicitly as extrapolated.
phi_binomial_limit: consider the simpler binomial family.
After changing the design or family, refit and rerun check_tls(). For an
existing data set that cannot be augmented, use the warning to limit the
scientific claim; vignette("profile-likelihood") explains the strict
fallback = FALSE diagnostic and the default bootstrap recovery attempt.
d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) check_tls(fit)d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) check_tls(fit)
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.
clock_to_minutes(x)clock_to_minutes(x)
x |
Time value(s). |
Numeric vector of minutes.
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 minclock_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
confint() returns confidence intervals for the natural-scale parameters of a
fit_tls() model. Three methods are available:
## S3 method for class 'profile_tls' confint( object, parm = NULL, level = 0.95, method = c("profile", "wald", "bootstrap"), npoints = 30L, trace = FALSE, fallback = TRUE, nboot = 1000L, boot_seed = NULL, cores = 1L, ... )## S3 method for class 'profile_tls' confint( object, parm = NULL, level = 0.95, method = c("profile", "wald", "bootstrap"), npoints = 30L, trace = FALSE, fallback = TRUE, nboot = 1000L, boot_seed = NULL, cores = 1L, ... )
object |
A |
parm |
Character vector of target names (for example |
level |
Confidence level (default |
method |
One of |
npoints |
Number of grid points used per profile (default |
trace |
Logical; print inner-optimisation progress (profile/bootstrap). |
fallback |
Logical; when |
nboot |
Number of bootstrap replicates for |
boot_seed |
Optional integer seed making the bootstrap reproducible
without disturbing the caller's random stream (default |
cores |
Number of CPU cores for the bootstrap refits (default |
... |
Reserved; must be empty. |
method = "profile" (default) computes profile-likelihood confidence
intervals by inverting the likelihood-ratio test: the interval is
{psi : D(psi) <= qt(1 - alpha/2, df)^2}, found by stats::uniroot() on each
side of the MLE on the unconstrained internal coordinate, with the endpoints
transformed to the natural scale. The cutoff is the squared profile-t
quantile on df = n - p residual degrees of freedom (Bates-Watts profile-t),
not qchisq(level, 1); the two coincide as df -> Inf. These intervals are
prior-free and respect asymmetry. They are equivariant under monotone
reparameterisation, so the z interval equals exp() of the internal
log_z interval.
method = "wald" reuses the Phase-2 Wald path: estimate +/- t * se (with
t = qt(1 - alpha/2, df) on df = n - p) on the internal (link) scale,
back-transformed.
method = "bootstrap" returns prior-free parametric-bootstrap percentile
intervals: survival counts are regenerated at the observed design from the
fitted 4PL, the model is refitted nboot times, and the interval is the
percentile range of the replicate estimates. This is the likelihood-path
analogue of the bayesTLS posterior interval. It returns a finite interval
only when enough stable, non-degenerate refits remain.
When a profile does not close on one side, or the fitted Hessian is not
positive definite (pdHess = FALSE), confint() falls back to the parametric
bootstrap for the affected parameters (with a message). The fallback can still
return NA when too few valid refits remain. Set fallback = FALSE to keep
the strict profile behaviour, which returns NA on
the open side (never a fabricated bound) with a warning that the parameter is
weakly identified (see vignette("profile-likelihood")). The upper asymptote up has its own
coordinate beta_up under disjoint bounds but is not yet profiled, so it is
reported with the delta-method Wald interval under the profile/Wald methods,
with a message.
For a fit with a random intercept (CTmax ~ <fixed> + (1 | group)),
method = "profile" profiles the fixed-effect coordinates by re-running the
Laplace approximation at each grid point, which is slower than a fixed-effects
profile. Variance components keep their log-scale Wald intervals under the
profile method, and a non-closing random-effects profile falls back to Wald.
method = "bootstrap" instead redraws every active random-intercept block and
refits with the Laplace approximation, returning percentile intervals when
enough stable refits remain.
A tibble with one row per target and columns
parameter, conf.low, conf.high, estimate, level, method,
scale, and conf.status.
d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) confint(fit, "CTmax", method = "profile") confint(fit, "z", method = "profile")d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) confint(fit, "CTmax", method = "profile") confint(fit, "z", method = "profile")
derive_ctmax() inverts the fitted 4PL for temperature: it returns the
assay temperature at which survival equals a target surv after exposure
duration. By default surv is the relative midpoint threshold
(low + up) / 2 and duration is tref, so derive_ctmax(fit) reproduces
the fitted CTmax. Supplying an absolute surv gives the absolute-
threshold critical temperature (the analogue of the bayesTLS
extract_tdt() absolute mode), with the asymmetry correction
qlogis((surv - low) / (up - low)) / k built in.
derive_ctmax(object, surv = NULL, duration = NULL, group = NULL)derive_ctmax(object, surv = NULL, duration = NULL, group = NULL)
object |
A |
surv |
Target survival probability in |
duration |
Exposure duration(s) (native time unit; strictly positive).
Defaults to the fit's |
group |
Optional single group level (grouped fits only). |
Solving surv = low + (up - low) * plogis(-k (log10(duration) - mid)) with
mid = log10(tref) - (temp - CTmax) / z for the temperature gives
The target surv must lie strictly between low and up.
For a random-effects fit this is a population-level derived quantity; it does
not add a group BLUP.
A numeric vector of temperatures (degrees C), one per duration.
d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) derive_ctmax(fit) # ~ CTmax (relative, at tref) derive_ctmax(fit, surv = 0.5, duration = c(1, 4)) # absolute 50% survivald <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) derive_ctmax(fit) # ~ CTmax (relative, at tref) derive_ctmax(fit, surv = 0.5, duration = c(1, 4)) # absolute 50% survival
derive_lt() solves the fitted 4PL for the duration at which survival
crosses a target probability p at a given temperature (an "LT" / lethal-
time-style quantity, e.g. p = 0.5 gives the median survival time). Because
the threshold is interpreted relative to the asymptotes, the default
p = 0.5 returns the curve's midpoint duration, where
log10(duration) = mid exactly.
derive_lt(object, p = 0.5, temp, group = NULL)derive_lt(object, p = 0.5, temp, group = NULL)
object |
A |
p |
Target survival probability in |
temp |
Numeric temperature(s) at which to solve. |
group |
Optional single group level (grouped fits only). Required when the fit is grouped. |
Survival follows
p = low + (up - low) * plogis(-k (log10(duration) - mid)), so the duration
at which survival equals a target p solves
The target must lie strictly between low and up for a finite crossing;
otherwise the survival curve never reaches p and derive_lt() aborts with
an explanatory message (confidence-language, never silent).
For a random-effects fit this is a population-level derived quantity; it does
not add a group BLUP.
A numeric vector of durations (same length as temp) on the data's
native time unit.
d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) # Median survival duration at 36 C: derive_lt(fit, p = 0.5, temp = 36)d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) # Median survival duration at 36 C: derive_lt(fit, p = 0.5, temp = 36)
derive_tcrit() returns the rate-multiplier critical temperature T_crit: the
temperature at which the thermal-damage rate falls to a chosen low floor
rate. It is the maximum-likelihood analogue of the bayesTLS
extract_tdt() absolute-family T_crit, and follows directly from the fitted
CTmax and z:
Because rate < 100 makes log10(rate / 100) < 0 and z > 0, T_crit sits
below CTmax: it is the lower thermal threshold at which damage becomes
negligible (the temperature cutoff a heat-injury accumulation model treats as
"no damage").
derive_tcrit(object, rate = 1, group = NULL)derive_tcrit(object, rate = 1, group = NULL)
object |
A |
rate |
Damage-rate floor(s), a percentage of the lethal dose per hour
(strictly positive). A scalar or a vector; default |
group |
Optional single group level (grouped fits only; required when the fit is grouped). |
rate is a damage-rate floor expressed as a percentage of the lethal dose
per hour; bayesTLS brackets observed breakpoints with a default range of
0.1–1 %/hour. Unlike the Bayesian path, which samples rate to fold an
operational choice into the posterior, freqTLS treats rate as a fixed
input and returns the deterministic transform of the fitted CTmax and z
(combine their confidence intervals if you need to propagate uncertainty).
For a random-effects fit this is a population-level derived quantity; it does
not add a group BLUP.
T_crit assumes a lethal endpoint: it is a damage-accumulation concept, so
for sublethal endpoints (knockdown, photosynthetic failure) the steeper z
drives it implausibly low. derive_tcrit() says so, once per call.
A numeric vector of critical temperatures (degrees C), one per rate.
derive_ctmax() for the absolute-threshold critical temperature.
d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) derive_tcrit(fit, rate = c(0.1, 1)) # lower thermal thresholdsd <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) derive_tcrit(fit, rate = c(0.1, 1)) # lower thermal thresholds
diagnose_tdt_fit)The maximum-likelihood analogue of bayesTLS::diagnose_tdt_fit(): where the
Bayesian version reports Rhat / ESS / divergences, the freqTLS version reports
optimiser convergence, a positive-definite Hessian, and the gradient norm at
the optimum, with a single all_pass flag.
diagnose_tdt_fit(object)diagnose_tdt_fit(object)
object |
A |
A one-row tibble of convergence diagnostics.
raw <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) dat <- standardize_data( raw, temp = "temp", duration = "duration", n_total = "total", n_surv = "survived" ) fit <- fit_4pl(dat, family = "binomial", t_ref = 1, quiet = TRUE) diagnose_tdt_fit(fit)raw <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) dat <- standardize_data( raw, temp = "temp", duration = "duration", n_total = "total", n_surv = "survived" ) fit <- fit_4pl(dat, family = "binomial", t_ref = 1, quiet = TRUE) diagnose_tdt_fit(fit)
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). Only
dead is a valid freqTLS response: aggregate it to counts for the
beta-binomial lethal fit. The t_coma and prod columns are
retained to preserve the deposited record and provide study context; they
require time-to-event and reproductive-response models that freqTLS does not
fit. 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.
dsuzukiidsuzukii
A data frame with 1407 rows and 9 variables:
Unique individual identifier (temp-lvl-sex-rep).
Assay temperature (degrees C).
Exposure duration as a percentage of the estimated median time-to-coma from the authors' initial TDT curves.
Exposure duration (minutes).
Sex, a factor with levels F, M.
Replicate vial within a temperature x lvl x sex cell.
Reproductive productivity (offspring per female per day).
Mortality indicator: 1 = died, 0 = survived.
Time to heat coma (minutes); NA where no coma was
recorded for that individual.
Ø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 = "freqTLS").
# Lethal endpoint: aggregate per-individual deaths to cell counts, then # prepare the data for a beta-binomial 4PL. cells <- stats::aggregate( cbind(n_total = rep.int(1L, nrow(dsuzukii)), n_dead = dead) ~ temp + time + sex, data = dsuzukii, FUN = sum ) std <- standardize_data(cells, temp = "temp", duration = "time", n_total = "n_total", n_dead = "n_dead", duration_unit = "minutes")# Lethal endpoint: aggregate per-individual deaths to cell counts, then # prepare the data for a beta-binomial 4PL. cells <- stats::aggregate( cbind(n_total = rep.int(1L, nrow(dsuzukii)), n_dead = dead) ~ temp + time + sex, data = dsuzukii, FUN = sum ) std <- standardize_data(cells, temp = "temp", duration = "time", n_total = "n_total", n_dead = "n_dead", duration_unit = "minutes")
The frequentist twin of bayesTLS::extract_tdt(). Runs a parametric bootstrap
(via the freqTLS engine), derives the thermal-death-time quantities on each
replicate, and returns the same nested $z / $CTmax / $T_crit structure
(each a list of draws + summary). The per-replicate tables are the
frequentist analogue of posterior draws; *_median is the maximum-likelihood
point estimate and *_lower / *_upper are bootstrap percentiles.
extract_tdt( object, target_surv = "relative", lethal = FALSE, TC_rate_range = c(0.1, 1), nboot = 1000L, level = 0.95, seed = NULL, by = NULL )extract_tdt( object, target_surv = "relative", lethal = FALSE, TC_rate_range = c(0.1, 1), nboot = 1000L, level = 0.95, seed = NULL, by = NULL )
object |
A |
target_surv |
|
lethal |
If |
TC_rate_range |
Damage-rate floor range (percent of lethal dose per hour)
for |
nboot |
Number of bootstrap replicates (default 1000; smaller is faster). |
level |
Confidence level (default 0.95). |
seed |
Optional RNG seed for reproducible replicates / rate draws. |
by |
Optional name for the grouping column; defaults to the fit moderator. |
A list with $z, $CTmax, ($T_crit when lethal), and $meta.
Each quantity is list(draws = <tibble>, summary = <tibble>). Column names
follow bayesTLS: z_median/z_lower/z_upper for z; temp_median/temp_lower/ temp_upper for CTmax and T_crit; per-draw value columns are z / temp.
fit_4pl(), tls(), get_z_summary(), get_ctmax_summary()
raw <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) dat <- standardize_data( raw, temp = "temp", duration = "duration", n_total = "total", n_surv = "survived" ) fit <- fit_4pl(dat, family = "binomial", t_ref = 1, quiet = TRUE) tdt <- extract_tdt(fit, nboot = 10, seed = 1) tdt$CTmax$summaryraw <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) dat <- standardize_data( raw, temp = "temp", duration = "duration", n_total = "total", n_surv = "survived" ) fit <- fit_4pl(dat, family = "binomial", t_ref = 1, quiet = TRUE) tdt <- extract_tdt(fit, nboot = 10, seed = 1) tdt$CTmax$summary
The frequentist twin of bayesTLS::fit_4pl(). Consumes standardize_data()
output and fits the single-stage 4PL thermal death-time model, parameterised
directly in CTmax and thermal sensitivity (z), via the freqTLS TMB engine.
Returns a freq_tls workflow object; uncertainty (Wald / profile / bootstrap)
is computed on demand by the quantity twins (tls(), confint()).
fit_4pl( data, ctmax = NULL, z = NULL, up = NULL, low = NULL, k = NULL, by = NULL, threshold = c("relative", "absolute"), p = 0.5, t_ref = 60, bounds = c(0, 1), family = NULL, method = c("profile", "wald", "bootstrap"), start = NULL, control = list(), trace = FALSE, quiet = FALSE )fit_4pl( data, ctmax = NULL, z = NULL, up = NULL, low = NULL, k = NULL, by = NULL, threshold = c("relative", "absolute"), p = 0.5, t_ref = 60, bounds = c(0, 1), family = NULL, method = c("profile", "wald", "bootstrap"), start = NULL, control = list(), trace = FALSE, quiet = FALSE )
data |
Output of |
ctmax, z, up, low, k, by
|
Direct-mode formula interface; see
|
threshold |
|
p |
Survival level for the absolute threshold (default 0.5). |
t_ref |
Reference exposure time (in the data's |
bounds |
Asymptote range. Only |
family |
|
method |
Default interval method for downstream extraction
( |
start, control, trace, quiet
|
Passed to the engine |
A freq_tls object: a list with $fit (the engine fit), $data,
$formula, and $meta (threshold, t_ref, bounds, temp_mean, response_type,
family, grouped, moderators, method).
Before interpreting the fit, run check_tls(). Its help page gives a recovery
action for each data-adequacy warning; vignette("profile-likelihood") explains
strict open profiles and the default bootstrap fallback.
standardize_data(), make_4pl_formula(), fit_tls(), check_tls()
raw <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) dat <- standardize_data( raw, temp = "temp", duration = "duration", n_total = "total", n_surv = "survived" ) fit <- fit_4pl( dat, family = "binomial", t_ref = 1, method = "wald", quiet = TRUE ) coef(fit)raw <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) dat <- standardize_data( raw, temp = "temp", duration = "duration", n_total = "total", n_surv = "survived" ) fit <- fit_4pl( dat, family = "binomial", t_ref = 1, method = "wald", quiet = TRUE ) coef(fit)
fit_tls() fits the descending four-parameter logistic (4PL) thermal
death-time model to survival-count data, parameterised directly in CTmax
and z (thermal sensitivity) so that both headline quantities can be
profiled. Survival is modelled as a function of log10(duration); the
midpoint moves with temperature through CTmax and z (see
vignette("model-math")).
fit_tls( x, y, n, time, temp, group = NULL, family = c("beta_binomial", "binomial", "beta"), tref = 1, start = NULL, control = list(), trace = FALSE, quiet = FALSE, data = NULL )fit_tls( x, y, n, time, temp, group = NULL, family = c("beta_binomial", "binomial", "beta"), tref = 1, start = NULL, control = list(), trace = FALSE, quiet = FALSE, data = NULL )
x |
Either a data frame (column interface) or a |
y |
< |
n |
< |
time |
< |
temp |
< |
group |
< |
family |
One of |
tref |
Reference time at which |
start |
Optional named list of starting values on the internal
(unconstrained) scale, overriding the defaults. Names must match the
parameters in |
control |
List of optimiser controls; |
trace |
Logical; print optimiser progress. A shortcut for
|
quiet |
Logical; if |
data |
Used only in the formula interface: the data frame the
|
There are two equivalent interfaces. In the column interface, columns are
referenced with tidy evaluation: pass the bare column names of data (as in
dplyr), not strings. In the formula interface, pass a tls_bf() object
as x and the data frame as data; the brms/drmTMB-style grammar names the
response, the two axes, and the CTmax / log_z predictors. Both interfaces
feed the same likelihood engine, so a grouped formula fit and the matching
group = column fit are numerically identical.
An object of class c("profile_tls", "tls_fit"): a list with the
call, the resolved family, tref, group_levels, a data_summary, the
internal-scale MLE par, an estimates data frame of natural-scale
parameters with standard errors, the vcov of the internal coordinates,
the logLik, residual df, AIC, a convergence list
(code/pdHess/message), the name_map, and the underlying TMB obj,
optimiser opt, and sdreport.
Run check_tls() on the fitted object. Its help page maps every
data-adequacy warning to a concrete design or analysis response.
vignette("profile-likelihood") explains strict open profiles and the
default bootstrap recovery attempt.
check_tls(), confint.profile_tls()
d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) fit$estimates # The same fit through the formula interface: fit2 <- fit_tls( tls_bf(survived | trials(total) ~ time(duration) + temp(temp)), data = d, family = "binomial", tref = 1 )d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) fit$estimates # The same fit through the formula interface: fit2 <- fit_tls( tls_bf(survived | trials(total) ~ time(duration) + temp(temp)), data = d, family = "binomial", tref = 1 )
Format a point estimate plus confidence interval as a single string
format_interval(median, lower, upper, digits = 2)format_interval(median, lower, upper, digits = 2)
median, lower, upper
|
Numeric (scalar or vector). |
digits |
Integer rounding precision. |
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.
format_interval(5.123, 4.872, 5.401) format_interval(NA, 1, 2) # -> NAformat_interval(5.123, 4.872, 5.401) format_interval(NA, 1, 2) # -> NA
Extract the CTmax estimate(s)
get_ctmax(fit, conf.int = TRUE, conf.level = 0.95)get_ctmax(fit, conf.int = TRUE, conf.level = 0.95)
fit |
A |
conf.int |
Logical; include Wald |
conf.level |
Confidence level for the Wald interval (default |
A tibble of the CTmax row(s) from
tidy_parameters().
d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) get_ctmax(fit)d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) get_ctmax(fit)
Extract the shape parameters (low, up, k, and phi)
get_shape(fit, conf.int = TRUE, conf.level = 0.95)get_shape(fit, conf.int = TRUE, conf.level = 0.95)
fit |
A |
conf.int |
Logical; include Wald |
conf.level |
Confidence level for the Wald interval (default |
A tibble of the shape rows (low, up, k, and
phi for the beta-binomial family) from tidy_parameters().
d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) get_shape(fit)d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) get_shape(fit)
Extract the thermal-sensitivity (z) estimate(s)
get_z(fit, conf.int = TRUE, conf.level = 0.95)get_z(fit, conf.int = TRUE, conf.level = 0.95)
fit |
A |
conf.int |
Logical; include Wald |
conf.level |
Confidence level for the Wald interval (default |
A tibble of the z row(s) from
tidy_parameters().
d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) get_z(fit)d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) get_z(fit)
heat_injury_envelope() is the uncertainty counterpart of
predict_heat_injury(): it redraws the fitted curve parameters by parametric
bootstrap (the same machinery confint() uses), re-integrates the survival
trajectory under the temperature trace for each draw with the documented
dose-accumulation map, and returns a pointwise confidence band around the
point-estimate survival curve. The band is prior-free – it carries no
prior and makes no probability statement about the parameters; it is the
likelihood-path analogue of the bayesTLS posterior survival band, not a
credible band.
heat_injury_envelope( object, trace, group = NULL, target_surv = NULL, t_c = NULL, repair = NULL, irreversible = TRUE, nboot = 1000L, conf.level = 0.95, seed = NULL )heat_injury_envelope( object, trace, group = NULL, target_surv = NULL, t_c = NULL, repair = NULL, irreversible = TRUE, nboot = 1000L, conf.level = 0.95, seed = NULL )
object |
A |
trace |
A data frame with numeric columns |
group |
Optional single group level (grouped fits only; required when the fit is grouped). |
target_surv |
Optional absolute survival threshold defining one lethal
dose: a single probability strictly between the fitted lower and upper
asymptotes. |
t_c |
Optional damage-cutoff temperature (degrees C): at or below it the
damage rate is zero. |
repair |
Optional named list of Sharpe-Schoolfield repair parameters (see
Details); |
irreversible |
Logical; if |
nboot |
Number of bootstrap replicates (default |
conf.level |
Width of the pointwise confidence band (default |
seed |
Optional integer seed; when supplied the bootstrap is reproducible without disturbing the caller's random stream. |
A tibble with time, temp, survival (the
point-estimate trajectory from predict_heat_injury()), and conf.low /
conf.high (the pointwise parametric-bootstrap confidence band).
predict_heat_injury() for the point trajectory, plot_heat_injury()
to draw the band.
d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) trace <- data.frame(time = seq(0, 2, by = 0.1), temp = 34 + 6 * sin(seq(0, 2, by = 0.1))) heat_injury_envelope(fit, trace, nboot = 50, seed = 1)d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) trace <- data.frame(time = seq(0, 2, by = 0.1), temp = 34 + 6 * sin(seq(0, 2, by = 0.1))) heat_injury_envelope(fit, trace, nboot = 50, seed = 1)
Translates the bayesTLS-style direct-mode arguments (ctmax, z, up,
low, k, by) into the engine's tls_bf() tls_formula object. Supplying
a ctmax and/or z formula is the direct parameterisation; by is shorthand
for grouping CTmax and z by a single moderator (~ 0 + by).
make_4pl_formula( ctmax = NULL, z = NULL, up = NULL, low = NULL, k = NULL, by = NULL, family = "beta_binomial" )make_4pl_formula( ctmax = NULL, z = NULL, up = NULL, low = NULL, k = NULL, by = NULL, family = "beta_binomial" )
ctmax, z, up, low, k
|
One-sided formulas (or |
by |
Optional single moderator column name; shorthand for
|
family |
|
Following the freqTLS constant-shape invariant, the asymptotes and steepness
(up, low, k) default to shared (~ 1) so the temperature effect runs
through the midpoint (CTmax / z) only; pass an explicit formula to let a shape
vary. ctmax and z must produce the same fixed-effect model-matrix columns.
Supported random intercepts go inside the ctmax/z/low/k formulas;
up random effects are not supported. For example,
ctmax = ~ 1 + (1 | batch) keeps the same intercept-only fixed design as the
default z = ~ 1 while adding a CTmax random intercept.
A tls_formula object (as built by tls_bf()).
fit_4pl(), tls_bf(), standardize_data()
make_4pl_formula() make_4pl_formula(by = "population", family = "binomial")make_4pl_formula() make_4pl_formula(by = "population", family = "binomial")
plot_confidence_eye() draws the freqTLS Confidence Eye for one or more
headline parameters (CTmax, z, or any other confint.profile_tls()
target, including grouped names). It is a HORIZONTAL forest display: each
parameter (and group level) is a row, the parameter value runs along the
x-axis, and the confidence interval is a short, wide pale lens with a
hollow point estimate. The shallow horizontal lens reads as a confidence
interval, never a posterior density – freqTLS intervals are likelihood
confidence intervals, so the wording is "confidence", never
"posterior". The layout follows the gllvmTMB / drmTMB
Confidence-Eye contract.
plot_confidence_eye( fit, parm = c("CTmax", "z"), method = c("profile", "wald", "bootstrap"), level = 0.95, style = c("eye", "line"), raw_data = TRUE, fallback = TRUE, nboot = 1000L, boot_seed = NULL, cores = 1L, ... )plot_confidence_eye( fit, parm = c("CTmax", "z"), method = c("profile", "wald", "bootstrap"), level = 0.95, style = c("eye", "line"), raw_data = TRUE, fallback = TRUE, nboot = 1000L, boot_seed = NULL, cores = 1L, ... )
fit |
A |
parm |
Character vector of target parameter names. Defaults to
|
method |
One of |
level |
Confidence level (default |
style |
One of |
raw_data |
Logical; overlay observed assay temperatures as a rug on
temperature-scale rows (default |
fallback, nboot, boot_seed, cores
|
Forwarded to |
... |
Reserved; must be empty. |
When a profile does not close (conf.status is "open_lower",
"open_upper", "open_both", or a bound is NA), no lens is drawn for that
row: a hollow point marks the estimate and the subtitle flags the open
interval. The eye is never fabricated from an open profile. A
"wald_fallback" interval (e.g. up) still gets a lens, with the source
noted in the caption.
With raw_data = TRUE (default), the observed assay temperatures are drawn as
a rug beneath any temperature-scale row (CTmax), showing the data support
and flagging extrapolation when CTmax sits outside the assayed range.
Parameters on different scales (temperature for CTmax, a positive
multiplier for z) are stacked in separate panels with a free x-axis.
A ggplot object.
d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) plot_confidence_eye(fit, parm = c("CTmax", "z"))d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) plot_confidence_eye(fit, parm = c("CTmax", "z"))
plot_heat_injury() draws the point-estimate survival trajectory from
predict_heat_injury() inside the pointwise parametric-bootstrap confidence
band from heat_injury_envelope(). The band is prior-free – a confidence
band, never a posterior / credible band (the project's honest-uncertainty
contract).
plot_heat_injury( object, trace, group = NULL, target_surv = NULL, t_c = NULL, repair = NULL, irreversible = TRUE, nboot = 1000L, conf.level = 0.95, seed = NULL, time_div = 1, xlab = "Time", ylab = "Survival" )plot_heat_injury( object, trace, group = NULL, target_surv = NULL, t_c = NULL, repair = NULL, irreversible = TRUE, nboot = 1000L, conf.level = 0.95, seed = NULL, time_div = 1, xlab = "Time", ylab = "Survival" )
object |
A |
trace |
A data frame with numeric columns |
group |
Optional single group level (grouped fits only; required when the fit is grouped). |
target_surv |
Optional absolute survival threshold defining one lethal
dose: a single probability strictly between the fitted lower and upper
asymptotes. |
t_c |
Optional damage-cutoff temperature (degrees C): at or below it the
damage rate is zero. |
repair |
Optional named list of Sharpe-Schoolfield repair parameters (see
Details); |
irreversible |
Logical; if |
nboot |
Number of bootstrap replicates (default |
conf.level |
Width of the pointwise confidence band (default |
seed |
Optional integer seed; when supplied the bootstrap is reproducible without disturbing the caller's random stream. |
time_div |
Optional positive divisor applied to |
xlab, ylab
|
Axis labels. |
A ggplot object.
heat_injury_envelope() for the band data, predict_heat_injury()
for the point trajectory.
d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) trace <- data.frame(time = seq(0, 2, by = 0.1), temp = 34 + 6 * sin(seq(0, 2, by = 0.1))) plot_heat_injury(fit, trace, nboot = 50, seed = 1)d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) trace <- data.frame(time = seq(0, 2, by = 0.1), temp = 34 + 6 * sin(seq(0, 2, by = 0.1))) plot_heat_injury(fit, trace, nboot = 50, seed = 1)
plot_survival_curves() draws the fitted survival probability as a function
of exposure duration (on a log10 x-axis), one curve per temperature, with the
observed survival proportions overlaid as points. For a grouped fit the
curves are faceted by group.
plot_survival_curves(fit, temps = NULL, times = NULL, ...)plot_survival_curves(fit, temps = NULL, times = NULL, ...)
fit |
A |
temps |
Numeric vector of temperatures to draw curves for. Defaults to the distinct observed temperatures (capped at a readable number). |
times |
Numeric vector of durations to evaluate the smooth curve over. Defaults to a log-spaced sequence over the observed duration range. |
... |
Reserved; must be empty. |
A ggplot object.
d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) plot_survival_curves(fit)d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) plot_survival_curves(fit)
plot_survival_surface() draws the fitted survival probability as a filled
heatmap over a temperature-by-duration grid, with contour lines, using
predict_survival_surface(). Duration is on a log10 axis. For a grouped fit
the surface is faceted by group.
plot_survival_surface(fit, temps = NULL, times = NULL, contour = TRUE, ...)plot_survival_surface(fit, temps = NULL, times = NULL, contour = TRUE, ...)
fit |
A |
temps, times
|
Numeric grids passed to |
contour |
Logical; overlay contour lines (default |
... |
Reserved; must be empty. |
A ggplot object.
d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) plot_survival_surface(fit)d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) plot_survival_surface(fit)
plot_tdt_curve() draws the duration at which survival crosses a target
probability p (default the relative midpoint, p = 0.5) against
temperature – the classic thermal-death-time line, here read directly off the
fitted 4PL via derive_lt(). Time is shown on a log10 axis. For a grouped fit
a line is drawn per group.
plot_tdt_curve(fit, p = 0.5, temps = NULL, ...)plot_tdt_curve(fit, p = 0.5, temps = NULL, ...)
fit |
A |
p |
Target survival probability for the threshold (default |
temps |
Numeric vector of temperatures. Defaults to a sequence over the observed temperature range. |
... |
Reserved; must be empty. |
A ggplot object.
d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) plot_tdt_curve(fit)d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) plot_tdt_curve(fit)
plot() for a "profile_tls_profile" object draws the likelihood-ratio
deviance curve against the natural-scale parameter. A dotted horizontal line
marks the profile-t cutoff qt(1 - alpha/2, df)^2; a solid vertical line marks the
point estimate; dashed vertical lines mark the interval endpoints when they
are finite. The wording is deliberately "confidence" – this
is a likelihood curve, never a posterior. A non-closing side is
annotated rather than drawn as a closed bound.
## S3 method for class 'profile_tls_profile' plot(x, ...)## S3 method for class 'profile_tls_profile' plot(x, ...)
x |
A |
... |
Reserved; must be empty. |
This is the per-parameter profile curve; the full Confidence-Eye interval displays are added in Phase 4.
A ggplot object (invisibly when printed for its side effect).
d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) plot(profile(fit, "CTmax"))d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) plot(profile(fit, "CTmax"))
predict_heat_injury() is the deterministic, maximum-likelihood prediction
analogue of bayesTLS::predict_heat_injury(): given a fitted thermal-load-
sensitivity curve and a temperature time-series (a "trace"), it accumulates
thermal damage as a fraction of the lethal dose and reads survival back off the
fitted 4PL. It does not fit an injury or repair model – fitting injury /
repair dynamics remains a bayesTLS concern (the complementary boundary);
predict_heat_injury() only predicts injury from the already-fitted survival
curve. For a random-effects fit it uses the population curve and does not add
a fitted group BLUP.
predict_heat_injury( object, trace, group = NULL, target_surv = NULL, t_c = NULL, repair = NULL, irreversible = TRUE )predict_heat_injury( object, trace, group = NULL, target_surv = NULL, t_c = NULL, repair = NULL, irreversible = TRUE )
object |
A |
trace |
A data frame with numeric columns |
group |
Optional single group level (grouped fits only; required when the fit is grouped). |
target_surv |
Optional absolute survival threshold defining one lethal
dose: a single probability strictly between the fitted lower and upper
asymptotes. |
t_c |
Optional damage-cutoff temperature (degrees C): at or below it the
damage rate is zero. |
repair |
Optional named list of Sharpe-Schoolfield repair parameters (see
Details); |
irreversible |
Logical; if |
One lethal dose is the thermal load that drives survival to a target set by
target_surv. With the default target_surv = NULL the target is the
project-default relative threshold – the curve midpoint (low + up) / 2.
At temperature T the lethal time to the target is
LT(T) = tref * 10^((CTmax - T) / z - q / k), with
q = qlogis((target_surv - low) / (up - low)) (so q = 0 at the midpoint, the
same quantity as derive_lt() at survival target_surv). The instantaneous
damage rate is 1 / LT(T) (lethal doses per time unit). Cumulative dose is
accumulated by forward Euler over the trace, using the actual per-step time
increments:
where . Survival is read back from the 4PL by
treating the accumulated dose as an equivalent log10-time:
survival(D) = low + (up - low) * plogis(-k * log10(D) + q), so D = 1 (one
lethal dose) reaches exactly target_surv – the relative midpoint
(low + up) / 2 by default, or an absolute survival threshold when
target_surv is supplied.
The trace time and the fit's duration / tref must share a time unit
(the damage rate is per that unit). With irreversible = TRUE (default)
survival is monotone non-increasing. A damage cutoff t_c (for example from
derive_tcrit()) sets the damage rate to zero at or below t_c.
This integrator is forward Euler (left-endpoint, per actual step), not the
single-dt scheme some implementations use; irregular traces are integrated
with their real increments.
If repair is supplied, a Sharpe-Schoolfield repair rate is subtracted each
step (scaled by the current survival fraction, so repair shrinks as the
population dies). The repair parameters are a user-supplied scenario layer:
they are not identified by the survival data the model was fitted to, so
predict_heat_injury() warns when they are used. repair is a named list with
r_ref, t_a, t_al, t_ah, t_l, t_h, t_ref, with the four reference
temperatures in Kelvin.
A data frame with columns time, temp, dose (cumulative, as a
fraction of the lethal dose), injury (dose * 100, percent), and
survival.
derive_lt() for the lethal time, derive_tcrit() for a damage
cutoff temperature.
d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) trace <- data.frame(time = seq(0, 2, by = 0.05), temp = 34 + 6 * sin(seq(0, 2, by = 0.05))) head(predict_heat_injury(fit, trace))d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) trace <- data.frame(time = seq(0, 2, by = 0.05), temp = 34 + 6 * sin(seq(0, 2, by = 0.05))) head(predict_heat_injury(fit, trace))
The frequentist twin of bayesTLS::predict_survival_curves(). Evaluates the
fitted 4PL survival probability over a temperature-by-duration grid and adds
parametric-bootstrap confidence bands. For random-effects fits the curves are
population-level: random intercepts are integrated during bootstrap refits,
but no fitted group BLUP is added to the reported curve.
predict_survival_curves( object, temps = NULL, durations = NULL, nboot = 500L, level = 0.95, seed = NULL, by = NULL )predict_survival_curves( object, temps = NULL, durations = NULL, nboot = 500L, level = 0.95, seed = NULL, by = NULL )
object |
A |
temps |
Temperatures to predict at (default: the observed assay temps). |
durations |
Exposure durations (default: 100 points log-spaced over the observed range, in the data's duration unit). |
nboot |
Number of bootstrap replicates for the bands (default 500). |
level |
Confidence level (default 0.95). |
seed |
Optional RNG seed. |
by |
Optional name for the grouping column. |
A freq_surv_curves object: $summary (a tibble of
[<group>,] temp, duration, survival_lower, survival_median, survival_upper)
and $meta.
fit_4pl(), predict_survival_surface(), tls()
raw <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) dat <- standardize_data( raw, temp = "temp", duration = "duration", n_total = "total", n_surv = "survived" ) fit <- fit_4pl(dat, family = "binomial", t_ref = 1, quiet = TRUE) curves <- predict_survival_curves( fit, temps = c(34, 36), durations = c(1, 4), nboot = 10, seed = 1 ) curves$summaryraw <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) dat <- standardize_data( raw, temp = "temp", duration = "duration", n_total = "total", n_surv = "survived" ) fit <- fit_4pl(dat, family = "binomial", t_ref = 1, quiet = TRUE) curves <- predict_survival_curves( fit, temps = c(34, 36), durations = c(1, 4), nboot = 10, seed = 1 ) curves$summary
predict_survival_surface() evaluates the fitted survival probability on a
factorial grid of temperatures by durations, returning a long data frame
suitable for a heatmap or contour plot (see plot_survival_surface()).
For random-effects fits this helper returns population-level predictions
(random intercepts set to zero); use predict(..., re.form = "conditional")
for known-group conditional predictions. General continuous fixed designs
require predict() with their covariate columns supplied in newdata.
predict_survival_surface(object, temps = NULL, times = NULL, group = NULL)predict_survival_surface(object, temps = NULL, times = NULL, group = NULL)
object |
A |
temps |
Numeric vector of temperatures. Defaults to a length-60 sequence spanning the fit's observed temperature range. |
times |
Numeric vector of durations (strictly positive). Defaults to a length-60 log-spaced sequence spanning the fit's observed duration range. |
group |
Optional single group level (grouped fits only). When |
A long data.frame with columns temp, duration, survival
(and group when the fit is grouped).
d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) head(predict_survival_surface(fit, temps = c(34, 36, 38), times = c(1, 2, 4)))d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) head(predict_survival_surface(fit, temps = c(34, 36, 38), times = c(1, 2, 4)))
predict() evaluates the fitted four-parameter logistic (4PL) thermal-load-
sensitivity model at new temperature-by-duration cells, using exactly the
same forward map as the TMB engine in src/profile_tls.cpp:
## S3 method for class 'profile_tls' predict( object, newdata, type = c("survival", "link", "midpoint"), re.form = c("population", "conditional"), ... )## S3 method for class 'profile_tls' predict( object, newdata, type = c("survival", "link", "midpoint"), re.form = c("population", "conditional"), ... )
object |
A |
newdata |
A data frame with numeric columns |
type |
One of |
re.form |
How to handle fitted random intercepts. |
... |
Reserved; must be empty. |
Three response types are available:
"survival" (default) returns the fitted survival probability in (0, 1).
"link" returns the logit of the survival probability,
qlogis(survival).
"midpoint" returns the temperature-dependent 4PL midpoint mid on the
log10(duration) axis (constant within a temperature, so the duration
column is ignored for this type but a temp column is still required).
newdata must also contain every predictor used by the fitted fixed-effect
designs for CTmax, log_z, low, up, or log_k. For a grouped column-
interface fit, supply group with values from the fitted group_levels.
For a formula fit, a literal tls_bf() call preserves the fixed-design
formulas needed to rebuild transformed or interacted terms. If the model was
instead passed through a formula-object variable, predict() can rebuild
direct numeric design columns but asks the user to refit with a literal
tls_bf() call when a transformed or interacted design cannot be recovered
safely.
A numeric vector with one element per row of newdata. Survival
values lie in (0, 1).
d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) nd <- expand.grid(temp = c(34, 36, 38), duration = c(1, 2, 4)) predict(fit, nd, type = "survival") # A continuous predictor used by CTmax and log_z must also be in newdata. d$x <- rep(c(-1, 1), length.out = nrow(d)) fit_x <- fit_tls( tls_bf(survived | trials(total) ~ time(duration) + temp(temp), CTmax ~ x, log_z ~ x), data = d, family = "binomial", tref = 1 ) predict(fit_x, data.frame(temp = 36, duration = 2, x = c(-1, 1))) # Choose population or fitted-group prediction explicitly for an RE fit. dre <- simulate_tls(family = "binomial", CTmax = 36, z = 4, re_sd = 1, n_re_groups = 8, seed = 2) fit_re <- fit_tls( tls_bf(survived | trials(total) ~ time(duration) + temp(temp), CTmax ~ 1 + (1 | colony)), data = dre, family = "binomial", tref = 1 ) colony <- as.character(ranef(fit_re)$group[1]) nd_re <- data.frame(temp = 36, duration = 2, colony = colony) predict(fit_re, nd_re, re.form = "population") predict(fit_re, nd_re, re.form = "conditional")d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) nd <- expand.grid(temp = c(34, 36, 38), duration = c(1, 2, 4)) predict(fit, nd, type = "survival") # A continuous predictor used by CTmax and log_z must also be in newdata. d$x <- rep(c(-1, 1), length.out = nrow(d)) fit_x <- fit_tls( tls_bf(survived | trials(total) ~ time(duration) + temp(temp), CTmax ~ x, log_z ~ x), data = d, family = "binomial", tref = 1 ) predict(fit_x, data.frame(temp = 36, duration = 2, x = c(-1, 1))) # Choose population or fitted-group prediction explicitly for an RE fit. dre <- simulate_tls(family = "binomial", CTmax = 36, z = 4, re_sd = 1, n_re_groups = 8, seed = 2) fit_re <- fit_tls( tls_bf(survived | trials(total) ~ time(duration) + temp(temp), CTmax ~ 1 + (1 | colony)), data = dre, family = "binomial", tref = 1 ) colony <- as.character(ranef(fit_re)$group[1]) nd_re <- data.frame(temp = 36, duration = 2, colony = colony) predict(fit_re, nd_re, re.form = "population") predict(fit_re, nd_re, re.form = "conditional")
profile() computes the profile-likelihood deviance curve for one scalar
target of a fit_tls() model. For the target it fixes the corresponding
internal (unconstrained) coordinate on a grid, re-optimises the remaining
coordinates at each grid point, and returns the deviance
D = 2 * (logLik_hat - logLik_profile) together with the profile-t cutoff and
the profile-likelihood confidence interval. Because the profile is taken on the
unconstrained coordinate and the endpoints are then transformed by a monotone
function, the interval is exactly equivariant: the z interval equals exp()
of the internal log_z interval (the headline equivariance check).
## S3 method for class 'profile_tls' profile(fitted, parm, level = 0.95, npoints = 30L, trace = FALSE, ...) ## S3 method for class 'profile_tls_profile' print(x, digits = 4, ...)## S3 method for class 'profile_tls' profile(fitted, parm, level = 0.95, npoints = 30L, trace = FALSE, ...) ## S3 method for class 'profile_tls_profile' print(x, digits = 4, ...)
fitted |
A |
parm |
A single target name (see Targets). |
level |
Confidence level for the interval and the cutoff line (default
|
npoints |
Number of grid points for the deviance curve (default |
trace |
Logical; print inner-optimisation progress. |
... |
Reserved; must be empty. |
x |
A |
digits |
Number of significant digits for the printed summary (default
|
The algorithm is a map-refit profile: the target coordinate is fixed with
TMB's map mechanism and the rest re-optimised, mirroring the bracket-then-
stats::uniroot() endpoint solver in drmTMB::R/profile.R:2314-2373. See
vignette("profile-likelihood").
An object of class "profile_tls_profile": a list with parm,
profile_value (grid on the natural scale), deviance, estimate,
conf.low, conf.high, conf.status, cutoff, level, scale, and
transformation.
print(profile_tls_profile): Print a compact summary of the profile.
| Target | Profiled coordinate | Endpoint transform |
CTmax, CTmax:<grp> |
beta_CT[g] |
identity |
z, z:<grp> |
beta_logz[g] |
exp |
log_z, log_z:<grp> |
beta_logz[g] |
identity |
low |
beta_low |
plogis |
k |
beta_logk |
exp |
phi |
log_phi |
exp |
up |
(Wald/delta fallback) | -- |
dCTmax:<a>-<b>, dlog_z:<a>-<b> |
contrast recoding | identity |
Under the disjoint-bounds parameterisation up = up_min + up_w * plogis(beta_up)
has its own coordinate beta_up, but freqTLS does not yet profile it (the profile
path is wired for low but not up — symmetric work, simply not implemented).
freqTLS falls back to the delta-method Wald interval for up
and says so. Group contrasts (dCTmax, dlog_z) are profiled
directly by recoding the design so the contrast is itself a coordinate.
d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) pc <- profile(fit, "CTmax") pc$conf.low pc$conf.highd <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) pc <- profile(fit, "CTmax") pc$conf.low pc$conf.high
ranef() returns the predicted random intercepts (the conditional modes /
BLUPs) with their conditional standard errors, for a fit with a random
intercept on any of CTmax, log_z, low, or log_k
(<param> ~ <fixed> + (1 | group)). It errors for a fixed-effects-only fit.
Each BLUP is a deviation on its coordinate's internal scale: CTmax in degrees
C, log_z on log(z), low on logit(low), log_k on log(k). When several
REs are present the rows are stacked in CTmax, log_z, low, log_k order.
ranef(object, ...) ## S3 method for class 'profile_tls' ranef(object, ...)ranef(object, ...) ## S3 method for class 'profile_tls' ranef(object, ...)
object |
A |
... |
Reserved; must be empty. |
A tibble with one row per group level (per RE term):
group, term ("CTmax", "log_z", "low", or "log_k"), estimate (the
BLUP), and std.error (the conditional SE).
d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, re_sd = 1.5, n_re_groups = 12, seed = 42) fit <- fit_tls( tls_bf(survived | trials(total) ~ time(duration) + temp(temp), CTmax ~ 1 + (1 | colony)), data = d, family = "binomial", tref = 1) ranef(fit)d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, re_sd = 1.5, n_re_groups = 12, seed = 42) fit <- fit_tls( tls_bf(survived | trials(total) ~ time(duration) + temp(temp), CTmax ~ 1 + (1 | colony)), data = d, family = "binomial", tref = 1) ranef(fit)
Replicate lethal-TDT trials for brown shrimp (Crangon crangon). Each row is one tank of individuals exposed to a fixed assay temperature for a fixed duration; the response is the proportion that died. The model-ready frame for Case Study 1 (lethal endpoint).
shrimp_lethalshrimp_lethal
A data frame with 148 rows and 6 variables:
Experiment date (use as a grouping factor).
Holding-tank identifier (use as a grouping factor).
Assay temperature (degrees C).
Exposure duration (hours).
Number of individuals in the trial.
Proportion that died during the trial
(deaths / N_individuals_after_trial), in the unit interval.
Consumed by
standardize_data(mortality = "Mortality_after_trial").
Brown shrimp lethal-TDT assay (Case Study 1), obtained from the
bayesTLS package distribution by Noble, Arnold, and Pottier (2026),
licensed CC BY 4.0. freqTLS retains the mortality proportion and documents
its count reconstruction above. Raw file:
system.file("extdata", "data_lethal_TDT_brown_shrimp.csv", package = "freqTLS").
std <- standardize_data(shrimp_lethal, temp = "Temperature_assay", duration = "Duration_exposure_hours", n_total = "N_individuals_after_trial", mortality = "Mortality_after_trial", random_effects = c("Date", "Tank"), duration_unit = "hours")std <- standardize_data(shrimp_lethal, temp = "Temperature_assay", duration = "Duration_exposure_hours", n_total = "N_individuals_after_trial", mortality = "Mortality_after_trial", random_effects = c("Date", "Tank"), duration_unit = "hours")
Sublethal TDT trials for brown shrimp (Crangon crangon): each cup of individuals contributes the elapsed time to loss of response to touch (knockdown) at a fixed assay temperature. Cleaned from the raw clock-time records (excluded rows dropped; start/stop times parsed to elapsed minutes).
shrimp_sublethalshrimp_sublethal
A data frame with 299 rows and 5 variables:
Assay temperature (degrees C).
Time to knockdown (minutes).
Experiment date (grouping factor).
Holding-tank identifier (grouping factor).
Cup identifier, Trial_ID_Sample (grouping factor).
Brown shrimp sublethal time-to-knockdown assay (Case Study 1,
sublethal endpoint), obtained from the bayesTLS package distribution
by Noble, Arnold, and Pottier (2026), licensed CC BY 4.0. freqTLS dropped
excluded rows and converted the clock times to elapsed minutes. Raw file:
system.file("extdata", "data_sublethal_TDT_brown_shrimp.csv", package = "freqTLS").
simulate_tls() draws survival counts from the locked data-generating
process used throughout the freqTLS test suite and benchmarks. It builds a
factorial grid of temperatures by durations by replicates, computes the true
survival probability under the direct-CTmax/z 4PL (the same forward map as
the TMB engine in src/profile_tls.cpp), and draws binomial or beta-binomial
counts. The simulating truth is attached as attr(, "truth").
simulate_tls( temps = seq(30, 42, by = 2), times = c(0.5, 1, 2, 4, 8), reps = 3, n = 20, low = 0.02, up = 0.98, k = 5, CTmax = 36, z = 3, phi = NULL, family = c("binomial", "beta_binomial", "beta"), group = NULL, re_sd = NULL, re_sd_z = NULL, re_sd_low = NULL, re_sd_logk = NULL, n_re_groups = NULL, re_group_name = "colony", tref = 1, seed = NULL )simulate_tls( temps = seq(30, 42, by = 2), times = c(0.5, 1, 2, 4, 8), reps = 3, n = 20, low = 0.02, up = 0.98, k = 5, CTmax = 36, z = 3, phi = NULL, family = c("binomial", "beta_binomial", "beta"), group = NULL, re_sd = NULL, re_sd_z = NULL, re_sd_low = NULL, re_sd_logk = NULL, n_re_groups = NULL, re_group_name = "colony", tref = 1, seed = NULL )
temps |
Numeric vector of assay temperatures (degrees C). |
times |
Numeric vector of exposure durations (native unit, e.g. hours). |
reps |
Number of replicate observations per temperature-by-duration cell (per group). |
n |
Number of individuals per observation (binomial size). |
low, up
|
Lower and upper survival asymptotes ( |
k |
Steepness of the logistic on the |
CTmax |
Critical thermal maximum at |
z |
Thermal sensitivity ( |
phi |
Dispersion (sum of Beta shapes); |
family |
One of |
group |
Optional atomic vector of group labels (character, factor,
or numeric). When supplied, |
re_sd |
Optional standard deviation of a random intercept on |
re_sd_z |
Optional standard deviation of a random intercept on
|
re_sd_low, re_sd_logk
|
Optional standard deviations of random
intercepts on the lower asymptote |
n_re_groups |
Number of random-effect groups (required with any |
re_group_name |
Name of the grouping column added to the output for the
random-effect mode (default |
tref |
Reference time at which |
seed |
Optional integer seed for reproducibility. |
phi conventionFor the beta-binomial family, phi is the sum of the Beta shape
parameters: counts are drawn as prob <- rbeta(a = p * phi, b = (1 - p) * phi) followed by rbinom(n, prob). The Beta mean is p and its variance is
p (1 - p) / (phi + 1), so larger phi means less overdispersion and the
binomial is recovered as phi -> Inf. This matches the engine's
parameterisation in beta_binomial_tls().
A base data.frame with columns temp, duration, the true
probability p, and (when grouped) group. The count families
(binomial, beta_binomial) add total and survived; the beta family
instead adds a single continuous proportion column prop. The
data-generating parameters are attached as attr(, "truth").
d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) head(d) attr(d, "truth")$CTmaxd <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) head(d) attr(d, "truth")$CTmax
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.
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 )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 )
data |
Raw data frame or tibble. |
temp |
Column name of the assay temperature (°C). |
duration |
Column name of the exposure duration. The unit is
whatever is in the source data; record it via
|
n_total |
Column name for total individuals per replicate.
Required for count responses; leave |
n_surv |
Column name for survivor counts. |
n_dead |
Column name for death counts. Converted to |
survival |
Column name for survival proportions in |
mortality |
Column name for mortality proportions in |
proportion |
Column name for a continuous proportion response in
|
proportion_eps |
Boundary clamp applied to |
random_effects |
Optional character vector of grouping variables for
random effects, e.g. |
duration_unit |
Label for the unit of |
temp_mean |
Value to subtract from |
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
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.), retain the grouping column. Use fit_4pl(by = "group") or grouped
formulas in tls_bf() to estimate category-level effects, or filter first
when separate models are scientifically preferable.
A tibble with the standardised columns plus a "tdt_meta" attribute
storing temp_mean, duration_unit, random_effects,
response_type ("count" or "proportion"), and response_var
(the response column name for a proportion fit, else NULL).
# 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")# 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")
tdt_parameter_table)Returns the fitted 4PL parameters (low, up, k, CTmax, z, and phi
for over-dispersed families) as point estimates with confidence intervals, in
bayesTLS's parameter / [group] / median / lower / upper shape.
tdt_parameter_table(object, method = NULL, level = 0.95)tdt_parameter_table(object, method = NULL, level = 0.95)
object |
A |
method |
Interval method: |
level |
Confidence level (default 0.95). |
A tibble with parameter, group, median, lower, upper.
raw <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) dat <- standardize_data( raw, temp = "temp", duration = "duration", n_total = "total", n_surv = "survived" ) fit <- fit_4pl(dat, family = "binomial", t_ref = 1, quiet = TRUE) tdt_parameter_table(fit, method = "wald")raw <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) dat <- standardize_data( raw, temp = "temp", duration = "duration", n_total = "total", n_surv = "survived" ) fit <- fit_4pl(dat, family = "binomial", t_ref = 1, quiet = TRUE) tdt_parameter_table(fit, method = "wald")
Quantile wrapper with TDT-friendly defaults
tdt_quantile(x, probs = c(0.025, 0.5, 0.975))tdt_quantile(x, probs = c(0.025, 0.5, 0.975))
x |
Numeric vector. |
probs |
Numeric vector of quantile probabilities. |
Numeric vector of length length(probs).
tdt_quantile(rnorm(100))tdt_quantile(rnorm(100))
Twins of the bayesTLS get_*_summary / get_*_draws accessors. *_summary
returns the median + interval tibble; *_draws returns the per-replicate
(bootstrap) tibble — the frequentist analogue of posterior draws.
get_z_summary(et) get_z_draws(et) get_ctmax_summary(et) get_ctmax_draws(et) get_tcrit_summary(et) get_tcrit_draws(et)get_z_summary(et) get_z_draws(et) get_ctmax_summary(et) get_ctmax_draws(et) get_tcrit_summary(et) get_tcrit_draws(et)
et |
An |
A tibble (see extract_tdt() for the column contract).
raw <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) dat <- standardize_data( raw, temp = "temp", duration = "duration", n_total = "total", n_surv = "survived" ) fit <- fit_4pl(dat, family = "binomial", t_ref = 1, quiet = TRUE) tdt <- extract_tdt(fit, nboot = 10, seed = 1) get_z_summary(tdt) get_ctmax_draws(tdt)raw <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) dat <- standardize_data( raw, temp = "temp", duration = "duration", n_total = "total", n_surv = "survived" ) fit <- fit_4pl(dat, family = "binomial", t_ref = 1, quiet = TRUE) tdt <- extract_tdt(fit, nboot = 10, seed = 1) get_z_summary(tdt) get_ctmax_draws(tdt)
tidy_parameters() returns a broom-style tibble of the natural-scale
parameter estimates with optional Wald confidence intervals. The intervals
are computed on the internal (unconstrained / link) scale as
estimate +/- z * std.error and then back-transformed to the natural scale,
so they respect each parameter's bounds (for example z > 0, 0 < low < up)
and are equivariant under the link. For CTmax (identity link) this is the
usual symmetric Wald interval.
tidy_parameters( fit, conf.int = TRUE, conf.level = 0.95, method = c("wald", "profile") )tidy_parameters( fit, conf.int = TRUE, conf.level = 0.95, method = c("wald", "profile") )
fit |
A |
conf.int |
Logical; include |
conf.level |
Confidence level for the interval (default |
method |
Either |
With method = "profile" the intervals are profile-likelihood confidence
intervals (see confint.profile_tls()); with method = "wald" (default) they
are the back-transformed internal-link Wald intervals. The returned shape is
identical; only interval_type and the interval values differ. A profile that
does not close returns NA on the open side (never a fabricated bound).
A tibble with one row per natural-scale parameter
and the columns parameter, group, estimate, std.error,
conf.low, conf.high, interval_type, and scale. scale is the link
on which the interval was constructed ("identity", "log", or
"logit"); interval_type is "wald" or "profile".
d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) tidy_parameters(fit) tidy_parameters(fit, method = "profile")d <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) fit <- fit_tls(d, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1) tidy_parameters(fit) tidy_parameters(fit, method = "profile")
The frequentist twin of bayesTLS::tls(). Reads a fit_4pl() (freq_tls)
fit and returns the headline thermal-death-time quantities — thermal
sensitivity z and CTmax — as point estimates with confidence intervals,
one row per group when the fit is grouped. Uncertainty uses the engine's
profile-likelihood intervals by default (or Wald / bootstrap via method).
tls( object, by = NULL, params = c("all", "z", "ctmax"), target_surv = "relative", lethal = FALSE, method = NULL, level = 0.95, nboot = 1000L, TC_rate_range = c(0.1, 1), seed = NULL ) tls_z(object, ...) tls_ctmax(object, ...) tls_tcrit(object, ...)tls( object, by = NULL, params = c("all", "z", "ctmax"), target_surv = "relative", lethal = FALSE, method = NULL, level = 0.95, nboot = 1000L, TC_rate_range = c(0.1, 1), seed = NULL ) tls_z(object, ...) tls_ctmax(object, ...) tls_tcrit(object, ...)
object |
A |
by |
Optional name for the grouping column in |
params |
|
target_surv |
Survival threshold for CTmax: |
lethal |
If |
method |
Interval method for the relative path: |
level |
Confidence level (default 0.95). |
nboot, TC_rate_range, seed
|
Passed to |
... |
Passed from |
A tls object: a list with $summary (a tibble of
[<group>,] quantity, median, lower, upper) and $meta.
fit_4pl(), tls_z(), tls_ctmax(), confint.profile_tls()
raw <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) dat <- standardize_data( raw, temp = "temp", duration = "duration", n_total = "total", n_surv = "survived" ) fit <- fit_4pl( dat, family = "binomial", t_ref = 1, method = "wald", quiet = TRUE ) tls(fit) tls_z(fit) tls_ctmax(fit)raw <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1) dat <- standardize_data( raw, temp = "temp", duration = "duration", n_total = "total", n_surv = "survived" ) fit <- fit_4pl( dat, family = "binomial", t_ref = 1, method = "wald", quiet = TRUE ) tls(fit) tls_z(fit) tls_ctmax(fit)
tls_bf() captures the per-sub-parameter formulas that define a freqTLS
model and returns them, unevaluated, as a tls_formula object. It is the
formula complement to the column interface of fit_tls(): instead of passing
bare column names, you write one response formula plus a formula per model
sub-parameter.
tls_bf(...)tls_bf(...)
... |
The response-and-axes formula first (unnamed, with a left-hand side), then sub-parameter formulas keyed by their left-hand side. See the grammar above. |
The first argument must be the unnamed response-and-axes formula. Its
left-hand side names the survival counts, in either the brms idiom
successes | trials(total) or the glm idiom cbind(successes, failures).
Its right-hand side names the two thermal-load-sensitivity axes with the
tagged markers time(<duration>) and temp(<temperature>) (order does not
matter).
The remaining arguments are sub-parameter formulas keyed by their left-hand
side, one of low, up, log_k, CTmax, or log_z. Any sub-parameter you
omit defaults to ~ 1. Each of low, up, and log_k may carry its
own design independently — a grouping factor (low ~ group), a
continuous covariate (log_k ~ body_size), or an intercept — and need not
share one factor or match the headline-parameter grouping. CTmax and log_z
accept fixed-effect formulas but must produce the same model-matrix columns
(for example, use CTmax ~ group, log_z ~ group); their supported random-
intercept groupings may differ. A single random intercept,
<param> ~ <fixed> + (1 | group), is accepted on CTmax, log_z, low, and
log_k (one grouping factor each, intercept only) – but not on the upper
asymptote up, for which the compiled objective has no random-intercept term. Putting the same
grouping factor on two or more of them fits independent variances (no
correlation term) and warns.
A tls_formula object: a list with the captured response formula,
the named sub-parameter formulas, and the calling environment.
The shape of the parser (variadic capture via substitute(), a per-entry
formula walk, and random-bar detection) is adapted from drmTMB's
drm_formula() / parse_drm_formula_entry() (GPL-3); see inst/COPYRIGHTS.
freqTLS writes its own grammar (the time() / temp() axis markers, the
five fixed sub-parameter handles, and the package's supported random-effect
grammar).
fit_tls(), which accepts either a tls_formula or the column
interface.
tls_bf( survived | trials(total) ~ time(duration) + temp(temp), CTmax ~ life_stage, log_z ~ life_stage ) # cbind() response idiom, ungrouped: tls_bf(cbind(survived, died) ~ time(duration) + temp(temp))tls_bf( survived | trials(total) ~ time(duration) + temp(temp), CTmax ~ life_stage, log_z ~ life_stage ) # cbind() response idiom, ungrouped: tls_bf(cbind(survived, died) ~ time(duration) + temp(temp))
binomial_tls() and beta_binomial_tls() describe the count response
distribution for fit_tls(), and beta_tls() the continuous-proportion
response in (0, 1) (e.g. PSII operating efficiency or relative chlorophyll
fluorescence). All three model survival as a four-parameter logistic function
of log10 duration; the beta-binomial and beta families add a dispersion
parameter phi.
binomial_tls() beta_binomial_tls() beta_tls()binomial_tls() beta_binomial_tls() beta_tls()
The phi convention for the beta-binomial family is the sum of the Beta
shape parameters: for fitted survival probability p, counts are
Beta-Binomial with shapes a = p * phi and b = (1 - p) * phi. Larger phi
means less overdispersion (the binomial is recovered as phi grows). This
matches the simulation convention in simulate_tls() and differs from the
precision/size parameterisations used by some other packages. The beta
family uses the same shapes for the continuous proportion, y ~ Beta(p * phi, (1 - p) * phi), so phi carries the identical meaning and a
larger phi again means a tighter response around the fitted curve.
A tls_family object: a list with family, family_code
(0 binomial, 1 beta-binomial, 2 beta), and links for the natural-scale
parameters.
binomial_tls() beta_binomial_tls() beta_tls()binomial_tls() beta_binomial_tls() beta_tls()
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.
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 )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 )
stage2 |
Output of |
method |
|
level |
Confidence level. Default 0.95. |
t_ref, time_multiplier, TC_rate_range
|
As in |
temp_grid |
Temperatures for the line band ( |
n_sim |
MVN draws ( |
seed |
RNG seed ( |
For "delta", a list with z and CTmax_1hr, each
list(point, lower, upper, lower_t, upper_t, se), plus df_resid.
For "mvn", a list with summary_ci (z/CTmax/T_crit bounds) and
curve_ci (per-temp_grid line band).
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")$zd <- 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
Median LT-vs-temperature line from a two-stage fit
ts_curve(stage2, temp_grid, time_multiplier = 1)ts_curve(stage2, temp_grid, time_multiplier = 1)
stage2 |
Output of |
temp_grid |
Temperatures (°C) to evaluate. |
time_multiplier |
Multiplier to minutes. Default 1. |
A tibble with temp and duration_median (minutes).
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))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))
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).
ts_stage1( data, temp = "temp", duration = "duration", n_surv = "n_surv", n_total = "n_total", family = c("binomial", "betabinomial") )ts_stage1( data, temp = "temp", duration = "duration", n_surv = "n_surv", n_total = "n_total", family = c("binomial", "betabinomial") )
data |
Data frame with one row per (temperature, duration) replicate. |
temp, duration, n_surv, n_total
|
Column names (strings) for assay temperature (°C), exposure duration, survivors, and trials. |
family |
|
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.
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.
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")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")
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 + z * mean(log10(TC_rate_range/100)).
ts_stage2( stage1, t_ref = 60, time_multiplier = 1, TC_rate_range = c(0.1, 1), rows = c("stage1_ok", "finite_ok") )ts_stage2( stage1, t_ref = 60, time_multiplier = 1, TC_rate_range = c(0.1, 1), rows = c("stage1_ok", "finite_ok") )
stage1 |
Output of |
t_ref |
Reference exposure duration for CTmax (minutes). Default 60. |
time_multiplier |
Multiplier from the Stage-1 duration unit to minutes (e.g. 60 if durations are in hours). Default 1. |
TC_rate_range |
Length-2 HI-rate range (% per hour) for T_crit. |
rows |
Which Stage-1 rows to keep: |
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.
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)$summaryd <- 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
Lethal-TDT trials for zebrafish (Danio rerio) at three life stages. Built from the raw daily survival sheet by summing the per-day morning/afternoon mortality counts into one death count per trial and dropping excluded rows. One row per assay trial. The model-ready frame for Case Study 2.
zebrafish_lethalzebrafish_lethal
A data frame with 323 rows and 7 variables:
Assay temperature (degrees C).
Exposure duration (hours).
Number of individuals in the trial.
Number that survived.
Number that died (n_total - n_surv).
Life stage, a factor with levels young_embryos,
old_embryos, larvae.
Experiment date (grouping factor).
Zebrafish lethal-TDT assay across life stages (Case Study 2),
obtained from the bayesTLS package distribution by Noble, Arnold, and
Pottier (2026), licensed CC BY 4.0. freqTLS removed excluded trials,
aggregated daily mortality counts, and derived survivors as documented
above. Raw file:
system.file("extdata", "data_lethal_TDT_zebrafish.csv", package = "freqTLS").
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 profile-likelihood
confidence intervals on those direct parameters.
zebrafish_o2zebrafish_o2
A data frame with 905 rows and 10 variables:
Larval cohort identifier.
Ploidy, a factor with levels diploid, triploid.
Oxygen treatment, a factor with levels hypoxia,
normoxia, hyperoxia (the modelling moderator).
Nominal oxygen target as percent air saturation (25 / 100 / 225).
Measured oxygen level (percent air saturation).
Target assay temperature (degrees C).
Measured assay temperature (degrees C).
Exposure duration (minutes).
Number of larvae in the assay group.
Number surviving after exposure and recovery.
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 = "freqTLS").
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 treatmentstd <- 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