---
title: "Case study: Drosophila suzukii lethal thermal death time across sexes"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Case study: Drosophila suzukii lethal thermal death time across sexes}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6,
  fig.height = 4,
  dpi = 96
)
```

This article fits a lethal thermal-death-time (TDT) model to the vinegar fly
*Drosophila suzukii*, asking the question Ørsted, Hoffmann, Sgrò and colleagues
posed in their 2024 study: **do the two sexes differ in their thermal limits?**
It mirrors the lethal subset of Case Study 4 in the `bayesTLS` supplement, but it
is a `freqTLS`-only treatment: the fit runs live by maximum likelihood (TMB,
**no Stan**), uncertainty is summarised by profile-likelihood **confidence
intervals**, and the sex difference is a **frequentist contrast**, not a
posterior.

**This vignette builds without Stan.** The `freqTLS` side runs live by
maximum likelihood (TMB); the Bayesian (`bayesTLS`) and classical two-stage
columns of the three-way comparison are read from the maintainer-built benchmark
cache (`inst/extdata/bayesTLS_benchmark_cache.rds`), so the full per-sex table
renders without a Stan toolchain.

```{r setup}
library(freqTLS)
```

```{r load-data, include = FALSE}
data(dsuzukii)
# Per-individual records -> survival counts per (temp, time, sex) cell.
.nd  <- aggregate(list(n_dead  = dsuzukii$dead), dsuzukii[c("temp", "time", "sex")], sum)
.nt  <- aggregate(list(n_total = dsuzukii$dead), dsuzukii[c("temp", "time", "sex")], length)
mort <- merge(.nd, .nt)
```

## The dataset

`dsuzukii` is the per-individual mortality data from the Ørsted *D. suzukii*
assays, vendored from [`bayesTLS`](https://github.com/daniel1noble/bayesTLS)
under CC BY 4.0 (primary deposit
[10.5281/zenodo.10602268](https://doi.org/10.5281/zenodo.10602268); cite
`@orsted_suzukii_2024` and `citation("freqTLS")`). We aggregate it to survival
counts per `(temp, time, sex)` cell — `n_dead` of `n_total` died — exactly as the
`bayesTLS` case study does.

```{r data-shape}
str(mort)
table(mort$sex)
```

There are `r nrow(mort)` cells (a temperature × duration × sex grid),
`r sum(mort$sex == "F")` for females and
`r sum(mort$sex == "M")` for males, totalling
`r sum(mort$n_total)` flies. Two columns deserve emphasis:

* **`time` is in minutes.** The reference exposure for this study is **4 hours =
  240 minutes**, so we set `tref = 240`.
* We interpret the threshold as **absolute** (a fixed 50% survival), following
  Ørsted. The model itself is still parameterised by the **relative midpoint**
  (`CTmax` is the 4PL midpoint, the `(low + up)/2` survival level). The relative
  and absolute thresholds coincide only when the asymptotes sit near 0 and 1 —
  which they do for these data (`low` ≈ 0, `up` ≈ 1, confirmed in the fit below),
  so `derive_*` at `surv = 0.5` recovers the absolute LT50 with no extra argument.
  If your own data have asymptotes far from 0/1, the two thresholds diverge and you
  must pass the `surv` you want explicitly.

## Fitting the model: one fit, grouped by sex

A single grouped call fits a beta-binomial 4PL with sex-specific thermal
location (`CTmax`, `z`) and a shared curve shape (`low`, `up`, `k`) — the
constant-shape configuration that matches the `bayesTLS` "sex on the midpoint
only" model. The temperature effect runs through the midpoint, so each sex gets
its own `CTmax` and `z` while the asymptotes and steepness are pooled.

```{r fit}
std <- standardize_data(
  mort, temp = "temp", duration = "time",
  n_total = "n_total", n_dead = "n_dead", duration_unit = "minutes"
)
# fit_4pl(by = "sex") gives a per-sex CTmax and z (the "sex on the midpoint"
# model); $fit is the engine fit the extractors below read directly.
fit <- fit_4pl(std, by = "sex", t_ref = 240, family = "beta_binomial",
               quiet = TRUE)$fit
