## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6,
  fig.height = 3.6,
  dpi = 96
)

## ----setup--------------------------------------------------------------------
library(freqTLS)
library(ggplot2)

## ----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))

## ----tcrit--------------------------------------------------------------------
tcrit <- derive_tcrit(fit, rate = 1)   # 1% of a lethal dose per hour
round(tcrit, 2)

## ----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))

## ----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)

## ----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)

## ----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)

## ----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)

## ----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)

## ----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")

## ----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)

## ----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")

## ----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)
)