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

library(dplyr)
library(ggplot2)
library(readxl)

## ----load-packages------------------------------------------------------------
library(EpiExposure)
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"
  )
)

## ----simulate-temperature-----------------------------------------------------
temperature_profile <- simulate_exposures(
  max_lag = 85,
  mode = "profile",
  seed = 123,
  background = list(
    dist = "ar1",
    mean = 25,
    sd = 4,
    phi = 0.90
  )
)

## ----inspect-temperature-profile, echo=FALSE----------------------------------
temperature_profile$profile

## -----------------------------------------------------------------------------

df_tmean <- data.frame(
  time = 0:85,
  value = as.numeric(temperature_profile$profile)
)

 df_tmean%>% 
  ggplot(aes(time,value))+
  geom_line(color = "#C44E52", size = 1.4)+
  theme_bw()+
  labs(x = "Days after planting the crop",
       y = "Exposure")+
  scale_x_continuous(
  limits = c(0, 85),
  breaks = seq(0, 85, by = 10)
)+
  theme(text = element_text(size = 10, face = "bold"))

## ----simulate-rainfall--------------------------------------------------------
rain_profile <- simulate_exposures(
  max_lag = 85,
  mode = "pattern",
  pattern = "block_times",
  block_times = c(
    0,
    40
  ),
  fixed_value = 15,
  background = list(
    dist = "fixed",
    value = 0
  )
)

## ----inspect-rain-profile, echo=FALSE-----------------------------------------
rain_profile$profile

## -----------------------------------------------------------------------------
df_rain <- data.frame(
  time = 0:85,
  value = as.numeric(rain_profile$profile)
)

 df_rain%>% 
  ggplot(aes(time,value))+
  geom_line(color = "#4C72B0", size = 1.4)+
  theme_bw()+
  labs(x = "Days after planting the crop",
       y = "Exposure")+
  scale_x_continuous(
  limits = c(0, 85),
  breaks = seq(0, 85, by = 10)
)+
  theme(text = element_text(size = 10, face = "bold"))

## ----simulate-wetness---------------------------------------------------------
wetness_profile <- simulate_exposures(
  max_lag = 85,
  mode = "profile",
  seed = 456,
  background = list(
    dist = "ar1",
    mean = 10,
    sd = 4,
    phi = 0.90
  )
)

## ----inspect-wetness-profile, echo=FALSE--------------------------------------
wetness_profile$profile

## -----------------------------------------------------------------------------
df_wetness <- data.frame(
  time = 0:85,
  value = as.numeric(wetness_profile$profile)
)


 df_wetness%>% 
  ggplot(aes(time,value))+
  geom_line(color = "#55A868", size = 1.4)+
  theme_bw()+
  labs(x = "Days after planting the crop",
       y = "Exposure")+
  scale_x_continuous(
  limits = c(0, 85),
  breaks = seq(0, 85, by = 10)
)+
  theme(text = element_text(size = 10, face = "bold"))

## ----deterministic-prediction-------------------------------------------------
prediction <- predict_outcomes(
  fit = fit,
  profiles = list(
    tmean = temperature_profile$profile,
    rain = rain_profile$profile,
    wetness = wetness_profile$profile
  ),
  uncertainty = FALSE
)

## -----------------------------------------------------------------------------
prediction

## ----prediction-uncertainty---------------------------------------------------
prediction_samples <- predict_outcomes(
  fit = fit,
  profiles = list(
    tmean = temperature_profile$profile,
    rain = rain_profile$profile,
    wetness = wetness_profile$profile
  ),
  uncertainty = TRUE,
  n_samples = 100,
  output = "samples",
  seed = 123
)

## -----------------------------------------------------------------------------
head(prediction_samples)

## ----inspect-prediction-samples, echo=FALSE-----------------------------------
knitr::kable(
  utils::head(
    prediction_samples,
    10L
  ),
  digits = 4,
  align = "c"
)

## ----plot-prediction-uncertainty, fig.width=7, fig.height=5-------------------
ggplot(
  prediction_samples,
  aes(
    x = prediction * 100
  )
) +
  geom_histogram(
    bins = 25,
    fill = "#4C72B0",
    color = "white"
  ) +
  theme_bw() +
  labs(
    x = "Predicted disease response (%)",
    y = "Number of uncertainty realizations",
    title = "Parameter uncertainty in the predicted response"
  )

