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

## ----estimate-lag-effects-----------------------------------------------------
epi_l <- summarise_effects(
  fit = fit,
  data = epi_data,
  scale = "lag",
  effect_measure = "percent",
  uncertainty = TRUE,
  n_samples = 10,
  output = "summary",
  at = list(
    tmean = seq(
      15,
      35,
      length.out = 100
    ),
    rain = seq(
      0,
      15,
      length.out = 100
    ),
    wetness = seq(
      0,
      24,
      length.out = 100
    )
  )
)

## ----inspect-lag-effects, echo=FALSE------------------------------------------
knitr::kable(
  utils::head(
    epi_l,
    10L
  ),
  digits = 3,
  align = "c"
)

## ----plot-lag-effects, fig.width=9, fig.height=7------------------------------
plot_effects(
  data = epi_l,
  scale = "lag",
  metric = "effect",
  vars = c(
    "tmean",
    "rain",
    "wetness"
  )
)

## ----plot-delta-effects, fig.width=9, fig.height=7----------------------------
plot_effects(
  data = epi_l,
  scale = "lag",
  metric = "delta",
  delta_multiplier = 100,
  vars = c(
    "tmean",
    "rain",
    "wetness"
  )
)

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

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

## ----estimate-period-effects--------------------------------------------------
epi_p_summary <- summarise_effects(
  fit = fit,
  data = epi_data,
  scale = "period",
  lag_periods = periods,
  effect_measure = "percent",
  uncertainty = TRUE,
  n_samples = 20,
  output = "summary"
)

## ----inspect-period-effects, echo=FALSE---------------------------------------
knitr::kable(
  utils::head(
    epi_p_summary,
    10L
  ),
  digits = 3,
  align = "c"
)

## ----plot-period-effects, fig.width=9, fig.height=7---------------------------
plot_effects(
  data = epi_p_summary,
  scale = "period",
  vars = c(
    "tmean",
    "rain",
    "wetness"
  ),
  period_order = c(
    "Period4",
    "Period3",
    "Period2",
    "Period1"
  ),
  vars_order = c(
    "tmean",
    "rain",
    "wetness"
  ),
  period_response = "eta",
  period_ylab = "Effect"
)

## ----period-effect-samples----------------------------------------------------
epi_p_samples <- summarise_effects(
  fit = fit,
  data = epi_data,
  scale = "period",
  lag_periods = periods,
  effect_measure = "percent",
  uncertainty = TRUE,
  n_samples = 20,
  output = "samples"
)

## ----plot-period-effects2, fig.width=9, fig.height=7--------------------------
plot_effects(
  data = epi_p_samples,
  scale = "period",
  vars = c(
    "tmean",
    "rain",
    "wetness"
  ),
  period_order = c(
    "Period4",
    "Period3",
    "Period2",
    "Period1"
  ),
  vars_order = c(
    "tmean",
    "rain",
    "wetness"
  ),
  period_response = "eta",
  period_ylab = "Effect"
)

## ----inspect-period-samples, echo=FALSE---------------------------------------
knitr::kable(
  utils::head(
    epi_p_samples,
    10L
  ),
  digits = 3,
  align = "c"
)

## ----compare-periods, eval=FALSE----------------------------------------------
# epi_p_samples <- summarise_effects(
#   fit = fit,
#   data = epi_data,
#   scale = "period",
#   lag_periods = periods,
#   effect_measure = "percent",
#   uncertainty = TRUE,
#   n_samples = 20,
#   output = "samples"
# )
# 
# period_comparison <- compare_periods(
#   period_df = epi_p_samples,
#   period_ref = "Period1",
#   periods = periods,
#   output = "summary"
# )

## ----reduce-overall-temperature-----------------------------------------------
reduced_tmean <- reduce_effects(
  fit = fit,
  data = epi_data,
  vars = "tmean",
  group = "epi_id",
  time = "time",
  type = "overall",
  scale = "percent",
  uncertainty = TRUE,
  output = "summary",
  n_samples = 10,
  seed = 123,
  ref = list(
    method = "median",
    value = NULL
  ),
  at = seq(
    15,
    35,
    length.out = 20
  )
)

## ----inspect-reduced-temperature, echo=FALSE----------------------------------
knitr::kable(
  utils::head(
    reduced_tmean,
    10L
  ),
  digits = 3,
  align = "c"
)

## ----plot-reduced-temperature, fig.width=7, fig.height=5----------------------
ggplot(
  reduced_tmean,
  aes(
    x = x,
    y = effect
  )
) +
  geom_ribbon(
    aes(
      ymin = low_eff,
      ymax = high_eff
    ),
    fill = "#C44E52",
    alpha = 0.20
  ) +
  geom_hline(
    yintercept = 0,
    color = "grey50",
    linetype = "dashed"
  ) +
  geom_line(
    color = "#C44E52",
    linewidth = 1
  ) +
  theme_bw() +
  labs(
    x = "Mean temperature (°C)",
    y = "Change in modeled odds (%)",
    title = "Overall cumulative temperature association"
  ) +
  theme(
    text = element_text(
      size = 10,
      face = "bold"
    )
  )

## ----reduce-rainfall-selected-lag---------------------------------------------
reduced_rain_lag10 <- reduce_effects(
  fit = fit,
  data = epi_data,
  vars = "rain",
  group = "epi_id",
  time = "time",
  type = "lag",
  value = 10,
  scale = "percent",
  uncertainty = TRUE,
  output = "summary",
  n_samples = 10,
  seed = 123,
  ref = list(
    method = "median",
    value = NULL
  ),
  at = seq(
    0,
    15,
    length.out = 20
  )
)

