--- title: "Heat injury under a temperature trace" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Heat injury under a temperature trace} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 3.6, dpi = 96 ) ``` This vignette is for an applied thermal biologist who has fitted a thermal-load-sensitivity (TLS) curve and now wants to ask a forward question: if an organism experiences *this* temperature time-series, how much lethal damage accumulates, and how does survival fall? It mirrors the heat-injury panel of the `bayesTLS` manuscript (their "Heat injury and survival" figure) but is computed entirely from the `freqTLS` maximum-likelihood fit — no Stan, no Bayesian comparison, because heat injury here is a **deterministic prediction off the fitted curve**, not a new model. ```{r setup} library(freqTLS) library(ggplot2) ``` ## What heat-injury prediction is (and is not) `predict_heat_injury()` accumulates thermal damage along a temperature trace and reads survival back off the already-fitted 4PL. At each temperature `T` the fitted curve gives a lethal time `LT(T)`; the instantaneous damage rate is its reciprocal, `1 / LT(T)`, in **lethal doses per time unit**. Cumulative dose is the forward-Euler integral of that rate over the trace, and survival is the 4PL evaluated at the accumulated dose. One **lethal dose** (`dose = 1`, or 100% on the injury axis) is the thermal load that drives survival to a chosen target — by default the project-default **relative** threshold, the curve midpoint `(low + up) / 2`. This is the maximum-likelihood analogue of `bayesTLS::predict_heat_injury()`. The boundary is deliberate and matches the sibling-package split: * `freqTLS` **predicts** injury from a curve it has already fitted. * Fitting injury or repair *dynamics* as a model — estimating damage and recovery rates from data — is a `bayesTLS` concern. `freqTLS` does not fit them. Everything below is a transform of one fitted `CTmax` / `z` (plus the shape `low`, `up`, `k`), so the uncertainty in the prediction is exactly the confidence-interval uncertainty in those estimates — which we propagate with a prior-free bootstrap, never a posterior. ## Fit the curve We use the vendored brown-shrimp lethal assay (`shrimp_lethal`), the same dataset as the `bayesTLS` shrimp case study: a beta-binomial 4PL with constant shape and a one-hour reference time. ```{r fit} data(shrimp_lethal) shrimp_std <- standardize_data( shrimp_lethal, temp = "Temperature_assay", duration = "Duration_exposure_hours", n_total = "N_individuals_after_trial", mortality = "Mortality_after_trial", duration_unit = "hours" ) # The heat-injury helpers take the engine fit; fit_4pl()$fit is that engine fit. fit <- fit_4pl(shrimp_std, t_ref = 1, family = "beta_binomial", quiet = TRUE)$fit ctmax <- get_ctmax(fit)$estimate z <- get_z(fit)$estimate up <- tidy_parameters(fit)$estimate[tidy_parameters(fit)$parameter == "up"] c(CTmax = round(ctmax, 2), z = round(z, 2), upper_asymptote = round(up, 3)) ``` The fitted survival ceiling is `up` $\approx$ `r round(up, 2)`: even brief benign exposures do not return survival to a full 100%. Heat-injury trajectories will therefore **start at that fitted ceiling**, not at 1 — an honesty point the `bayesTLS` figure also respects. The rate-multiplier critical temperature `T_crit` — the temperature below which the damage rate is treated as negligible — follows directly from `CTmax` and `z`: ```{r tcrit} tcrit <- derive_tcrit(fit, rate = 1) # 1% of a lethal dose per hour round(tcrit, 2) ``` So damage accrues meaningfully only above `r round(tcrit, 1)` °C. We will pass `t_c = tcrit` to `predict_heat_injury()` as the damage cutoff and draw it as a reference line. ## A temperature trace We construct a realistic five-day diel heatwave: a sinusoidal day–night cycle whose afternoon peaks rise into the mid-heatwave and then ease. The trace is a plain `data.frame(time, temp)` with `time` in **hours** — the same unit as the fit's `tref`, which is required (the damage rate is per that unit). Like the `bayesTLS` shrimp panel, which used a warming projection because the literal fjord never reaches the shrimp's `T_crit`, this trace is a near-future scenario: its peaks sit above `T_crit` but stay **below** the fitted `CTmax`. ```{r trace} hours <- seq(0, 5 * 24, by = 0.25) # quarter-hourly, 5 days day_peak <- approx(c(0:5) * 24, c(25, 27, 29, 30, 29.5, 28), xout = hours)$y night_min <- approx(c(0:5) * 24, c(20, 21, 22.5, 23.5, 23, 21.5), xout = hours)$y mid_t <- (day_peak + night_min) / 2 amp_t <- (day_peak - night_min) / 2 trace <- data.frame(time = hours, temp = mid_t - amp_t * cos(2 * pi * (hours %% 24) / 24)) c(max_temp = round(max(trace$temp), 2), hours_above_Tcrit = round(sum(trace$temp > tcrit) * 0.25, 1)) ``` The trace peaks at `r round(max(trace$temp), 1)` °C and spends about `r round(sum(trace$temp > tcrit) * 0.25)` hours above `T_crit` over the five days. ## Predict injury and survival One call accumulates the dose and reads survival off the curve. We pass the `T_crit` cutoff so sub-`T_crit` temperatures contribute no damage. ```{r predict} hi <- predict_heat_injury(fit, trace, t_c = tcrit) str(hi, give.attr = FALSE) final <- hi[nrow(hi), c("dose", "injury", "survival")] round(final, 3) ``` Under this trace the cumulative dose reaches **`r round(final$dose, 2)` lethal doses** (`r round(final$injury)`% on the injury axis) and predicted survival falls from the fitted ceiling of `r round(hi$survival[1], 2)` to **`r round(final$survival, 2)`**. Crossing `dose = 1` is the moment the accumulated load equals one lethal dose; survival at that point is the relative midpoint by construction. The figure below stacks the three quantities the `bayesTLS` panel shows: the temperature trace (with `CTmax` and `T_crit` reference lines), cumulative injury (100% = one lethal dose), and predicted survival. ```{r plot-trace, fig.height = 2.6, fig.alt = "Diel temperature trace over five days: a daily cycle oscillating between about 20 and 30 degrees C, with a dashed line at the fitted CTmax (31.8 C, above all peaks) and a dotted line at T_crit (27.4 C, crossed by the afternoon peaks)."} panel_trace <- ggplot(hi, aes(time / 24, temp)) + geom_hline(yintercept = ctmax, linetype = "dashed", colour = "grey40") + geom_hline(yintercept = tcrit, linetype = "dotted", colour = "grey40") + geom_line(colour = "#B2182B", linewidth = 0.6) + annotate("text", x = 0.1, y = ctmax, label = "CTmax", vjust = -0.4, hjust = 0, size = 3, colour = "grey30") + annotate("text", x = 0.1, y = tcrit, label = "T_crit", vjust = 1.3, hjust = 0, size = 3, colour = "grey30") + labs(x = NULL, y = "Temp (°C)") + theme_bw(base_size = 11) print(panel_trace) ``` ```{r plot-injury, fig.height = 2.6, fig.alt = "Cumulative heat injury over five days: a stepwise rise in injury percentage, crossing the 100%-of-one-lethal-dose dashed reference line near day three and reaching about 160% by day five."} panel_injury <- ggplot(hi, aes(time / 24, injury)) + geom_hline(yintercept = 100, linetype = "dashed", colour = "grey40") + geom_line(colour = "#762A83", linewidth = 0.6) + annotate("text", x = 0.1, y = 100, label = "one lethal dose", vjust = -0.4, hjust = 0, size = 3, colour = "grey30") + labs(x = NULL, y = "Injury (%)") + theme_bw(base_size = 11) print(panel_injury) ``` ```{r plot-survival, fig.height = 2.6, fig.alt = "Predicted survival over five days: declining in steps from the fitted ceiling near 0.94, tracking each afternoon heating episode, to about 0.29 by day five."} panel_surv <- ggplot(hi, aes(time / 24, survival)) + geom_line(colour = "#1B7837", linewidth = 0.6) + labs(x = "Day", y = "Survival") + ylim(0, 1) + theme_bw(base_size = 11) print(panel_surv) ``` (The three panels share the day axis; they are printed separately to avoid a layout dependency. The survival curve is monotone because `predict_heat_injury()` defaults to `irreversible = TRUE`: mortality does not reverse even when the dose rate drops at night.) ## A repair scenario (a what-if, not an estimate) Optionally, `predict_heat_injury()` can subtract a Sharpe–Schoolfield repair rate each step, modelling sub-lethal recovery between heat exposures. **These repair parameters are not identified by the survival data the curve was fitted to.** They are a user-supplied scenario layer — a what-if — and the function warns to say so. We surface that warning rather than hide it. ```{r repair} # Illustrative Sharpe-Schoolfield kernel (reference temperatures in Kelvin). # NOT fitted: a scenario, in the same spirit as the bayesTLS illustration. repair <- list( r_ref = 0.03, t_a = 14065, t_al = 50000, t_ah = 120000, t_l = 10.5 + 273.15, t_h = 22.5 + 273.15, t_ref = 17 + 273.15 ) hi_repair <- predict_heat_injury(fit, trace, t_c = tcrit, repair = repair) round(hi_repair[nrow(hi_repair), c("dose", "survival")], 3) ``` With this illustrative repair the final dose drops from `r round(final$dose, 2)` to `r round(hi_repair$dose[nrow(hi_repair)], 2)` lethal doses and final survival rises from `r round(final$survival, 2)` to `r round(hi_repair$survival[nrow(hi_repair)], 2)`. The repaired trajectory is strictly a sensitivity scenario: a different, equally-unidentified repair kernel would move it elsewhere. Read it as "what if recovery looked like this", never as an estimate. ```{r plot-repair, fig.alt = "Predicted survival over five days for the same trace, with and without the illustrative repair kernel. The no-repair curve (solid) declines to about 0.29; the with-repair curve (dashed) declines less steeply to about 0.40. A caption notes the repair kernel is a non-identified what-if."} surv_compare <- rbind( data.frame(time = hi$time, survival = hi$survival, scenario = "no repair"), data.frame(time = hi_repair$time, survival = hi_repair$survival, scenario = "with repair (what-if)") ) ggplot(surv_compare, aes(time / 24, survival, linetype = scenario)) + geom_line(colour = "#1B7837", linewidth = 0.6) + scale_linetype_manual(values = c("no repair" = "solid", "with repair (what-if)" = "longdash")) + labs(x = "Day", y = "Survival", linetype = NULL, caption = "Repair kernel is a user-supplied scenario, not identified by the data.") + ylim(0, 1) + theme_bw(base_size = 11) + theme(legend.position = "top") ``` ## A bootstrap envelope on survival A single deterministic trajectory hides the uncertainty in the fitted `CTmax` and `z`. The honest frequentist counterpart of the `bayesTLS` posterior band is a **parametric bootstrap**: redraw the curve parameters from the fitted model, re-predict the trace for each draw, and shade the spread of survival trajectories. This carries no prior and makes no probability statement about the parameters — it is a confidence band, not a credible band. `heat_injury_envelope()` does this for us: it redraws the curve parameters by parametric bootstrap (the same machinery `confint()` uses) and re-integrates each draw through the **same** dose-accumulation map `predict_heat_injury()` documents (lethal time $LT(T) = t_{ref}\,10^{(CTmax - T)/z}$, damage rate $1/LT$, survival $= low + (up - low)\,\mathrm{logit}^{-1}(-k\log_{10}\!D)$), so the band is consistent with the point trajectory by construction. It returns a tidy `time` / `temp` / `survival` / `conf.low` / `conf.high` table. ```{r bootstrap} # Prior-free parametric-bootstrap confidence band on survival. env <- heat_injury_envelope(fit, trace, t_c = tcrit, nboot = 300L, seed = 1) round(env[nrow(env), c("survival", "conf.low", "conf.high")], 3) ``` ```{r bootstrap-plot, fig.alt = "Predicted survival over five days as a dark green point-estimate line inside a pale green shaded 95% bootstrap confidence band that widens as survival declines."} plot_heat_injury(fit, trace, t_c = tcrit, nboot = 300L, seed = 1, time_div = 24, xlab = "Day") ``` At the end of the trace the point estimate is `r round(env$survival[nrow(env)], 2)`, with a 95% bootstrap confidence band of **[`r round(env$conf.low[nrow(env)], 2)`, `r round(env$conf.high[nrow(env)], 2)`]**. The band widens as survival falls because the damage rate $1/LT(T) = 10^{(T - CTmax)/z}/t_{ref}$ is exponential in temperature: a one-degree shift in `CTmax` multiplies the hourly dose by $10^{1/z}$, so small parameter uncertainty fans out into a wide survival envelope late in a heatwave. That is the same exponential-sensitivity caveat the `bayesTLS` figure flags — here it falls out of the likelihood, not a posterior. ## Absolute versus relative: what "one lethal dose" means By default, `dose = 1` is the **relative** threshold — the load that drives survival to the curve midpoint `(low + up) / 2`. Supplying `target_surv` redefines one lethal dose as an **absolute** survival level. With a fitted ceiling of `r round(up, 2)`, the relative midpoint sits near `r round((up + 0) / 2, 2)`, so asking for an absolute 50% survival is a slightly stricter target and shifts where `dose = 1` lands. ```{r target} hi_abs <- predict_heat_injury(fit, trace, t_c = tcrit, target_surv = 0.5) data.frame( definition = c("relative midpoint (default)", "absolute 50% survival"), final_dose = round(c(final$dose, hi_abs$dose[nrow(hi_abs)]), 3), final_surv = round(c(final$survival, hi_abs$survival[nrow(hi_abs)]), 3) ) ``` The trajectories are the same physical exposure; only the *yardstick* changes. The absolute target rescales the dose axis (the final dose moves from `r round(final$dose, 2)` to `r round(hi_abs$dose[nrow(hi_abs)], 2)` lethal doses) without changing the underlying survival prediction much, because both thresholds read off the same fitted curve. Use the relative threshold to stay consistent with the benchmark configuration; use an absolute `target_surv` when a specific survival level (a regulatory or ecological floor) defines "lethal" for your question. ## Extensions and boundaries * **Group-specific traces.** For a grouped fit, pass `group = ""` to predict under one group's curve; loop over levels for a multi-panel figure. * **Irregular traces.** The integrator uses each step's real time increment, so a raw logger series with uneven sampling integrates correctly — no need to resample to a fixed `dt`. * **What stays in `bayesTLS`.** Estimating injury and repair *rates* from data — fitting a repair kernel rather than supplying one — is outside `freqTLS`. The repair layer here is illustrative only, and `predict_heat_injury()` says so on every call. All numbers on this page are computed live as it renders, from the fitted `freqTLS` curve, with no Stan and no posterior. The uncertainty is a prior-free bootstrap confidence band; we use "confidence" language throughout and never describe a `freqTLS` interval as a posterior or credible interval.