fit$convergence$code   # 0 = converged
```

The optimiser converged (code `0`, positive-definite Hessian). The per-sex
thermal limits, with **profile-likelihood confidence intervals**:

```{r per-sex}
confint(fit, c("CTmax:F", "CTmax:M", "z:F", "z:M"), method = "profile")[
  , c("parameter", "estimate", "conf.low", "conf.high", "method")]
```

So the maximum-likelihood estimates are:

| Quantity | Female | Male |
| --- | --- | --- |
| `CTmax` at 4 h (°C) | 35.23 [35.13, 35.32] | 35.25 [35.16, 35.34] |
| `z` (°C / decade) | 3.01 [2.86, 3.18] | 3.18 [3.01, 3.36] |

Both sexes have an upper thermal limit near 35.2 °C at a four-hour exposure, and
a thermal sensitivity `z` of about 3 °C per tenfold change in survival time. The
intervals are confidence intervals from the profile likelihood — the range of
parameter values the data do not reject at the 95% level — not credible
intervals and not posterior summaries.

### The Confidence Eye

The `freqTLS` uncertainty display is the **Confidence Eye**: a confidence
lens with a hollow point estimate, carrying no prior and making no probability
statement about the parameter. Here it shows `CTmax` and `z` for both sexes
side by side.

```{r eye, fig.height = 5, fig.alt = "Confidence Eye for Drosophila suzukii CTmax and z, one lens per sex (female and male). The lenses for the two sexes overlap almost completely for both CTmax (near 35.2 C) and z (near 3 C per decade), each drawn as a pale confidence lens with a hollow point estimate."}
plot_confidence_eye(fit, parm = c("CTmax", "z"), method = "profile")
```

The female and male lenses overlap heavily for both quantities — the first
visual hint that the sexes are not strongly separated. Each lens is a closed
confidence interval (the profiles closed on both sides), so the eyes are
drawn closed rather than hollow-and-open.

### Survival curves

The fitted survival curves, with the observed cell proportions overlaid, show
the 4PL decline of survival with exposure time at each temperature, per sex.

```{r survival, fig.height = 4.5, fig.alt = "Fitted Drosophila suzukii survival curves by temperature (34 to 38 C) and sex, with observed cell proportions overlaid. Survival falls from near one to near zero as exposure time increases, and the curves shift to shorter survival times at higher temperatures; the female and male curves are nearly indistinguishable at each temperature."}
plot_survival_curves(fit)
```

### A lower thermal threshold: T_crit

Because this is a **lethal** endpoint, the rate-multiplier critical temperature
`T_crit` is meaningful here (it is a damage-accumulation concept and would be
misleading for sublethal endpoints such as heat-coma). `derive_tcrit()` returns
the temperature at which the thermal-damage rate falls to a chosen low floor; at
the `bayesTLS` default floor of 1% of the lethal dose per hour:

```{r tcrit}
c(
  T_crit_F = round(derive_tcrit(fit, rate = 1, group = "F"), 2),
  T_crit_M = round(derive_tcrit(fit, rate = 1, group = "M"), 2)
)
```

`T_crit` sits near **29 °C** for both sexes (about 6 °C below `CTmax`), the lower
thermal threshold a heat-injury model would treat as the onset of negligible
damage. It is a deterministic transform of the fitted `CTmax` and `z`, not a new
fit.

## Do the sexes differ? A frequentist contrast

The biological question is whether the small female-minus-male gaps in `CTmax`
and `z` are distinguishable from zero. `freqTLS` answers this with a profile
or bootstrap confidence interval **on the difference itself** — a frequentist
contrast, *not* a posterior probability statement. The contrast is requested by
name: `dCTmax:F-M` (the `CTmax` difference, female minus male) and `dz:F-M` (the
sensitivity difference, taken on the log-`z` scale so the `z` ratio is its
exponential).

```{r contrast}
sex_contrasts <- confint(
  fit, c("dCTmax:F-M", "dz:F-M"), method = "profile",
  fallback = TRUE, nboot = 1000L, boot_seed = 20260713L
)
sex_contrasts[, c("parameter", "estimate", "conf.low", "conf.high", "method")]
```

Both intervals **span zero**:

* `dCTmax:F-M` = **-0.024 °C, 95% CI [-0.154, 0.092]** — the sexes' four-hour
  `CTmax` values are indistinguishable.
* `dz:F-M` = **-0.054 [-0.133, 0.026]** on the log-`z` scale (a `z` ratio of
  `exp(-0.054)` ≈ 0.95, a point estimate of about a 5% difference in `z` between
  the sexes), with a confidence interval that includes 0 (ratio 1).

This agrees with the published finding: Ørsted reported a sex difference in `z`
whose interval spans zero, i.e. **no clear sex difference in thermal sensitivity
or limit**. The likelihood path reaches the same conclusion as the Bayesian path
did, by a contrast that makes no posterior claim. The `method` column above
records when the requested profile used the asymmetry-respecting parametric-
bootstrap fallback without a prior.

## Validation: does ML recover the published values?

A useful cross-check is whether the maximum-likelihood fit lands on the same
parameter values Ørsted (2024) reported in their Table 1, which were obtained by
a different route. The published values are `z` ≈ 3.03 (F) / 3.28 (M) and a
four-hour `CTmax` ≈ 35.2 °C for both sexes.

```{r validation}
ml <- confint(fit, c("CTmax:F", "CTmax:M", "z:F", "z:M"), method = "profile")
data.frame(
  parameter   = ml$parameter,
  ML_estimate = round(ml$estimate, 2),
  published   = c(35.2, 35.2, 3.03, 3.28)
)
```

The match is close: `CTmax` ≈ 35.23 / 35.25 °C against a published ≈ 35.2 °C, and
`z` ≈ 3.01 / 3.18 against a published 3.03 / 3.28. The likelihood path reproduces
the published thermal-limit estimates on the same data — a validation success
that the `freqTLS` engine targets the same model and the same fitted curve as
the analyses that established these numbers.

## The three-way comparison, per sex

The cache holds the maintainer-built `bayesTLS` (posterior) and classical
two-stage summaries; the `freqTLS` column is computed live as this page
renders. All three use the matched configuration (beta-binomial, constant shape,
`tref = 240` minutes). The `CTmax` compared is the relative-midpoint parameter;
for these near-0/near-1 lethal asymptotes it coincides with the absolute LT50 the
study reports above, so the benchmark and the headline per-sex table agree.

```{r three-way-suzukii, echo = FALSE}
cache_path <- system.file("extdata", "bayesTLS_benchmark_cache.rds", package = "freqTLS")
has_cache  <- nzchar(cache_path) && file.exists(cache_path)

