## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.align = "center",
  message = FALSE,
  warning = FALSE
)

## ----load-packages------------------------------------------------------------
library(EpiExposure)
library(dplyr)
library(ggplot2)

## ----load-data----------------------------------------------------------------
data(
  "epi_data")

## ----define-exposures---------------------------------------------------------
cb <- define_exposures(
  data = epi_data,
  vars = c(
    "tmean",
    "rain",
    "wetness"
  ),
  max_lag = 85,
  df_var = 3,
  df_lag = 2
)

## ----build-design-------------------------------------------------------------
dat <- build_design(
  data = epi_data,
  cb_templates = cb
)

## ----prepare-response---------------------------------------------------------
dat <- prepare_response(
  data = dat,
  response = "y",
  family = "beta"
)

## ----fit-model----------------------------------------------------------------
fit <- fit_epidlnm(
  data = dat,
  model_engine = "glmmTMB",
  family = "beta",
  random_effect = NULL,
  epiexposure_spec = attr(
    cb,
    "spec"
  )
)

## ----define-periods-----------------------------------------------------------
periods <- define_periods(
  max_lag = 85,
  cuts = c(
    20,
    40,
    60
  ),
  prefix = "Period"
)

periods

## ----temperature-gradient-----------------------------------------------------
temp_range <- simulate_ranges(
  periods = periods,
  vary = list(
    rain = 5,
    wetness = 10,
    tmean = seq(
      18,
      34,
      by = 0.5
    )
  ),
  scenario_type = "grid",
  scenario_var = "tmean"
)

## ----predict-temperature-gradient---------------------------------------------
temp_predictions <- simulate_scenarios(
  fit = fit,
  scenarios = temp_range,
  data = epi_data,
  uncertainty = TRUE,
  output = "summary",
  n_samples = 100,
  seed = 123
)

## ----prepare-temperature-response---------------------------------------------
response_curve <- temp_predictions |>
  dplyr::transmute(
    value = tmean,
    prediction = prediction * 100,
    lower = prediction_lower * 100,
    upper = prediction_upper * 100
  ) |>
  dplyr::arrange(
    value
  )

## ----inspect-temperature-response, echo=FALSE---------------------------------
knitr::kable(
  utils::head(
    response_curve,
    10L
  ),
  digits = 3,
  caption = "First 10 predicted responses along the temperature gradient.",
  align = "c"
)

## ----plot-temperature-response, fig.width=8, fig.height=5---------------------
ggplot(
  response_curve,
  aes(
    x = value,
    y = prediction
  )
) +
  geom_ribbon(
    aes(
      ymin = lower,
      ymax = upper
    ),
    fill = "#C44E52",
    alpha = 0.18
  ) +
  geom_line(
    linewidth = 1,
    color = "#C44E52"
  ) +
  theme_bw() +
  labs(
    x = "Mean temperature (°C)",
    y = "Expected disease response (%)",
    title = "Predicted response along the temperature gradient"
  )

## ----compute-temperature-sensitivity------------------------------------------
temp_sensitivity <- epi_sensitivity(
  data = response_curve,
  x = "value",
  y = "prediction",
  scenario_var = NULL,
  elasticity = TRUE,
  critical = TRUE,
  k = 3
)

## ----inspect-temperature-sensitivity, echo=FALSE------------------------------
knitr::kable(
  utils::head(
    temp_sensitivity,
    10L
  ),
  digits = 4,
  caption = "First 10 temperature sensitivity results.",
  align = "c"
)

## ----plot-temperature-sensitivity, fig.width=8, fig.height=5------------------
ggplot(
  temp_sensitivity,
  aes(
    x = value,
    y = sensitivity
  )
) +
  geom_hline(
    yintercept = 0,
    color = "grey60",
    linetype = "dashed"
  ) +
  geom_line(
    linewidth = 1,
    color = "#C44E52"
  ) +
  theme_bw() +
  labs(
    x = "Mean temperature (°C)",
    y = "Local sensitivity",
    title = "Sensitivity along the temperature gradient"
  )

## ----plot-temperature-elasticity, fig.width=8, fig.height=5-------------------
ggplot(
  temp_sensitivity,
  aes(
    x = value,
    y = elasticity
  )
) +
  geom_hline(
    yintercept = 0,
    color = "grey60",
    linetype = "dashed"
  ) +
  geom_line(
    linewidth = 1,
    color = "#4C72B0"
  ) +
  theme_bw() +
  labs(
    x = "Mean temperature (°C)",
    y = "Elasticity",
    title = "Elasticity along the temperature gradient"
  )

## ----inspect-critical-temperature-points--------------------------------------
critical_temperature <- temp_sensitivity |>
  dplyr::filter(
    critical
  )

critical_temperature

