Case study: Drosophila suzukii lethal
thermal death time across sexes
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.
The dataset
dsuzukii is the per-individual mortality data from the
Ørsted D. suzukii assays, vendored from bayesTLS
under CC BY 4.0 (primary deposit 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.
str(mort)
#> 'data.frame': 94 obs. of 5 variables:
#> $ temp : num 34 34 34 34 34 34 34 34 34 34 ...
#> $ time : num 1022 106 1168 146 212 ...
#> $ sex : Factor w/ 2 levels "F","M": 1 2 1 1 2 1 2 2 1 2 ...
#> $ n_dead : int 15 1 15 0 0 2 0 0 3 5 ...
#> $ n_total: int 15 15 15 15 15 15 15 15 15 15 ...
table(mort$sex)
#>
#> F M
#> 45 49
There are 94 cells (a temperature × duration × sex grid), 45 for
females and 49 for males, totalling 1407 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.
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
#> [1] 0
The optimiser converged (code 0, positive-definite
Hessian). The per-sex thermal limits, with profile-likelihood
confidence intervals:
confint(fit, c("CTmax:F", "CTmax:M", "z:F", "z:M"), method = "profile")[
, c("parameter", "estimate", "conf.low", "conf.high", "method")]
#> # A tibble: 4 × 5
#> parameter estimate conf.low conf.high method
#> <chr> <dbl> <dbl> <dbl> <chr>
#> 1 CTmax:F 35.2 35.1 35.3 profile
#> 2 CTmax:M 35.3 35.2 35.3 profile
#> 3 z:F 3.01 2.86 3.18 profile
#> 4 z:M 3.18 3.01 3.36 profile
So the maximum-likelihood estimates are:
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.
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.
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:
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` assumes a lethal endpoint; for sublethal data its steeper `z` makes it
#> implausibly low.
#> `T_crit` assumes a lethal endpoint; for sublethal data its steeper `z` makes it
#> implausibly low.
#> T_crit_F T_crit_M
#> 29.21 28.90
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).
sex_contrasts <- confint(
fit, c("dCTmax:F-M", "dz:F-M"), method = "profile",
fallback = TRUE, nboot = 1000L, boot_seed = 20260713L
)
#> Warning: The profile likelihood for "dCTmax:F-M" did not close on the lower and upper
#> sides: "dCTmax:F-M" is weakly identified.
#> ℹ Returning "NA" on the open side rather than a fabricated bound.
#> ℹ Consider bayesTLS or a bootstrap for this parameter.
#> Warning: The profile likelihood for "dz:F-M" did not close on the lower and upper sides:
#> "dz:F-M" is weakly identified.
#> ℹ Returning "NA" on the open side rather than a fabricated bound.
#> ℹ Consider bayesTLS or a bootstrap for this parameter.
#> ! Using a parametric bootstrap for 2 parameters where the profile did not
#> close.
#> ℹ Set `fallback = FALSE` to keep the profile-only behaviour ("NA" on a
#> non-closing side).
sex_contrasts[, c("parameter", "estimate", "conf.low", "conf.high", "method")]
#> # A tibble: 2 × 5
#> parameter estimate conf.low conf.high method
#> <chr> <dbl> <dbl> <dbl> <chr>
#> 1 dCTmax:F-M -0.0244 -0.154 0.0920 bootstrap
#> 2 dz:F-M -0.0538 -0.133 0.0263 bootstrap
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.
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)
)
#> parameter ML_estimate published
#> 1 CTmax:F 35.23 35.20
#> 2 CTmax:M 35.25 35.20
#> 3 z:F 3.01 3.03
#> 4 z:M 3.18 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.
D. suzukii CTmax and z per sex from three estimators (matched
constant-shape configuration, tref = 240 min).
| Female |
CTmax (°C) |
34.80 [34.56, 35.04] |
35.20 [35.08, 35.30] |
35.23 [35.12, 35.32] |
| Female |
z (°C / decade) |
2.99 [2.60, 3.38] |
3.02 [2.88, 3.18] |
3.01 [2.86, 3.18] |
| Male |
CTmax (°C) |
34.98 [34.34, 35.62] |
35.22 [35.10, 35.32] |
35.25 [35.16, 35.34] |
| Male |
z (°C / decade) |
3.40 [2.15, 4.65] |
3.18 [2.98, 3.39] |
3.18 [3.01, 3.36] |
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
sessionInfo()
#> R version 4.6.0 (2026-04-24)
#> Platform: aarch64-apple-darwin23
#> Running under: macOS Tahoe 26.5.2
#>
#> Matrix products: default
#> BLAS: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
#>
#> locale:
#> [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
#>
#> time zone: America/Edmonton
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] freqTLS_0.1.0
#>
#> loaded via a namespace (and not attached):
#> [1] Matrix_1.7-5 gtable_0.3.6 jsonlite_2.0.0 dplyr_1.2.1
#> [5] compiler_4.6.0 tidyselect_1.2.1 jquerylib_0.1.4 scales_1.4.0
#> [9] yaml_2.3.12 fastmap_1.2.0 lattice_0.22-9 ggplot2_4.0.3
#> [13] R6_2.6.1 labeling_0.4.3 generics_0.1.4 knitr_1.51
#> [17] tibble_3.3.1 bslib_0.11.0 pillar_1.11.1 RColorBrewer_1.1-3
#> [21] TMB_1.9.21 rlang_1.2.0 utf8_1.2.6 cachem_1.1.0
#> [25] xfun_0.59 sass_0.4.10 S7_0.2.2 otel_0.2.0
#> [29] viridisLite_0.4.3 cli_3.6.6 withr_3.0.3 magrittr_2.0.5
#> [33] digest_0.6.39 grid_4.6.0 lifecycle_1.0.5 vctrs_0.7.3
#> [37] evaluate_1.0.5 glue_1.8.1 farver_2.1.2 rmarkdown_2.31
#> [41] tools_4.6.0 pkgconfig_2.0.3 htmltools_0.5.9