fmt <- function(est, lo, hi) {
  ifelse(is.na(est), "--", sprintf("%.2f [%.2f, %.2f]", est, lo, hi))
}

sex_levels <- fit$group_levels
sex_label  <- c(F = "Female", M = "Male")

# Live freqTLS per-sex CTmax and z as profile-likelihood confidence
# intervals (with explicit open-profile status; no Stan).
pro <- suppressMessages(confint(
  fit,
  c(paste0("CTmax:", sex_levels), paste0("z:", sex_levels)),
  method = "profile"
))
pr_of <- function(parm, g) {
  r <- pro[pro$parameter == paste0(parm, ":", g), , drop = FALSE]
  if (nrow(r) == 0L) return(stats::setNames(rep(NA_real_, 3), c("est", "lo", "hi")))
  stats::setNames(c(r$estimate[1L], r$conf.low[1L], r$conf.high[1L]),
                  c("est", "lo", "hi"))
}

# The cache keys each sex as "dsuzukii:<sex>".
cache_has_suzukii <- FALSE
if (has_cache) {
  cache <- readRDS(cache_path)
  cache_has_suzukii <- !is.null(cache$bayesian) &&
    "dataset" %in% names(cache$bayesian) &&
    all(paste0("dsuzukii:", sex_levels) %in% cache$bayesian$dataset)
}