## ----prediction-summary-------------------------------------------------------
prediction_summary <- predict_outcomes(
  fit = fit,
  profiles = list(
    tmean = temperature_profile$profile,
    rain = rain_profile$profile,
    wetness = wetness_profile$profile
  ),
  uncertainty = TRUE,
  n_samples = 100,
  output = "summary",
  seed = 123
)

## -----------------------------------------------------------------------------
prediction_summary

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

## -----------------------------------------------------------------------------
periods

## ----generate-scenarios-------------------------------------------------------
scenarios <- simulate_ranges(
  periods = periods,
  vary = list(
    tmean = c(
      20,
      34
    ),
    rain = seq(
      0,
      15,
      by = 1
    ),
    wetness = c(
      2,
      10
    )
  ),
  scenario_type = "grid",
  scenario_var = "rain"
)

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

## ----inspect-scenario-summary, echo=FALSE-------------------------------------
knitr::kable(
  utils::head(
    pred_scenarios,
    10L
  ),
  digits = 4,
  align = "c"
)

## ----plot-scenario-summary, fig.width=9, fig.height=6-------------------------
plot_scenarios(
  data = pred_scenarios,
  x = "rain",
  facet = c(
    "wetness",
    "tmean"
  ),
  x_limits = c(
    0,
    15
  ),
  output = "summary"
)

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

## ----inspect-scenario-samples, echo=FALSE-------------------------------------
knitr::kable(
  utils::head(
    pred_scenarios_s,
    10L
  ),
  digits = 4,
  align = "c"
)

## ----plot-scenario-samples, fig.width=9, fig.height=6-------------------------
plot_scenarios(
  data = pred_scenarios_s,
  x = "rain",
  facet = c(
    "wetness",
    "tmean"
  ),
  x_limits = c(
    0,
    15
  ),
  output = "samples"
)

## ----plot-scenario-grid, fig.width=9, fig.height=6----------------------------
ggplot(
  pred_scenarios,
  aes(
    x = rain,
    y = factor(
      tmean
    ),
    fill = prediction * 100
  )
) +
  geom_tile() +
  facet_wrap(
    vars(
      wetness
    ),
    labeller = label_both
  ) +
  scale_fill_viridis_c(
    name = "Expected\ndisease (%)"
  ) +
  theme_bw() +
  labs(
    x = "Rainfall",
    y = "Mean temperature",
    title = "Predicted response across the evaluated scenario grid"
  )

## ----define-comparison-scenario-a---------------------------------------------
exposure_tmean_A <- simulate_exposures(
  max_lag = 85,
  mode = "profile",
  seed = 101,
  background = list(
    dist = "ar1",
    mean = 25,
    sd = 2,
    phi = 0.90
  ),
  n = 1
)

exposure_rain_A <- simulate_exposures(
  max_lag = 85,
  mode = "profile",
  seed = 102,
  background = list(
    dist = "ar1",
    mean = 5,
    sd = 2,
    phi = 0.50
  ),
  n = 1
)

exposure_wetness_A <- simulate_exposures(
  max_lag = 85,
  mode = "profile",
  seed = 103,
  background = list(
    dist = "ar1",
    mean = 15,
    sd = 3,
    phi = 0.70
  ),
  n = 1
)

scenario_A <- list(
  tmean = exposure_tmean_A,
  rain = exposure_rain_A,
  wetness = exposure_wetness_A
)

## ----define-comparison-scenario-b---------------------------------------------
exposure_tmean_B <- simulate_exposures(
  max_lag = 85,
  mode = "profile",
  seed = 201,
  background = list(
    dist = "ar1",
    mean = 34,
    sd = 2,
    phi = 0.95
  ),
  n = 1
)

exposure_rain_B <- simulate_exposures(
  max_lag = 85,
  mode = "profile",
  seed = 202,
  background = list(
    dist = "ar1",
    mean = 5,
    sd = 1.5,
    phi = 0.40
  ),
  n = 1
)

exposure_wetness_B <- simulate_exposures(
  max_lag = 85,
  mode = "profile",
  seed = 203,
  background = list(
    dist = "ar1",
    mean = 2,
    sd = 0.7,
    phi = 0.60
  ),
  n = 1
)

