## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.align = "center",
  message = FALSE,
  warning = FALSE
)

## ----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"
  )
)

## ----define-loss-periods------------------------------------------------------
periods <- define_periods(
  max_lag = 85,
  cuts = c(
    20,
    40,
    60
  ),
  prefix = "Period"
)

## ----generate-loss-scenarios--------------------------------------------------
range_losses <- 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-loss-scenarios---------------------------------------------------
scenarios_losses <- simulate_scenarios(
  fit = fit,
  scenarios = range_losses,
  data = epi_data,
  uncertainty = TRUE,
  output = "summary",
  n_samples = 5,
  seed = 123
)

## ----convert-disease-scale----------------------------------------------------
scenarios_losses$prediction <- (
  scenarios_losses$prediction * 100
)

scenarios_losses$prediction_lower <- (
  scenarios_losses$prediction_lower * 100
)

scenarios_losses$prediction_upper <- (
  scenarios_losses$prediction_upper * 100
)

## ----simulate-losses----------------------------------------------------------
losses <- simulate_losses(
  data = scenarios_losses,
  lower = "prediction_lower",
  upper = "prediction_upper",
  y = "prediction",
  slope = 49.3,
  parameter_mode = "values",
  intercept = 9691.0,
  attainable_yield = seq(
    4000,
    10000,
    by = 50
  ),
  price = 300
)

## ----inspect-loss-columns-----------------------------------------------------
names(
  losses
)

## ----inspect-loss-results, echo=FALSE-----------------------------------------
knitr::kable(
  utils::head(
    losses,
    10L
  ),
  digits = 3,
  align = "c"
)

## ----plot-loss-heatmap, fig.width=10, fig.height=6----------------------------
losses_plot <- plot_losses(
  data = losses,
  x = "rain",
  y = "att_yield",
  facet = c(
    "tmean",
    "wetness"
  ),
  type = "heatmap",
  x_breaks = seq(
    0,
    15,
    by = 3
  ),
  y_breaks = seq(
    4000,
    10000,
    by = 2000
  ),
  x_label = "Daily precipitation (mm)",
  nrow = 2,
  facet_ncol = 4
)

losses_plot

## ----simulate-loss-sensitivity------------------------------------------------
losses_sensitivity <- simulate_losses(
  data = scenarios_losses,
  lower = "prediction_lower",
  upper = "prediction_upper",
  y = "prediction",
  slope = c(
    49.3,
    80
  ),
  parameter_mode = "values",
  intercept = c(
    9691.0,
    12000
  ),
  attainable_yield = seq(
    4000,
    10000,
    by = 50
  ),
  price = 300
)

## ----plot-loss-regression, fig.width=10, fig.height=6-------------------------
plot_losses(
  data = losses_sensitivity,
  x = "rain",
  y = "att_yield",
  facet = c(
    "tmean",
    "wetness"
  ),
  type = "regression",
  reg_y = "ap_yield",
  att_yield = 10000,
  x_breaks = seq(
    0,
    15,
    by = 3
  ),
  x_label = "Daily precipitation (mm)"
)