## ----plot-critical-temperature-points, fig.width=8, fig.height=5--------------
ggplot(
  temp_sensitivity,
  aes(
    x = value,
    y = sensitivity
  )
) +
  geom_hline(
    yintercept = 0,
    color = "grey60",
    linetype = "dashed"
  ) +
  geom_line(
    linewidth = 0.9,
    color = "#C44E52"
  ) +
  geom_point(
    data = critical_temperature,
    color = "#B22222",
    size = 2.8
  ) +
  theme_bw() +
  labs(
    x = "Mean temperature (°C)",
    y = "Local sensitivity",
    title = "Model-derived critical regions along the temperature gradient"
  )

## ----rainfall-gradient--------------------------------------------------------
rain_range <- simulate_ranges(
  periods = periods,
  vary = list(
    tmean = 25,
    wetness = 10,
    rain = seq(
      0,
      15,
      by = 0.5
    )
  ),
  scenario_type = "grid",
  scenario_var = "rain"
)

## ----predict-rainfall-gradient------------------------------------------------
rain_predictions <- simulate_scenarios(
  fit = fit,
  scenarios = rain_range,
  data = epi_data,
  uncertainty = TRUE,
  output = "samples",
  n_samples = 100,
  seed = 123
)

## ----compute-rainfall-sensitivity---------------------------------------------
sens_rain <- epi_sensitivity(
  data = rain_predictions,
  x = "rain",
  y = "prediction",
  scenario_var = "sample",
  elasticity = TRUE,
  critical = TRUE,
  method = "finite",
  k = 3
)

## ----wetness-gradient---------------------------------------------------------
wetness_range <- simulate_ranges(
  periods = periods,
  vary = list(
    tmean = 25,
    rain = 5,
    wetness = seq(
      0,
      24,
      by = 0.5
    )
  ),
  scenario_type = "grid",
  scenario_var = "wetness"
)

## ----predict-wetness-gradient-------------------------------------------------
wetness_predictions <- simulate_scenarios(
  fit = fit,
  scenarios = wetness_range,
  data = epi_data,
  uncertainty = TRUE,
  output = "samples",
  n_samples = 100,
  seed = 456
)

## ----compute-wetness-sensitivity----------------------------------------------
sens_wetness <- epi_sensitivity(
  data = wetness_predictions,
  x = "wetness",
  y = "prediction",
  scenario_var = "sample",
  elasticity = TRUE,
  critical = TRUE,
  method = "finite",
  k = 3
)

## ----combine-sensitivity-results----------------------------------------------
sens_rain_plot <- sens_rain |>
  dplyr::mutate(
    variable = "Rainfall",
    exposure_value = rain
  )

sens_wetness_plot <- sens_wetness |>
  dplyr::mutate(
    variable = "Leaf wetness",
    exposure_value = wetness
  )

combined_sensitivity <- dplyr::bind_rows(
  sens_rain_plot,
  sens_wetness_plot
)

## ----plot-combined-sensitivity-samples, fig.width=9, fig.height=5-------------
ggplot(
  combined_sensitivity,
  aes(
    x = exposure_value,
    y = sensitivity,
    group = sample
  )
) +
  geom_hline(
    yintercept = 0,
    color = "grey60",
    linetype = "dashed"
  ) +
  geom_line(
    alpha = 0.10,
    linewidth = 0.35,
    color = "#4C72B0"
  ) +
  facet_wrap(
    vars(
      variable
    ),
    scales = "free_x"
  ) +
  theme_bw() +
  labs(
    x = "Exposure value",
    y = "Local sensitivity",
    title = "Sensitivity across environmental gradients"
  )

## ----summarize-sensitivity-samples--------------------------------------------
sensitivity_summary <- combined_sensitivity |>
  dplyr::group_by(
    variable,
    exposure_value
  ) |>
  dplyr::summarise(
    sensitivity_median = stats::median(
      sensitivity,
      na.rm = TRUE
    ),
    sensitivity_lower = stats::quantile(
      sensitivity,
      probs = 0.025,
      na.rm = TRUE
    ),
    sensitivity_upper = stats::quantile(
      sensitivity,
      probs = 0.975,
      na.rm = TRUE
    ),
    .groups = "drop"
  )

## ----inspect-sensitivity-summary, echo=FALSE----------------------------------
knitr::kable(
  utils::head(
    sensitivity_summary,
    10L
  ),
  digits = 4,
  caption = "First 10 summarized sensitivity estimates.",
  align = "c"
)

## ----plot-sensitivity-summary, fig.width=9, fig.height=5----------------------
ggplot(
  sensitivity_summary,
  aes(
    x = exposure_value,
    y = sensitivity_median
  )
) +
  geom_hline(
    yintercept = 0,
    color = "grey60",
    linetype = "dashed"
  ) +
  geom_ribbon(
    aes(
      ymin = sensitivity_lower,
      ymax = sensitivity_upper
    ),
    fill = "#4C72B0",
    alpha = 0.18
  ) +
  geom_line(
    linewidth = 0.9,
    color = "#4C72B0"
  ) +
  facet_wrap(
    vars(
      variable
    ),
    scales = "free_x"
  ) +
  theme_bw() +
  labs(
    x = "Exposure value",
    y = "Local sensitivity",
    title = "Summarized sensitivity across environmental gradients"
  )