scenario_B <- list(
  tmean = exposure_tmean_B,
  rain = exposure_rain_B,
  wetness = exposure_wetness_B
)

## ----define-comparison-scenario-c---------------------------------------------
exposure_tmean_C <- simulate_exposures(
  max_lag = 85,
  mode = "profile",
  seed = 301,
  background = list(
    dist = "ar1",
    mean = 15,
    sd = 2,
    phi = 0.90
  ),
  n = 1
)

exposure_rain_C <- simulate_exposures(
  max_lag = 85,
  mode = "profile",
  seed = 302,
  background = list(
    dist = "ar1",
    mean = 0.5,
    sd = 0.5,
    phi = 0.40
  ),
  n = 1
)

exposure_wetness_C <- simulate_exposures(
  max_lag = 85,
  mode = "profile",
  seed = 303,
  background = list(
    dist = "ar1",
    mean = 2,
    sd = 0.8,
    phi = 0.60
  ),
  n = 1
)

scenario_C <- list(
  tmean = exposure_tmean_C,
  rain = exposure_rain_C,
  wetness = exposure_wetness_C
)

## ----define-comparison-scenario-d---------------------------------------------
exposure_tmean_D <- simulate_exposures(
  max_lag = 85,
  mode = "profile",
  seed = 401,
  background = list(
    dist = "ar1",
    mean = 10,
    sd = 2,
    phi = 0.90
  ),
  n = 1
)

exposure_rain_D <- simulate_exposures(
  max_lag = 85,
  mode = "profile",
  seed = 402,
  background = list(
    dist = "ar1",
    mean = 15,
    sd = 4,
    phi = 0.60
  ),
  n = 1
)

exposure_wetness_D <- simulate_exposures(
  max_lag = 85,
  mode = "profile",
  seed = 403,
  background = list(
    dist = "ar1",
    mean = 20,
    sd = 3,
    phi = 0.70
  ),
  n = 1
)

scenario_D <- list(
  tmean = exposure_tmean_D,
  rain = exposure_rain_D,
  wetness = exposure_wetness_D
)

## ----prepare-comparison-profiles----------------------------------------------
profile_to_data <- function(
    exposure,
    variable,
    scenario
) {
  data.frame(
    scenario = scenario,
    var = variable,
    time = 0:85,
    lag = 85:0,
    value = as.numeric(exposure$profile),
    stringsAsFactors = FALSE
  )
}

df_simulations <- dplyr::bind_rows(
  profile_to_data(
    exposure_tmean_A,
    variable = "tmean",
    scenario = "A"
  ),
  profile_to_data(
    exposure_rain_A,
    variable = "rain",
    scenario = "A"
  ),
  profile_to_data(
    exposure_wetness_A,
    variable = "wetness",
    scenario = "A"
  ),
  profile_to_data(
    exposure_tmean_B,
    variable = "tmean",
    scenario = "B"
  ),
  profile_to_data(
    exposure_rain_B,
    variable = "rain",
    scenario = "B"
  ),
  profile_to_data(
    exposure_wetness_B,
    variable = "wetness",
    scenario = "B"
  ),
  profile_to_data(
    exposure_tmean_C,
    variable = "tmean",
    scenario = "C"
  ),
  profile_to_data(
    exposure_rain_C,
    variable = "rain",
    scenario = "C"
  ),
  profile_to_data(
    exposure_wetness_C,
    variable = "wetness",
    scenario = "C"
  ),
  profile_to_data(
    exposure_tmean_D,
    variable = "tmean",
    scenario = "D"
  ),
  profile_to_data(
    exposure_rain_D,
    variable = "rain",
    scenario = "D"
  ),
  profile_to_data(
    exposure_wetness_D,
    variable = "wetness",
    scenario = "D"
  )
)