## ----plot-reduced-rainfall-lag, fig.width=7, fig.height=5---------------------
ggplot(
  reduced_rain_lag10,
  aes(
    x = x,
    y = effect
  )
) +
  geom_ribbon(
    aes(
      ymin = low_eff,
      ymax = high_eff
    ),
    fill = "#4C72B0",
    alpha = 0.20
  ) +
  geom_hline(
    yintercept = 0,
    color = "grey50",
    linetype = "dashed"
  ) +
  geom_line(
    color = "#4C72B0",
    linewidth = 1
  ) +
  theme_bw() +
  labs(
    x = "Rainfall",
    y = "Change in modeled odds (%)",
    title = "Rainfall association at lag 10"
  ) +
  theme(
    text = element_text(
      size = 10,
      face = "bold"
    )
  )

## ----reduce-wetness-selected-value--------------------------------------------
reduced_wetness_10 <- reduce_effects(
  fit = fit,
  data = epi_data,
  vars = "wetness",
  group = "epi_id",
  time = "time",
  type = "var",
  value = 10,
  scale = "percent",
  uncertainty = TRUE,
  output = "summary",
  n_samples = 10,
  seed = 123,
  ref = list(
    method = "median",
    value = NULL
  )
)

## ----plot-reduced-wetness-value, fig.width=7, fig.height=5--------------------
ggplot(
  reduced_wetness_10,
  aes(
    x = x,
    y = effect
  )
) +
  geom_ribbon(
    aes(
      ymin = low_eff,
      ymax = high_eff
    ),
    fill = "#55A868",
    alpha = 0.20
  ) +
  geom_hline(
    yintercept = 0,
    color = "grey50",
    linetype = "dashed"
  ) +
  geom_line(
    color = "#55A868",
    linewidth = 1
  ) +
  scale_x_reverse(
    breaks = seq(
      0,
      85,
      by = 15
    )
  ) +
  theme_bw() +
  labs(
    x = "Retrospective lag",
    y = "Change in modeled odds (%)",
    title = "Lag-response association for leaf wetness = 10"
  ) +
  theme(
    text = element_text(
      size = 10,
      face = "bold"
    )
  )

## ----reduce-effects-samples, eval=FALSE---------------------------------------
# reduced_tmean_samples <- reduce_effects(
#   fit = fit,
#   data = epi_data,
#   vars = "tmean",
#   group = "epi_id",
#   time = "time",
#   type = "overall",
#   scale = "percent",
#   uncertainty = TRUE,
#   output = "samples",
#   n_samples = 40,
#   seed = 123,
#   ref = list(
#     method = "median",
#     value = NULL
#   ),
#   at = seq(
#     15,
#     35,
#     length.out = 100
#   )
# )

## ----compute-eci--------------------------------------------------------------
eci <- compute_eci(
  data = epi_data,
  group = "epi_id",
  fit = fit,
  vars = c(
    "tmean",
    "rain",
    "wetness"
  ),
  uncertainty = TRUE,
  output = "summary",
  n_samples = 10
)

## ----inspect-eci, echo=FALSE--------------------------------------------------
knitr::kable(
  utils::head(
    eci,
    10L
  ),
  digits = 3,
  align = "c"
)

## ----compare-eci-tmean, fig.width=7, fig.height=5-----------------------------
eci[
  eci$var == "tmean",
  ,
  drop = FALSE
] |>
  ggplot(
    aes(
      x = ECI_raw,
      y = ECI_weighted
    )
  ) +
  geom_point(
    alpha = 0.65
  ) +
  theme_bw() +
  labs(
    x = "Raw ECI",
    y = "Weighted ECI",
    title = "tmean")

## ----compare-eci-wetness, fig.width=7, fig.height=5---------------------------
eci[
  eci$var == "wetness",
  ,
  drop = FALSE
] |>
  ggplot(
    aes(
      x = ECI_raw,
      y = ECI_weighted
    )
  ) +
  geom_point(
    alpha = 0.65
  ) +
  theme_bw() +
  labs(
    x = "Raw ECI",
    y = "Weighted ECI",
    title = "wetness")

## ----compare-eci-rain, fig.width=7, fig.height=5------------------------------
eci[
  eci$var == "rain",
  ,
  drop = FALSE
] |>
  ggplot(
    aes(
      x = ECI_raw,
      y = ECI_weighted
    )
  ) +
  geom_point(
    alpha = 0.65
  ) +
  theme_bw() +
  labs(
    x = "Raw ECI",
    y = "Weighted ECI",
    title = "rain")

## ----compute-eci-lag----------------------------------------------------------
eci_lag <- compute_ecilag(
  data = epi_data,
  group = "epi_id",
  group_level = 1,
  fit = fit,
  var = c(
    "tmean",
    "rain",
    "wetness"
  ),
  uncertainty = TRUE,
  output = "summary",
  n_samples = 10
)

## ----inspect-eci-lag, echo=FALSE----------------------------------------------
knitr::kable(
  utils::head(
    eci_lag$by_lag,
    10L
  ),
  digits = 3,
  align = "c"
)

## ----plot-eci, fig.width=9, fig.height=9--------------------------------------
plot_eci(
  eci_overall = eci,
  eci_lag = eci_lag,
  facet_scales = "free_y",
  overall_nrow = 3,
  lag_ncol = 1,
  sample_smooth = TRUE,
  linewidth = 0.7,
  alpha = 0.20,
  panel_labels = c(
    "(a)",
    "(b)",
    "(c)"
  ),
  overall_col = "ECI_weighted",
  lag_col = "ECI_percent"
)