if (cache_has_suzukii) {
  pick <- function(df, parm, g, cols) {
    r <- df[df$dataset == paste0("dsuzukii:", g) & df$parameter == parm, , drop = FALSE]
    if (nrow(r) == 0L) return(stats::setNames(rep(NA_real_, 3), c("est", "lo", "hi")))
    stats::setNames(as.numeric(r[1L, cols]), c("est", "lo", "hi"))
  }
  one_row <- function(g, parm, quantity) {
    ts <- pick(cache$two_stage, parm, g, c("estimate", "lower", "upper"))
    by <- pick(cache$bayesian,  parm, g, c("median",   "lower", "upper"))
    pr <- pr_of(parm, g)
    data.frame(
      Sex      = sex_label[[g]],
      Quantity = quantity,
      `Two-stage (delta CI)`    = fmt(ts["est"], ts["lo"], ts["hi"]),
      `bayesTLS (95% CrI)`      = fmt(by["est"], by["lo"], by["hi"]),
      `freqTLS (profile CI)` = fmt(pr["est"], pr["lo"], pr["hi"]),
      check.names = FALSE
    )
  }
  rows <- do.call(rbind, lapply(sex_levels, function(g) rbind(
    one_row(g, "CTmax", "CTmax (°C)"),
    one_row(g, "z",     "z (°C / decade)")
  )))
  knitr::kable(rows, row.names = FALSE,
               caption = "D. suzukii CTmax and z per sex from three estimators (matched constant-shape configuration, tref = 240 min).")
} else {
  prof_tbl <- do.call(rbind, lapply(sex_levels, function(g) {
    ct <- pr_of("CTmax", g); z <- pr_of("z", g)
    data.frame(
      Sex = sex_label[[g]],
      `CTmax (°C), profile CI`      = fmt(ct["est"], ct["lo"], ct["hi"]),
      `z (°C / decade), profile CI` = fmt(z["est"], z["lo"], z["hi"]),
      check.names = FALSE
    )
  }))
  knitr::kable(prof_tbl, row.names = FALSE,
               caption = "Live freqTLS per-sex estimates (no Stan). The bayesTLS and two-stage columns appear once the maintainer-built cache carries the D. suzukii sexes.")
}
```

Under the matched configuration the `freqTLS` profile-likelihood confidence
intervals and the `bayesTLS` credible intervals summarise the *same* fitted curve
per sex — one prior-free and by optimisation, the other with a prior and MCMC —
and the two model-based fits agree tightly, both reproducing the published thermal
limits (`CTmax` ≈ 35.2 °C; `z` ≈ 3.0 for females and 3.2 for males). The classical
two-stage estimate sits a few tenths of a degree lower on `CTmax` (≈ 34.8–35.0 °C)
with wider, less precise intervals — the expected behaviour of the cruder
per-temperature method, and a reminder that the headline equivalence here is
between the two model-based fits (see `vignette("comparing-to-bayesTLS")` and
`vignette("model-math")`).

## Boundary: lethal only

`freqTLS` is purpose-built for the 4PL survival/proportion model class. The
Ørsted study measured two further *D. suzukii* endpoints that are **non-goals**
here:

* **Heat-coma (knockdown)** is a right-censored time-to-event endpoint
  (`brms::bf(... | cens(cens) ~ ...)`), a model class outside the 4PL count
  engine. Its rate-multiplier `T_crit` is also endpoint-conditional (knockdown
  `z` < lethal `z`), so it must not be reported as a lethal threshold. For the
  heat-coma analysis, **see `bayesTLS`**.
* **Productivity (fertility)** is modelled with a hurdle-Gamma (a zero-inflation
  process plus a positive clutch-size process), two response processes, neither
  a 4PL survival curve, and the magnitude component has no LT50 or `z` at all.
  For productivity, **see `bayesTLS`**.

This article covers the **lethal-by-sex** subset only, which is the part of the
*D. suzukii* analysis that fits the `freqTLS` model class.

**Where to next:** `vignette("heat-injury")` walks through the `T_crit` /
heat-injury workflow that this absolute-threshold fit feeds.

## Session info

```{r session}
sessionInfo()
```