## ----plot-comparison-profiles, fig.width=10, fig.height=7---------------------
simulation_plot <- ggplot(
  df_simulations,
  aes(
    x = time,
    y = value,
    color = var
  )
) +
  geom_line(linewidth = 1) +
  facet_grid(
    rows = vars(var),
    cols = vars(scenario),
    scales = "free_y"
  ) +
  scale_x_continuous(
    limits = c(
      0,
      85
    ),
    breaks = seq(
      0,
      85,
      by = 15
    )
  ) +
  scale_color_manual(
    values = c(
      tmean = "#C44E52",
      rain = "#4C72B0",
      wetness = "#55A868"
    )
  ) +
  theme_bw() +
  labs(
    x = "Chronological time",
    y = "Exposure",
    color = "Variable"
  ) +
  theme(
    text = element_text(
      size = 10,
      face = "bold"
    ),
    legend.position = "none",
    strip.background = element_rect(
      color = "black",
      fill = "white"
    ),
    strip.text = element_text(
      face = "bold"
    )
  )

simulation_plot

## ----compare-temperature-histories--------------------------------------------
comparison_tmean_AD <- compare_exposures(
  exposure1 = exposure_tmean_A$profile,
  exposure2 = exposure_tmean_D$profile,
  mode = "timewise"
)

comparison_tmean_AD$var <- "tmean"
comparison_tmean_AD$scenario1 <- "A"
comparison_tmean_AD$scenario2 <- "D"

## ----compare-rainfall-histories-----------------------------------------------
comparison_rain_AD <- compare_exposures(
  exposure1 = exposure_rain_A$profile,
  exposure2 = exposure_rain_D$profile,
  mode = "timewise"
)

comparison_rain_AD$var <- "rain"
comparison_rain_AD$scenario1 <- "A"
comparison_rain_AD$scenario2 <- "D"

## ----compare-wetness-histories------------------------------------------------
comparison_wetness_AD <- compare_exposures(
  exposure1 = exposure_wetness_A$profile,
  exposure2 = exposure_wetness_D$profile,
  mode = "timewise"
)

comparison_wetness_AD$var <- "wetness"
comparison_wetness_AD$scenario1 <- "A"
comparison_wetness_AD$scenario2 <- "D"

## ----combine-exposure-comparisons---------------------------------------------
exposure_comparison_AD <- dplyr::bind_rows(
  comparison_tmean_AD,
  comparison_rain_AD,
  comparison_wetness_AD
)

exposure_comparison_AD$var <- factor(
  exposure_comparison_AD$var,
  levels = c(
    "tmean",
    "rain",
    "wetness"
  )
)

## ----inspect-exposure-comparisons, echo=FALSE---------------------------------
knitr::kable(
  utils::head(
    exposure_comparison_AD,
    12L
  ),
  digits = 3,
  align = "c"
)

## ----plot-exposure-comparisons, fig.width=9, fig.height=7---------------------
exposure_comparison_plot <- ggplot(
  exposure_comparison_AD,
  aes(
    x = time,
    y = diff,
    color = var
  )
) +
  geom_hline(
    yintercept = 0,
    color = "grey50",
    linetype = "dashed",
    linewidth = 0.5
  ) +
  geom_line(
    linewidth = 1
  ) +
  facet_wrap(
    vars(var),
    ncol = 1,
    scales = "free_y",
    labeller = as_labeller(
      c(
        tmean = "Mean temperature",
        rain = "Rainfall",
        wetness = "Leaf wetness"
      )
    )
  ) +
  scale_x_continuous(
    limits = c(
      0,
      85
    ),
    breaks = seq(
      0,
      85,
      by = 15
    )
  ) +
  scale_color_manual(
    values = c(
      tmean = "#C44E52",
      rain = "#4C72B0",
      wetness = "#55A868"
    )
  ) +
  theme_bw() +
  labs(
    x = "Days after planting the crop",
    y = "Scenario D - Scenario A",
    color = "Exposure"
  ) +
  theme(
    text = element_text(
      size = 10,
      face = "bold"
    ),
    legend.position = "none",
    strip.background = element_rect(
      fill = "white",
      color = "black"
    ),
    strip.text = element_text(
      face = "bold"
    )
  )

## -----------------------------------------------------------------------------
exposure_comparison_plot

## ----compare-simulated-predictions-summary------------------------------------
comp_pred <- compare_predictions(
  fit = fit,
  profiles = list(
    A = scenario_A,
    B = scenario_B,
    C = scenario_C,
    D = scenario_D
  ),
  uncertainty = TRUE,
  output = "summary",
  n_samples = 100,
  seed = 123
)

