---
title: "Temperature and a co-stressor"
vignette: >
  %\VignetteIndexEntry{Temperature and a co-stressor}
  %\VignetteEngine{quarto::html}
  %\VignetteEncoding{UTF-8}
format:
  html:
    toc: true
    toc-depth: 3
    toc-title: "Contents"
    html-math-method: mathml
    include-in-header:
      text: |
        <style>
        /* The contents list beside the text, in view while scrolling. It applies
           only to the plain page the package build makes, on a window wide enough
           to hold it; a themed Quarto page places its own contents list. */
        @media (min-width: 1380px) {
          body:not(:has(#quarto-content)) nav#TOC {
            position: fixed;
            top: 2em;
            left: calc(50% - 670px);
            width: 250px;
            max-height: 90vh;
            overflow-y: auto;
            font-size: 0.9em;
          }
          body:not(:has(#quarto-content)) nav#TOC ul { padding-left: 1.1em; }
        }
        /* An estimate and its interval stay on one line of a table; only the first
           column, which holds names, may wrap. */
        table td:not(:first-child), table th:not(:first-child) { white-space: nowrap; }
        </style>
knitr:
  opts_chunk:
    collapse: true
    comment: "#>"
    fig.width: 7
    fig.height: 4
---

```{r}
#| label: setup
#| include: false
results <- readRDS("results-04.rds")
models  <- names(results$models)
main    <- results$models[["Main effect"]]
inter   <- results$models[["Main effect and interaction"]]

# A value as a percentage, and an interval as "lower–upper", for the text.
percent  <- function(x) sprintf("%.0f", 100 * x)
interval <- function(x, digits = 1, scale = 1) {
  sprintf(paste0("%.", digits, "f–%.", digits, "f"), scale * x$lower, scale * x$upper)
}
# An estimate as "median (lower–upper)", kept on one line of a table.
estimate <- function(x, digits = 1, scale = 1) {
  f <- paste0("%.", digits, "f (%.", digits, "f–%.", digits, "f)")
  sprintf(f, scale * x$median, scale * x$lower, scale * x$upper)
}
```

Temperature does not always act alone. For example, fish exposed to high temperatures may at the same time encounter low dissolved oxygen, and below a species-specific oxygen level their upper thermal limit declines ([Ern et al. 2016](https://doi.org/10.1242/jeb.143495){target="_blank"}). In `rTFT`, a second environmental variable of this kind is declared as a co-stressor, and the failure rate then depends on both temperature and the co-stressor.

In this vignette, we fit two models to a ramping assay that was run at two levels of dissolved oxygen: one in which oxygen changes the failure rate by the same proportion at every temperature, and one in which its effect also changes with temperature (i.e., an interaction). We then: (i) derive the critical thermal maxima, the thermal death time curves and the thermal tolerance landscape at chosen levels of oxygen, and (ii) predict survival along a record of stream temperature and dissolved oxygen. We also show what an interaction between temperature and a co-stressor implies when a model is used to predict at temperatures below those of the assay. The steps shared with the previous vignettes are described in [*Getting started with rTFT*](getting-started.html) and [*Ramping assays*](ramping-assays.html).

```{r}
#| label: libraries
#| message: false
library(rTFT)
library(ggplot2)
```

## The assay

We use the 240 fish of `ramp_assay` that were warmed at 0.025, 0.1 or 0.3°C per min, which are the rates with the wider span in the vignette [*Ramping assays*](ramping-assays.html). Each rate was run in four trials of 10 fish in normoxia (10 mg/L of dissolved oxygen) and four in hypoxia (5 mg/L), so the rates and the levels of oxygen are fully crossed. The oxygen of each trial is in the column `do`.

```{r}
#| label: data-summary
rates <- c(0.025, 0.1, 0.3)
ramp  <- subset(ramp_assay, rate %in% rates)
ramp$ctmax <- 14 + ramp$rate * ramp$t_event
observed <- ramp |>
  dplyr::summarise(
    n_fish       = dplyr::n(),
    n_failed     = sum(event),
    median_ctmax = median(ctmax),
    .by = c(rate, do)
  ) |>
  dplyr::arrange(rate, dplyr::desc(do))
observed
```

```{r}
#| label: data-text
#| include: false
drop_in_ctmax <- observed$median_ctmax[observed$do == 10] - observed$median_ctmax[observed$do == 5]
```

At each of the three rates, the median CTmax was lower in hypoxia than in normoxia, by `r sprintf("%.1f", min(drop_in_ctmax))`–`r sprintf("%.1f", max(drop_in_ctmax))`°C.

::: {.callout-note}
The data were simulated (see `?ramp_assay`). In the model used to simulate them, low oxygen raises the failure rate, and the more so the warmer the water: each mg/L of oxygen below 10 mg/L adds 0.30 to the logarithm of the failure rate at 26°C, and this effect changes by 0.04 for each degree (i.e., it is 0.22 at 24°C and 0.38 at 28°C). The effect of oxygen is therefore a straight line in temperature, which is the form of the interaction fitted in this vignette, and it reaches zero at 18.5°C.
:::

## Dissolved oxygen as a co-stressor

A co-stressor is declared with the function `costressor()`, which requires three choices that have no safe default: the reference level at which the co-stressor is known or expected to add no stress (`ref`), the unit in which its effect is expressed (`scale`) and whether more or less of it is more stressful (`direction`). For dissolved oxygen we take normoxia as the reference (10 mg/L), express the effect per 1 mg/L, and set the direction to −1 because less oxygen is more stressful.

```{r}
#| label: costressor
DO <- costressor("DO", column = "do", ref = 10, scale = 1, direction = -1)
```

With these choices the model uses the oxygen deficit, $z = 10 - \mathrm{DO}$ (i.e., 0 in normoxia and 5 at 5 mg/L; this $z$ is the package's name for a standardised co-stressor, and not the $z$ of the thermal death time curves), and the failure rate becomes

$$
\log M(T, z) = \log a + b (T - T_r) + c \, z .
$$

At the reference ($z = 0$) the model is the one of the previous vignettes, so $a$ and $b$ keep their meaning and describe a fish in normoxia. The main effect $c$ (the parameter `cDO`) is the change in the logarithm of the failure rate for each mg/L of oxygen deficit, so each mg/L of deficit multiplies the failure rate by $\exp(c)$, at every temperature.

A co-stressor may also interact with temperature, which is asked for with `interaction = TRUE`.

```{r}
#| label: costressor-interaction
DO_interaction <- costressor("DO", column = "do", ref = 10, scale = 1, direction = -1,
                             interaction = TRUE)
```

The failure rate then has a further term,

$$
\log M(T, z) = \log a + b (T - T_r) + c \, z + d \, (T - T_r) \, z ,
$$

in which the interaction $d$ (the parameter `dTDO`) is the change in the effect of oxygen for each degree above the reference temperature $T_r$. With $d > 0$ the effect of low oxygen is stronger at higher temperatures, and $c$ is its effect at the reference temperature.

## The table of intervals at risk

The co-stressor is given to `build_hazard_data()`, which adds the column `z_DO` with the oxygen deficit of each row. Here oxygen was constant within a trial, so it is read from the column `do` of the fish. A co-stressor that changes during a trial is given as a table of its values through time, as temperature is (see the argument `trajectory` in `?costressor`).

```{r}
#| label: long
#| message: false
path <- ramp_path(c("r0.025", "r0.1", "r0.3"), T_start = 14, rate = rates,
                  t_max = 15.5 / rates)
long <- build_hazard_data(ramp, path, observation = "continuous",
                          covariates = c("mass", "trial"), costressors = list(DO))
unique(long[c("path_id", "z_DO")])
```

## Fitting

The same list of co-stressors is given to `fit_hazard()`. We fit the model of the previous vignettes (an exponential failure rate with a Weibull baseline, with $\log a$ depending on body mass and varying among trials) with oxygen entering in each of the two ways. The table of intervals at risk is the same for both, because it holds the oxygen deficit whichever way oxygen enters the model.

```{r}
#| label: fit
#| eval: false
formulas <- list(loga = ~ 1 + I(mass - 15) + (1 | trial))
fit_main <- fit_hazard(
  long,
  m_form        = "exp",
  baseline      = "weibull",
  dpar_formulas = formulas,
  costressors   = list(DO),
  chains = 4, iter = 2000, warmup = 1000, cores = 4, seed = 20261006,
  control = list(adapt_delta = 0.95), backend = "cmdstanr"
)
fit_interaction <- fit_hazard(
  long,
  m_form        = "exp",
  baseline      = "weibull",
  dpar_formulas = formulas,
  costressors   = list(DO_interaction),
  chains = 4, iter = 2000, warmup = 1000, cores = 4, seed = 20261006,
  control = list(adapt_delta = 0.95), backend = "cmdstanr"
)
```

The default priors of the parameters of the co-stressor are shown by `build_hazard_prior()` when it is given the same list (see `?build_hazard_prior`). The sampler settings are those chosen in the vignette [*Getting started with rTFT*](getting-started.html) and were not tuned for these models. Fitting the two models takes about an hour, so the code of this vignette that uses the fitted models was run once and its results were saved with the package (see the note at the end of this vignette).

```{r}
#| label: tbl-convergence
#| echo: false
#| tbl-cap: "Convergence of the two models (the largest $\\hat{R}$, the smallest bulk and tail effective sample sizes among the parameters and the number of divergent transitions) and their fit to the data (the number of observation times, among those of every ramp and level of oxygen, at which the proportion of fish observed to survive was outside the 95% band of the proportions simulated from the model, and the P-value of a Kolmogorov–Smirnov test of the residuals against the uniform distribution)."
checks <- results$fit_checks
knitr::kable(data.frame(
  Model = models,
  `Largest R-hat` = sprintf("%.3f", vapply(results$convergence, \(x) max(x$table$rhat), numeric(1))),
  `Smallest bulk ESS` = sprintf("%.0f", vapply(results$convergence, \(x) min(x$table$ess_bulk), numeric(1))),
  `Smallest tail ESS` = sprintf("%.0f", vapply(results$convergence, \(x) min(x$table$ess_tail), numeric(1))),
  `Divergent transitions` = vapply(results$convergence, \(x) x$n_divergent, numeric(1)),
  `Times outside the band` = sprintf("%d of %d", checks$n_outside, checks$n_times),
  `P (residuals)` = sprintf("%.2f", checks$ks_p),
  check.names = FALSE, row.names = NULL
), align = "lrrrrrr")
```

```{r}
#| label: checks-text
#| include: false
n_divergent <- vapply(results$convergence, \(x) x$n_divergent, numeric(1))
```

Both models converged`r if (all(n_divergent == 0)) ", with no divergent transitions" else ""` (@tbl-convergence). The two checks of fit were made as in the vignette [*Getting started with rTFT*](getting-started.html), with the proportion of fish surviving compared for each ramp at each level of oxygen (`check_survival(fit, by = "z_DO")`). The proportion observed was outside the simulated band at `r checks$n_outside[2]` of the `r format(checks$n_times[2], big.mark = ",")` observation times under the model with the interaction and at `r checks$n_outside[1]` under the model with the main effect alone; a difference of that size is no evidence that one model describes the data better than the other, and the two could be compared by cross-validation as in the vignette [*Comparing and averaging models*](comparing-and-averaging-models.html).

## The effect of oxygen

The parameters are obtained with `derive_tdt_ttl()`, which for a model with a co-stressor requires the level of the co-stressor at which the derived quantities are reported (`levels`), on its measurement scale.

```{r}
#| label: params-code
#| eval: false
tol <- derive_tdt_ttl(fit_interaction, newdata = data.frame(mass = 15), t_ref = 60,
                           levels = c(DO = 10), Ta = seq(25.5, 29.5, by = 0.1))
tol$params
```

```{r}
#| label: tbl-params
#| echo: false
#| tbl-cap: "Parameters of the two models for a fish of 15 g in normoxia (posterior medians and 95% credible intervals), with the values of the model used to simulate the data. The amplitude $a$ and the thermal sensitivity $b$ are on their natural scale (the models estimate their logarithms). The main effect of oxygen is per mg/L of oxygen deficit, and its interaction with temperature per mg/L and per °C. For the main effect, the simulated value is the effect of oxygen at the reference temperature of the fitted models."
truth   <- attr(ramp_assay, "truth")
Tr      <- unique(long$Tr)
tau_ref <- unique(long$tau_ref)
value <- function(model, quantity, digits = 2) {
  x <- model$params[model$params$quantity == quantity, ]
  if (!nrow(x)) return("")
  estimate(x, digits)
}
shown <- c(a = "a", b = "b (per °C)", shape = "k", cDO = "Main effect of oxygen",
           dTDO = "Interaction with temperature")
# the effect of oxygen in the model used to simulate the data, at the reference
# temperature of the fitted models
effect_at_Tr <- truth$b_do + truth$d_do * (Tr - truth$T_ref)
simulated <- c(
  # a at the reference temperature and time of the fitted models
  exp(truth$log_a + truth$b * (Tr - truth$T_ref)) * (tau_ref / truth$scale)^truth$shape,
  truth$b, truth$shape, effect_at_Tr, truth$d_do
)
knitr::kable(data.frame(
  Parameter = unname(shown),
  Simulated = sprintf("%.2f", simulated),
  `Main effect` = vapply(names(shown), \(q) value(main, q), character(1)),
  `Main effect and interaction` = vapply(names(shown), \(q) value(inter, q), character(1)),
  check.names = FALSE, row.names = NULL
), align = "lrrr")
```

```{r}
#| label: params-text
#| include: false
c_main  <- main$params[main$params$quantity == "cDO", ]
d_inter <- inter$params[inter$params$quantity == "dTDO", ]
factor_of <- function(x, deficit = 1) sprintf("%.1f", exp(deficit * x))
```

```{r}
#| label: params-text-2
#| include: false
c_inter <- inter$params[inter$params$quantity == "cDO", ]
```

In the model with the main effect alone, the effect of oxygen was `r sprintf("%.2f", c_main$median)` per mg/L (95% credible interval: `r interval(c_main, 2)`), which is on the scale of the logarithm of the failure rate. Its exponential is the factor by which each mg/L of oxygen deficit multiplied the failure rate: $\exp(`r sprintf("%.2f", c_main$median)`)$ = `r factor_of(c_main$median)` (`r factor_of(c_main$lower)`–`r factor_of(c_main$upper)`), at every temperature (@tbl-params). A deficit of several mg/L multiplies the failure rate by that factor once for each mg/L (e.g., by $\exp(5 \times `r sprintf("%.2f", c_main$median)`)$ = `r factor_of(c_main$median, 5)` for the 5 mg/L of deficit of the trials in hypoxia).

In the model with the interaction, the effect was `r sprintf("%.2f", c_inter$median)` per mg/L (`r interval(c_inter, 2)`) at the reference temperature (`r sprintf("%.1f", Tr)`°C) and increased by `r sprintf("%.2f", d_inter$median)` (`r interval(d_inter, 2)`) for each degree, with a posterior probability of `r sprintf("%.3f", results$p_interaction_positive)` that the interaction is positive. Low oxygen was therefore more harmful at higher temperatures, as it is in the model used to simulate the data. The estimate of the interaction was about twice the simulated value (`r sprintf("%.2f", truth$d_do)`), which is within its credible interval.

What each model implies for temperatures other than those of the assay is seen by computing the effect of oxygen on the failure rate across temperatures. The function `M_draws()` returns the posterior draws of the failure rate $M$ at given temperatures and levels of the co-stressor (the rate itself, and not its logarithm). The effect of 1 mg/L of oxygen deficit at a temperature is the difference between the logarithm of the failure rate at 9 mg/L and that at 10 mg/L (i.e., $\log M_9 - \log M_{10} = \log(M_9 / M_{10})$), which we compute for each posterior draw and summarise by its median and 95% credible interval.

```{r}
#| label: effect-code
#| eval: false
# The effect of 1 mg/L of oxygen deficit on the logarithm of the failure rate, at
# each temperature, from a fitted model.
oxygen_effect <- function(fit) {
  temperatures <- seq(16, 29, by = 0.5)
  fish     <- data.frame(mass = 15)
  normoxia <- M_draws(fit, Ta = temperatures, newdata = fish, levels = c(DO = 10),
                      n_draws = 1000)
  deficit  <- M_draws(fit, Ta = temperatures, newdata = fish, levels = c(DO = 9),
                      n_draws = 1000)
  normoxia$effect <- log(deficit$M / normoxia$M)
  normoxia |>
    dplyr::summarise(
      median = median(effect),
      lower  = quantile(effect, 0.025),
      upper  = quantile(effect, 0.975),
      .by = Ta
    )
}
effect <- rbind(
  data.frame(model = "Main effect", oxygen_effect(fit_main)),
  data.frame(model = "Main effect and interaction", oxygen_effect(fit_interaction))
)
```

```{r}
#| label: effect-load
#| include: false
effect <- rbind(
  data.frame(model = models[1], main$oxygen_effect),
  data.frame(model = models[2], inter$oxygen_effect)
)
```

```{r}
#| label: fig-effect
#| fig-cap: "Effect of 1 mg/L of oxygen deficit on the logarithm of the failure rate, across temperatures, under each of the two models. Lines are posterior medians and bands are 95% credible intervals. The dashed line is the effect in the model used to simulate the data, and the solid grey line marks no effect. The shaded area covers the temperatures at which fish failed in the assay."
#| fig-height: 3.4
effect$model <- factor(effect$model, levels = unique(effect$model))
# the temperatures at which fish failed in the assay
failed_at <- range(ramp$ctmax[ramp$event == 1])
# the effect in the model used to simulate the data
simulated_effect <- data.frame(Ta = seq(16, 29, by = 0.5))
simulated_effect$effect <- truth$b_do + truth$d_do * (simulated_effect$Ta - truth$T_ref)

ggplot(effect, aes(x = Ta)) +
  annotate("rect", xmin = failed_at[1], xmax = failed_at[2], ymin = -Inf, ymax = Inf,
           fill = "grey85", alpha = 0.5) +
  geom_hline(yintercept = 0, colour = "grey50") +
  geom_line(data = simulated_effect, aes(y = effect), linetype = "dashed", colour = "grey20") +
  geom_ribbon(aes(ymin = lower, ymax = upper), fill = "#225EA8", alpha = 0.2) +
  geom_line(aes(y = median), colour = "#225EA8") +
  facet_wrap(vars(model)) +
  labs(x = "Temperature (°C)", y = "Effect of 1 mg/L of oxygen deficit") +
  theme_light()
```

In @fig-effect, the horizontal axis is temperature and the vertical axis is what 1 mg/L of oxygen deficit adds to the logarithm of the failure rate at that temperature. A value above zero means that low oxygen raises the failure rate (e.g., a value of 0.4 multiplies it by $\exp(0.4)$ = 1.5), a value of zero that it has no effect, and a value below zero that it lowers the failure rate. Under the model with the main effect alone the effect is the same at every temperature, so its line is horizontal by construction. Under the model with the interaction it is a straight line whose slope is the interaction $d$.

```{r}
#| label: effect-text
#| include: false
at <- function(Ta) inter$oxygen_effect[inter$oxygen_effect$Ta == Ta, ]
crossing <- with(inter$oxygen_effect, max(Ta[median < 0]))
minus <- function(x, digits = 2) sub("-", "−", sprintf(paste0("%.", digits, "f"), x))
simulated_at <- function(Ta) truth$b_do + truth$d_do * (Ta - truth$T_ref)
within <- function(x, value) x$lower <= value & value <= x$upper
# what the two paragraphs below state about the estimates
stopifnot(
  within(at(28), simulated_at(28)), within(at(20), simulated_at(20)),
  within(d_inter, truth$d_do), d_inter$median > truth$d_do,
  at(20)$upper - at(20)$lower > 2 * (at(28)$upper - at(28)$lower),
  abs(c_main$median - effect_at_Tr) < 0.05
)
```

Within the temperatures at which most fish failed, the model with the interaction recovered the effect used to simulate the data: at 28°C it was `r sprintf("%.2f", at(28)$median)` (`r interval(at(28), 2)`), against a simulated `r sprintf("%.2f", simulated_at(28))` (@fig-effect). Below those temperatures, the credible band widens with the distance from them, because the interaction makes the effect of oxygen a straight line in temperature that is estimated where the fish failed and extended from there. In this assay the interaction was estimated to be steeper than the simulated one, so the fitted line crosses zero between `r sprintf("%.1f", crossing)` and `r sprintf("%.1f", crossing + 0.5)`°C, which the stream reaches during its heat wave (see below), whereas the simulated line does so at `r sprintf("%.1f", truth$T_ref - truth$b_do / truth$d_do)`°C. At 20°C the model gave an effect of `r minus(at(20)$median)` (95% credible interval: `r minus(at(20)$lower)` to `r minus(at(20)$upper)`), against a simulated `r sprintf("%.2f", simulated_at(20))`. The interval includes the simulated value, and it also includes no effect of oxygen and effects of the opposite sign, with a posterior probability of `r sprintf("%.2f", at(20)$p_negative)` that low oxygen lowers the failure rate at that temperature. These estimates come from a single simulated assay. Among many assays simulated from the same model, the estimates of the interaction are expected to be distributed around the simulated value, and the estimate from any one of them departs from it by chance, by an amount that its credible interval describes.

The model with the main effect alone applies one value to every temperature. Its effect of `r sprintf("%.2f", c_main$median)` is close to the simulated effect at the reference temperature (`r sprintf("%.2f", effect_at_Tr)` at `r sprintf("%.1f", Tr)`°C), around which most fish failed, and is larger than the simulated effect at the temperatures of a stream (e.g., `r sprintf("%.2f", simulated_at(20))` at 20°C). Its credible band has the same width at every temperature, including those at which no fish failed.

## Thermal tolerance at chosen levels of oxygen

The derived quantities are obtained at any level of the co-stressor by giving it in `levels`. `derive_dctmax()` takes the same argument.

```{r}
#| label: tolerance-code
#| eval: false
fish   <- data.frame(mass = 15)
tol_10 <- derive_tdt_ttl(fit_interaction, newdata = fish, t_ref = 60, levels = c(DO = 10),
                         Ta = seq(25.5, 29.5, by = 0.1))
tol_5  <- derive_tdt_ttl(fit_interaction, newdata = fish, t_ref = 60, levels = c(DO = 5),
                         Ta = seq(25.5, 29.5, by = 0.1))
dct_5  <- derive_dctmax(fit_interaction, rate = c(0.1, 0.3), T_start = 14, newdata = fish,
                        levels = c(DO = 5))
```

```{r}
#| label: tbl-tolerance
#| echo: false
#| tbl-cap: "Static critical thermal maximum for an exposure of 60 min (sCTmax) and dynamic critical thermal maximum under ramps from 14°C at two rates (dCTmax) of a 15 g fish in normoxia and in hypoxia, from each of the two models (posterior medians and 95% credible intervals, in °C)."
half <- function(x) x[x$p == 0.5, ]
tolerance_rows <- do.call(rbind, lapply(models, \(m) {
  x <- results$models[[m]]
  do.call(rbind, lapply(c("10", "5"), \(level) data.frame(
    Model = m,
    `Oxygen (mg/L)` = level,
    `sCTmax, 60 min` = estimate(half(x$sctmax[[level]])),
    `dCTmax, 0.1°C per min` = estimate(x$dctmax[[level]][1, ]),
    `dCTmax, 0.3°C per min` = estimate(x$dctmax[[level]][2, ]),
    check.names = FALSE
  )))
}))
knitr::kable(tolerance_rows, align = "lrrrr", row.names = FALSE)
```

```{r}
#| label: tolerance-text
#| include: false
sct <- function(model, level) half(model$sctmax[[level]])
```

Under the model with the interaction, the sCTmax for 60 min was `r sprintf("%.1f", sct(inter, "10")$median)`°C in normoxia and `r sprintf("%.1f", sct(inter, "5")$median)`°C at 5 mg/L (@tbl-tolerance). The two models differed by at most `r sprintf("%.1f", max(abs(c(sct(main, "10")$median - sct(inter, "10")$median, sct(main, "5")$median - sct(inter, "5")$median))))`°C in the sCTmax and by less in the dCTmax, because these quantities refer to the temperatures at which the fish failed in the assay.

The thermal death time curves at the two levels of oxygen are drawn together by joining the tables of the two results, with a column for the level of oxygen, which `plot_tdt()` then treats as it treats a covariate (@fig-tdt).

```{r}
#| label: tdt-code
#| eval: false
both <- list(tdt_curves = rbind(
  data.frame(oxygen = 10, tol_10$tdt_curves),
  data.frame(oxygen = 5, tol_5$tdt_curves)
))
```

```{r}
#| label: tdt-load
#| include: false
both <- list(tdt_curves = rbind(
  data.frame(oxygen = 10, inter$tdt_curves[["10"]]),
  data.frame(oxygen = 5, inter$tdt_curves[["5"]])
))
```

```{r}
#| label: fig-tdt
#| fig-cap: "Thermal death time curves for a fish of 15 g in normoxia (10 mg/L of dissolved oxygen) and in hypoxia (5 mg/L), from the model with the interaction of oxygen with temperature. Lines are posterior medians and bands are 95% credible intervals."
plot_tdt(both, time_unit = "min")
```

The thermal tolerance landscape at each level of oxygen is drawn from the corresponding result (@fig-ttl). At a given temperature, each contour of the probability of survival is reached after a shorter exposure in hypoxia than in normoxia.

```{r}
#| label: ttl-code
#| eval: false
plot_ttl(tol_10, time_unit = "min")
plot_ttl(tol_5, time_unit = "min")
```

```{r}
#| label: fig-ttl
#| echo: false
#| fig-cap: "Thermal tolerance landscape (posterior median of the probability of survival at each constant temperature and exposure time) for a fish of 15 g in normoxia (10 mg/L of dissolved oxygen) and in hypoxia (5 mg/L), from the model with the interaction of oxygen with temperature."
#| fig-subcap:
#|   - "10 mg/L"
#|   - "5 mg/L"
#| layout-ncol: 2
#| fig-width: 4.2
#| fig-height: 3.8
plot_ttl(list(ttl_surface = inter$ttl_surface[["10"]]), time_unit = "min")
plot_ttl(list(ttl_surface = inter$ttl_surface[["5"]]), time_unit = "min")
```

Any level of oxygen may be given, including levels that were not in the assay (e.g., `levels = c(DO = 7.5)`). The assay had two levels, so what the model gives between them or beyond them rests on the assumption that the effect of oxygen on the logarithm of the failure rate is proportional to the deficit, which two levels cannot test.

## Survival along a record of temperature and oxygen

For a prediction along a record, the co-stressor is given in one of two ways. When the record has a column with the name the co-stressor was declared with (here `do`), `predict_survival()` follows it through time as it follows temperature. Otherwise the co-stressor is held at the level given in `levels`.

The dataset `stream_record` holds temperature only, so we add to it a series of dissolved oxygen that falls as the stream warms: 10 mg/L at 15°C or below, and 0.5 mg/L less for each degree above 15°C.

```{r}
#| label: record
record    <- stream_record
record$do <- pmin(10, 10 - 0.5 * (record$temp - 15))
range(record$do)
```

```{r}
#| label: predict-code
#| eval: false
# oxygen following the record
pred_following <- predict_survival(fit_interaction, record, newdata = fish)
# oxygen held at 10 mg/L
pred_normoxia  <- predict_survival(fit_interaction, stream_record, newdata = fish,
                                   levels = c(DO = 10))
```

The record and the survival along it under each of the two models are shown in @fig-record. The survival predicted by a model is drawn with `plot_survival()`, as in the previous vignettes (e.g., `plot_survival(pred_following, temperature = TRUE)`).

```{r}
#| label: fig-record
#| echo: false
#| fig-cap: "Temperature and dissolved oxygen along the stream record, and the probability of survival of a 15 g fish along it, predicted by the model with the main effect of oxygen alone and by the model with its interaction with temperature (lines are posterior medians and bands are 95% credible intervals). The dashed line is the survival under the model used to simulate the data. The exposure is counted from the start of the record."
#| fig-height: 6
# Survival under the model used to simulate the data at each reading of the record,
# with oxygen following the record.
simulated_survival <- function(record) {
  n        <- nrow(record)
  midpoint <- (record$temp[-1] + record$temp[-n]) / 2
  deficit  <- 10 - (record$do[-1] + record$do[-n]) / 2
  oxygen   <- deficit * (truth$b_do + truth$d_do * (midpoint - truth$T_ref))
  accrued  <- exp(truth$log_a + oxygen + truth$b * (midpoint - truth$T_ref)) *
    diff((record$time / truth$scale)^truth$shape)
  c(1, exp(-cumsum(accrued)))
}
# The survival predicted by each model, with oxygen following the record.
predicted <- do.call(rbind, lapply(models, \(m) {
  s <- results$models[[m]]$prediction$following$survival
  stopifnot(identical(s$time, record$time))
  data.frame(model = m, s)
}))
predicted$model <- factor(predicted$model, levels = models)
panels <- c("Temperature (°C)", "Dissolved oxygen (mg/L)", "Survival")
along <- rbind(
  data.frame(panel = panels[1], day = record$time / 1440, value = record$temp),
  data.frame(panel = panels[2], day = record$time / 1440, value = record$do)
)
along$panel <- factor(along$panel, levels = panels)
survival <- data.frame(panel = factor(panels[3], levels = panels),
                       day = predicted$time / 1440, model = predicted$model,
                       value = predicted$median, lower = predicted$lower,
                       upper = predicted$upper)
simulated <- data.frame(panel = factor(panels[3], levels = panels),
                        day = record$time / 1440, value = simulated_survival(record))
# the panel of survival runs from 0 to 1
limits <- data.frame(panel = factor(panels[3], levels = panels), day = 0, value = c(0, 1))

ggplot(mapping = aes(x = day, y = value)) +
  geom_blank(data = limits) +
  geom_line(data = along, colour = "grey20") +
  geom_ribbon(data = survival, aes(ymin = lower, ymax = upper, fill = model), alpha = 0.2) +
  geom_line(data = survival, aes(colour = model)) +
  geom_line(data = simulated, linetype = "dashed", colour = "grey20") +
  scale_colour_manual(values = c("#E69F00", "#225EA8"), name = NULL) +
  scale_fill_manual(values = c("#E69F00", "#225EA8"), name = NULL) +
  facet_wrap(vars(panel), ncol = 1, scales = "free_y", strip.position = "left") +
  scale_x_continuous(breaks = seq(0, 30, by = 5)) +
  labs(x = "Time (days)", y = NULL) +
  theme_light() +
  theme(strip.placement = "outside", strip.background = element_blank(),
        strip.text = element_text(colour = "black"), legend.position = "top")
```

As in the previous vignettes, we also predict the mortality over the five days of the heat wave and over its warmest day, with the exposure counted from the start of each period.

```{r}
#| label: tbl-survival
#| echo: false
#| tbl-cap: "Probability (%) that a 15 g fish fails by the end of the whole stream record, of the five days of its heat wave (days 16 to 20) and of its warmest day (day 18), with dissolved oxygen held at 10 mg/L and with dissolved oxygen falling as the stream warms, predicted by each of the two models (posterior medians and 95% credible intervals, with posterior means in brackets), with the values given by the model used to simulate the data. For the two shorter periods, the exposure is counted from their first reading."
# Mortality under the model used to simulate the data over the days given, with the
# exposure counted from the first reading, and oxygen at 10 mg/L or following the
# record.
simulated_mortality <- function(days, following) {
  part     <- record[record$day %in% days, ]
  n        <- nrow(part)
  midpoint <- (part$temp[-1] + part$temp[-n]) / 2
  deficit  <- if (following) 10 - (part$do[-1] + part$do[-n]) / 2 else 0
  elapsed  <- part$time - part$time[1]
  oxygen   <- deficit * (truth$b_do + truth$d_do * (midpoint - truth$T_ref))
  accrued  <- exp(truth$log_a + oxygen + truth$b * (midpoint - truth$T_ref)) *
    diff((elapsed / truth$scale)^truth$shape)
  1 - exp(-sum(accrued))
}
cell <- function(x) {
  sprintf("%.0f (%.0f–%.0f) [%.0f]", 100 * x$median, 100 * x$lower, 100 * x$upper,
          100 * x$mean)
}
scenarios <- c(normoxia = "10 mg/L", following = "Falling as the stream warms")
rows <- do.call(rbind, lapply(names(scenarios), \(sc) {
  fitted <- do.call(rbind, lapply(models, \(m) {
    x <- results$models[[m]]
    data.frame(Oxygen = scenarios[[sc]], Model = m,
               `30 days` = cell(x$prediction[[sc]]$mortality),
               `Heat wave, 5 days` = cell(x$warm_periods[[sc]]$heat_wave),
               `Warmest day, 24 h` = cell(x$warm_periods[[sc]]$warmest_day),
               check.names = FALSE)
  }))
  follows <- sc == "following"
  rbind(fitted, data.frame(
    Oxygen = scenarios[[sc]], Model = "Model used to simulate the data",
    `30 days` = percent(simulated_mortality(1:30, follows)),
    `Heat wave, 5 days` = percent(simulated_mortality(16:20, follows)),
    `Warmest day, 24 h` = percent(simulated_mortality(18, follows)),
    check.names = FALSE))
}))
knitr::kable(rows, align = "llrrr", row.names = FALSE)
```

```{r}
#| label: survival-text
#| include: false
wave <- function(model, scenario) model$warm_periods[[scenario]]$heat_wave
```

```{r}
#| label: survival-text-2
#| include: false
b_of <- function(model) model$params[model$params$quantity == "b", ]
simulated_wave <- c(normoxia = simulated_mortality(16:20, FALSE),
                    following = simulated_mortality(16:20, TRUE))
hypoxic_failed <- min(ramp$ctmax[ramp$do == 5 & ramp$event == 1])
# The factor by which falling oxygen multiplies the cumulative hazard over the heat
# wave, from the mortality with oxygen at 10 mg/L and with oxygen falling.
hazard_ratio <- function(x) sprintf("%.1f", log(1 - x[["following"]]) / log(1 - x[["normoxia"]]))
fitted_wave  <- function(model) c(normoxia  = wave(model, "normoxia")$median,
                                  following = wave(model, "following")$median)
# what the paragraphs below state about the predictions
stopifnot(
  within(wave(main, "normoxia"), simulated_wave[["normoxia"]]),
  within(wave(main, "following"), simulated_wave[["following"]]),
  within(wave(inter, "normoxia"), simulated_wave[["normoxia"]]),
  within(wave(inter, "following"), simulated_wave[["following"]]),
  wave(main, "following")$median > 2 * wave(main, "normoxia")$median,
  wave(inter, "following")$median - wave(inter, "normoxia")$median < 0.1,
  wave(inter, "normoxia")$median > wave(main, "normoxia")$median,
  b_of(inter)$median < b_of(main)$median,
  max(stream_record$temp) < hypoxic_failed
)
```

Under the model used to simulate the data, the mortality over the heat wave is `r percent(simulated_wave[["normoxia"]])`% with oxygen at 10 mg/L and `r percent(simulated_wave[["following"]])`% with oxygen falling as the stream warms (@tbl-survival). The model with the main effect alone predicted `r percent(wave(main, "normoxia")$median)`% (95% credible interval: `r interval(wave(main, "normoxia"), 0, 100)`) and `r percent(wave(main, "following")$median)`% (`r interval(wave(main, "following"), 0, 100)`), and the model with the interaction `r percent(wave(inter, "normoxia")$median)`% (`r interval(wave(inter, "normoxia"), 0, 100)`) and `r percent(wave(inter, "following")$median)`% (`r interval(wave(inter, "following"), 0, 100)`). The credible intervals of both models include the simulated values, but the two models differ in how much of the mortality they attribute to oxygen: under the first, falling oxygen more than doubles the mortality over the heat wave, and under the second it changes the mortality little. In terms of the cumulative hazard over the heat wave, falling oxygen multiplies it by `r hazard_ratio(simulated_wave)` under the model used to simulate the data, by `r hazard_ratio(fitted_wave(main))` under the model with the main effect alone and by `r hazard_ratio(fitted_wave(inter))` under the model with the interaction (on the posterior medians), so the first attributes more to oxygen than the simulating model does and the second less, while predicting a higher mortality in normoxia.

The difference follows from the effect of oxygen that each model gives at the temperatures of the stream, which never exceeded `r sprintf("%.1f", max(stream_record$temp))`°C (@fig-effect). The model with the main effect alone applies there the effect that it estimated at the temperatures at which the fish failed, which is larger than the simulated effect at the temperatures of the stream. The model with the interaction gives there an effect near zero or negative, with a wide credible band. The model with the interaction also predicted a higher mortality with oxygen at 10 mg/L, because its thermal sensitivity $b$ was lower (`r sprintf("%.2f", b_of(inter)$median)` against `r sprintf("%.2f", b_of(main)$median)` per °C; @tbl-params) and a lower $b$ gives a higher failure rate at temperatures below the reference temperature. The two models give similar failure rates at the reference temperature, around which the fish failed, and at a temperature $\Delta T$ degrees below it the failure rate is $a \exp(-b \, \Delta T)$, which falls less steeply with a lower $b$.

Both predictions extend the effect of oxygen to temperatures below those at which fish failed in hypoxia in the assay (`r sprintf("%.1f", hypoxic_failed)`°C or above), where the assay does not show how the effect of oxygen changes with temperature.

```{r}
#| label: ceiling-text
#| include: false
# The cumulative hazard H that corresponds to a mortality of 1 - exp(-H).
cumulative_hazard <- function(x) {
  data.frame(lower = -log(1 - x$lower), upper = -log(1 - x$upper))
}
month   <- inter$prediction$following$mortality
H_month <- cumulative_hazard(month)
H_wave  <- cumulative_hazard(wave(inter, "following"))
# what the paragraph below states
stopifnot(
  is.finite(H_month$upper),
  month$upper - month$lower < wave(inter, "following")$upper - wave(inter, "following")$lower,
  H_month$upper / H_month$lower > H_wave$upper / H_wave$lower
)
```

Under the model with the interaction and with oxygen falling, the 95% credible interval of the mortality was narrower over the 30 days (`r interval(month, 0, 100)`%) than over the heat wave (`r interval(wave(inter, "following"), 0, 100)`%). The mortality over the 30 days is therefore the more certain of the two predictions, and the reason is the level of the hazard, not the precision with which it was estimated. Mortality is obtained from the cumulative hazard $H$ as $1 - \exp(-H)$ and cannot exceed 100%. Over the 30 days, the 95% credible interval of $H$ was `r interval(H_month, 1)`. Its upper limit was `r sprintf("%.0f", H_month$upper / H_month$lower)` times the lower one, but most fish fail at any value within it, so the interval of the mortality spanned `r sprintf("%.0f", 100 * (month$upper - month$lower))` percentage points. Over the heat wave, the interval of $H$ was `r interval(H_wave, 1)`. Its upper limit was only `r sprintf("%.0f", H_wave$upper / H_wave$lower)` times the lower one, but mortality changes the most with $H$ at such intermediate values, and the interval of the mortality spanned `r sprintf("%.0f", 100 * (wave(inter, "following")$upper - wave(inter, "following")$lower))` percentage points.

These predictions are extrapolations in the two ways described in the previous vignettes: the stream is cooler than the temperatures at which most fish failed in the assay, and the record is much longer than the trials. They also assume that there is no recovery during the cooler hours of each night ([Ørsted et al. 2022](https://doi.org/10.1242/jeb.244514){target="_blank"}). The model used to simulate the data has the form of the fitted models at every temperature and exposure time and makes the same assumption, so the agreement of a prediction with the simulated value says nothing about whether the extrapolations or the assumption hold for real organisms (see the vignettes [*Ramping assays*](ramping-assays.html) and [*Comparing and averaging models*](comparing-and-averaging-models.html)).

## Design of an assay with a co-stressor

The effect of a co-stressor, and how it changes with temperature, are informed only by the temperatures at which individuals fail, and in a ramping assay those cover a few degrees. Within them the two models of this vignette gave similar answers, and outside them they are two different extrapolations: a constant effect, and a straight line that sooner or later changes sign. We therefore recommend: (i) fitting the model with the main effect alone as well as the model with the interaction, and comparing what the two imply at the temperatures for which predictions are needed (@fig-effect), in addition to comparing them by cross-validation (see the vignette [*Comparing and averaging models*](comparing-and-averaging-models.html)); (ii) where the two disagree at those temperatures, reporting the predictions of both and treating the effect of the co-stressor below the temperatures at which individuals failed as an extrapolation; (iii) where predictions at lower temperatures are the aim, including in the assay trials in which individuals fail at such temperatures at each level of the co-stressor (e.g., static trials, or ramps slow enough for failures to occur at lower temperatures); (iv) running three or more levels of the co-stressor where its effect at intermediate levels matters, because with two levels the main effect is the difference between two treatments and its proportionality to the level cannot be assessed; and (v) running several trials of each combination of ramp rate and level of the co-stressor, because the co-stressor is applied to whole trials and its effect is estimated among them.

A model may have more than one co-stressor, each declared with `costressor()` and all given in one list. The derived quantities and the predictions then require a level for each.

```{r}
#| label: two-costressors
#| eval: false
# dissolved copper, in micrograms per litre, with its effect expressed per 10 units
copper <- costressor("Cu", column = "copper", ref = 0, scale = 10, direction = 1)
long <- build_hazard_data(ind, path, observation = "continuous",
                          costressors = list(DO, copper))
fit  <- fit_hazard(long, m_form = "exp", baseline = "weibull",
                   costressors = list(DO, copper))
derive_tdt_ttl(fit, t_ref = 60, levels = c(DO = 7, Cu = 20))
```

The help page of `costressor()` describes the further options: a co-stressor that is harmful in both directions from an optimum (e.g., pH), a treatment with categories in place of a measured level, and an interaction of a different kind, in which each unit of the co-stressor multiplies the thermal sensitivity $b$ by a constant factor (obtained by adding the column of the co-stressor to the formula of `logb` in `dpar_formulas`), whereas the interaction of this vignette adds a constant amount to $b$.

The two kinds of interaction are nearly equivalent when the co-stressor changes $b$ little. The interaction of this vignette has the form of an interaction in a regression and applies to every failure rate function. The proportional kind keeps $b$ positive at every level of the co-stressor, which may matter when the additive kind would make $b$ zero or negative at a level for which predictions are needed, and it is the kind that applies to a treatment with categories. Both models can be fitted and compared by cross-validation (see the vignette [*Comparing and averaging models*](comparing-and-averaging-models.html)).

::: {.callout-note}
## Precomputed results

The models in this vignette were fitted on `r results$stamp$date` with the development version of `rTFT` that became 0.1.0 (`r results$stamp$version`), by the script `vignettes/_precompute.R` in the package's repository. The vignette displays the saved summaries, and all other code runs when the vignette is built.
:::

## References

Ern, R., Norin, T., Gamperl, A. K. and Esbaugh, A. J. (2016). Oxygen dependence of upper thermal limits in fishes. *Journal of Experimental Biology* 219, 3376–3383. [https://doi.org/10.1242/jeb.143495](https://doi.org/10.1242/jeb.143495){target="_blank"}

Ørsted, M., Jørgensen, L. B. and Overgaard, J. (2022). Finding the right thermal limit: a framework to reconcile ecological, physiological and methodological aspects of CTmax in ectotherms. *Journal of Experimental Biology* 225, jeb244514. [https://doi.org/10.1242/jeb.244514](https://doi.org/10.1242/jeb.244514){target="_blank"}