comp_pred <- comp_pred |>
  dplyr::mutate(
    comparison = paste(
      scenario1,
      scenario2,
      sep = " vs. "
    )
  )

## ----inspect-simulated-predictions-summary, echo=FALSE------------------------
comparison_summary <- comp_pred |>
  dplyr::select(
    comparison,
    pred1,
    pred1_lower,
    pred1_upper,
    pred2,
    pred2_lower,
    pred2_upper,
    diff,
    diff_lower,
    diff_upper,
    percentage_point_change,
    percentage_point_change_lower,
    percentage_point_change_upper
  )

knitr::kable(
  comparison_summary,
  digits = 4,
  align = "c"
)

## ----compare-simulated-predictions-samples------------------------------------
comp_pred_samples <- compare_predictions(
  fit = fit,
  profiles = list(
    A = scenario_A,
    B = scenario_B,
    C = scenario_C,
    D = scenario_D
  ),
  uncertainty = TRUE,
  output = "samples",
  n_samples = 100,
  seed = 123
)

comp_pred_samples <- comp_pred_samples |>
  dplyr::mutate(
    comparison = paste(
      scenario1,
      scenario2,
      sep = " vs. "
    )
  )

## ----inspect-simulated-predictions-samples, echo=FALSE------------------------
knitr::kable(
  utils::head(
    comp_pred_samples,
    10L
  ),
  digits = 4,
  align = "c"
)

## ----plot-simulated-prediction-comparisons, fig.width=9, fig.height=6---------
comparison_order <- comp_pred |>
  dplyr::arrange(diff) |>
  dplyr::pull(comparison)

comp_pred <- comp_pred |>
  dplyr::mutate(
    comparison = factor(
      comparison,
      levels = comparison_order
    )
  )

comp_pred_samples <- comp_pred_samples |>
  dplyr::mutate(
    comparison = factor(
      comparison,
      levels = comparison_order
    )
  )

sc_simu_plot <- ggplot(
  comp_pred,
  aes(
    x = comparison,
    y = diff
  )
) +
  geom_jitter(
    data = comp_pred_samples,
    aes(
      x = comparison,
      y = diff,
      color = sample
    ),
    alpha = 0.4,
    size = 2,
    width = 0.20,
    height = 0,
    inherit.aes = FALSE
  ) +
  scale_color_viridis_b() +
  geom_pointrange(
    aes(
      ymin = diff_lower,
      ymax = diff_upper
    ),
    color = "black",
    linewidth = 0.6
  ) +
  geom_hline(
    yintercept = 0,
    linetype = 2,
    color = "gray40",
    linewidth = 0.5
  ) +
  theme_bw() +
  labs(
    x = "Scenario comparison",
    y = "Difference in expected response"
  ) +
  theme(
    text = element_text(
      size = 10,
      face = "bold"
    ),
    legend.position = "none",
    axis.text.x = element_text(
      angle = 45,
      hjust = 1
    )
  )

## -----------------------------------------------------------------------------
sc_simu_plot

## ----define-rainfall-only-comparison------------------------------------------
scenario_low_rain <- list(
  tmean = exposure_tmean_A,
  rain = simulate_exposures(
    max_lag = 85,
    mode = "profile",
    seed = 501,
    background = list(
      dist = "ar1",
      mean = 2,
      sd = 1,
      phi = 0.50
    ),
    n = 1
  ),
  wetness = exposure_wetness_A
)

scenario_high_rain <- list(
  tmean = exposure_tmean_A,
  rain = simulate_exposures(
    max_lag = 85,
    mode = "profile",
    seed = 502,
    background = list(
      dist = "ar1",
      mean = 15,
      sd = 3,
      phi = 0.50
    ),
    n = 1
  ),
  wetness = exposure_wetness_A
)

## ----compare-rainfall-only----------------------------------------------------
rainfall_comparison <- compare_predictions(
  fit = fit,
  profiles = list(
    low_rain = scenario_low_rain,
    high_rain = scenario_high_rain
  ),
  uncertainty = TRUE,
  output = "summary",
  n_samples = 100,
  seed = 123
)

## -----------------------------------------------------------------------------
rainfall_comparison