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

## ----message=FALSE, warning=FALSE---------------------------------------------

library(EpiExposure)
library(dplyr)
library(ggplot2)
library(patchwork)

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

data("epi_data")

head(epi_data, 10L )

## ----echo=TRUE, eval=FALSE----------------------------------------------------
# library(future)
# 
# future::plan(
#   future::multisession,
#   workers = 4
# )

## ----echo=TRUE, eval=FALSE----------------------------------------------------
# 
# future::plan(future::sequential )
# 

## ----echo=TRUE, eval=FALSE----------------------------------------------------
# best_models <- find_bestfit(
#   data = epi_data,
#   response = "y",
#   group = "epi_id",
#   time = "time",
#   vars = c("tmean","rain","wetness"),
#   max_lag = 85,
#   df_var_grid = c(2,3,4),
#   df_lag_grid = c(2,3,4),
#   model_engine = "glmmTMB",
#   family = "beta",
#   keep_fits = TRUE)

## ----include=FALSE------------------------------------------------------------
#best_models = readRDS("best_models.rds")
best_models = readxl::read_xlsx("data/ensemble/best_models.xlsx")


## ----echo=TRUE, eval=FALSE----------------------------------------------------
# plot_performance(
#   object = best_models,
#   model_id = c(5,12,19,26,28,33,35,40,42,47,49,54,56,61,63),
#   metrics = c("CCC"),
#   x = "observed",
#   y = "predicted",
#   scale_factor = 1,
#   ncol = 5,
#   x_lab = "Observed ",
#   y_lab = "Predicted",
#   size = 1.4
# ) +
#   ggplot2::theme(
#     text = ggplot2::element_text(
#       size = 10,
#       face = "bold"
#     )
#   )

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

## -----------------------------------------------------------------------------
best_models %>%
  arrange(desc(CCC)) %>%
  select(
    model_id,
    vars,
    CCC,
    RMSE,
    MAE
  ) %>%
  head(10)

## -----------------------------------------------------------------------------
best_models %>%
  mutate(
    splines = paste(
      df_var,
      df_lag,
      sep = " × "
    )
  ) %>%
  ggplot(
    aes(
      Cb,
      CCC,
      color = splines
    )
  ) +
  geom_point(size = 3) +
  theme_bw()

## -----------------------------------------------------------------------------
top_model <- best_models |>
  dplyr::arrange(
    dplyr::desc(
      CCC
    )
  ) |>
  dplyr::slice(
    1L
  )

## -----------------------------------------------------------------------------
candidate_models <- c(
  63,
  23,
  44
)

## ----echo=TRUE, eval=FALSE----------------------------------------------------
# ens_unweighted <- ensemble_bestfit(
#   bestfit = best_models,
# 
#   ensemble_scope = "model",
# 
#   method = "unweighted",
# 
#   model_id = candidate_models)

## ----eval = FALSE-------------------------------------------------------------
# ens_unweighted$ensemble_summary

## ----include=FALSE------------------------------------------------------------
ens_unweighted <- readxl::read_xlsx("data/ensemble/ens_unweighted.xlsx")

## ----echo=FALSE---------------------------------------------------------------
knitr::kable(
  ens_unweighted,
  digits = 4,
  align = "c"
)

## ----echo=TRUE, eval=FALSE----------------------------------------------------
# ens_weighted <- ensemble_bestfit(
#   bestfit = best_models,
# 
#   ensemble_scope = "model",
# 
#   method = "weighted",
# 
#   model_id = candidate_models,
# 
#   weight_metric = "CCC",
# 
#   weight_transform = "rank_inverse"
# )

## ----eval = FALSE-------------------------------------------------------------
# ens_weighted$ensemble_summary

## ----include=FALSE------------------------------------------------------------
ens_weighted <- readxl::read_xlsx("data/ensemble/ens_weighted.xlsx")

## ----echo=FALSE---------------------------------------------------------------
knitr::kable(
  ens_weighted,
  digits = 4,
  align = "c"
)

## ----echo=TRUE, eval=FALSE----------------------------------------------------
# ens_stacked <- ensemble_bestfit(
#   bestfit = best_models,
# 
#   ensemble_scope = "model",
# 
#   method = "stacked",
# 
#   model_id = candidate_models,
# 
#   stacking_model = "ridge",
# 
#   stack_objective = "regularized"
# )

## ----eval = FALSE-------------------------------------------------------------
# ens_stacked$ensemble_summary

## ----include=FALSE------------------------------------------------------------
ens_stacked <- readxl::read_xlsx("data/ensemble/ens_stacked.xlsx")

## ----echo=FALSE---------------------------------------------------------------
knitr::kable(
  ens_stacked,
  digits = 4,
  align = "c"
)

## ----echo=TRUE, eval=FALSE----------------------------------------------------
# ensemble_metrics <- dplyr::bind_rows(
#   Unweighted = ens_unweighted$ensemble_summary,
#   Weighted = ens_weighted$ensemble_summary,
#   Stacked = ens_stacked$ensemble_summary,
#   .id = "method"
# ) |>
#   dplyr::select(
#     method,
#     CCC,
#     RMSE,
#     MAE
#   )

## ----include=FALSE------------------------------------------------------------
ensemble_metrics <- readxl::read_xlsx("data/ensemble/ensemble_metrics.xlsx")

## -----------------------------------------------------------------------------
ensemble_metrics

## ----echo=TRUE, eval=FALSE----------------------------------------------------
# 
# ens_unweighted$ensemble_predictions$method = "unweighted"
# ens_weighted$ensemble_predictions$method = "weighted"
# ens_stacked$ensemble_predictions$method = "stacked"
# 
# ensemble_all = rbind(ens_unweighted$ensemble_predictions,ens_weighted$ensemble_predictions, ens_stacked$ensemble_predictions)

## ----include=FALSE------------------------------------------------------------
ensemble_all <- readxl::read_xlsx("data/ensemble/ensemble_all.xlsx")

## -----------------------------------------------------------------------------
knitr::kable(utils::head(ensemble_all, 10L),
              digits = 4, align = "c" )

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

ensemble_metrics$method = c("unweighted","weighted","stacked")

ordem_methods <- c(
  "unweighted",
  "weighted",
  "stacked"
)

ensemble_plot <- ensemble_all %>%
  filter(!is.na(method)) %>%
  mutate(
    method = as.character(method)
  ) %>%
  filter(method %in% ordem_methods) %>%
  mutate(
    method = factor(
      method,
      levels = ordem_methods
    )
  )



CCC_plot <- ensemble_metrics %>%
  mutate(
    method = factor(
      method,
      levels = ordem_methods
    )
  )



pos_CCC <- ensemble_plot %>%
  group_by(method) %>%
  summarise(
    x = min(observed, na.rm = TRUE) +
      0.05 * diff(range(observed, na.rm = TRUE)),
   
    y = max(predicted_ensemble, na.rm = TRUE) -
      0.05 * diff(range(predicted_ensemble, na.rm = TRUE)),
   
    .groups = "drop"
  ) %>%
  left_join(
    CCC_plot %>% select(method, CCC),
    by = "method"
  )



 ensemble_gg = ggplot(ensemble_plot,
            aes(observed, predicted_ensemble)) +
  geom_point()+
  geom_smooth(method = "lm",se = FALSE, color = "orange") +
  geom_abline(
    intercept = 0,
    slope = 1,
    color = "red",
    linetype = "dashed",
    linewidth = 1) +
  geom_text(
    data = pos_CCC,
    aes(
      x = 0.02,
      y = 0.97,
      label = sprintf("CCC = %.2f", CCC)
    ),
    inherit.aes = FALSE,
    hjust = 0,
    vjust = 1,
    size = 4) +
  facet_wrap(
    ~method,
    ncol = 1) +
  theme_bw() +
    scale_x_continuous(
    limits = c(0, 1),
    breaks = seq(0, 1, 0.20),
    expand = expansion(mult = 0.03)
  ) +
 
  scale_y_continuous(
    limits = c(0, 1),
    breaks = seq(0, 1, 0.20),
    expand = expansion(mult = 0.03)
  ) +
  labs(
    x = "Observed",
    y = "Predicted") +
  theme(text = element_text(size = 10,face = "bold"),legend.position = "none")
 
 ensemble_gg

## -----------------------------------------------------------------------------
best_models_plot <- best_models %>%
   mutate(
    splines = paste(df_lag, df_var, sep = " x ")
  )


best_plot <- best_models_plot %>%
  ggplot(aes(
    Cb, CCC,
    color = as.factor(splines),
    shape = as.factor(vars)
  )) +
 
  geom_point(size = 4) +
 
  scale_color_viridis_d() +
 
  labs(
    x = "Cb",
    y = "CCC",
    color = "Splines",
    shape = "Predictors"
  ) +
 
  guides(
    color = guide_legend(
      position = "top",
      nrow = 1,
      byrow = TRUE,
      title.position = "left",
      theme = theme(
        legend.key.height = unit(0.25, "cm"),
        legend.key.width = unit(0.45, "cm"),
        legend.spacing.y = unit(0, "cm"),
        legend.margin = margin(
          t = -3,
          r = 0,
          b = -3,
          l = 0
        )
      )
    ),
   
    shape = guide_legend(
      position = "right"
    )
  ) +
 
  theme_bw() +
 
  theme(
    text = element_text(
      size = 10,
      face = "bold"
    )
  )

## -----------------------------------------------------------------------------
final_plot <- ensemble_gg + best_plot +
  plot_layout(ncol = 2) +
  plot_annotation(
    tag_levels = "a",
    tag_prefix = "(",
    tag_suffix = ")"
  ) &
  theme(
    plot.tag = element_text(
      size = 12,
      face = "bold"
    ),
    plot.tag.position = c(0.02, 0.94)
  )

final_plot


## -----------------------------------------------------------------------------
epi_data1 = epi_data |> 
  filter(epi_id == 1)

## ----fit-lag-ensembles, eval = FALSE------------------------------------------
# lag_ensemble <- ensemble_bestfit(
#   bestfit = best_models,
#   data = epi_data1,
#   group = "epi_id",
#   var = c(
#     "tmean",
#     "rain",
#     "wetness"
#   ),
#   ensemble_scope = "lag",
#   method = c(
#     "unweighted",
#     "weighted"
#   ),
#   top_n = 3,
#   weight_metric = NULL,
#   weight_transform = "softmax",
#   lag_group_cols = c(
#     "epi_id",
#     "var",
#     "lag"
#   ),
#   lg_strategy = "requested_available",
#   rl_weights = TRUE,
#   compute_ecilag_args = list(
#     uncertainty = FALSE
#   ),
#   verbose = TRUE
# )
# 

## ----eval = FALSE-------------------------------------------------------------
# lag_ensemble$selected_models

## ----echo=FALSE---------------------------------------------------------------
selected_models = readxl::read_xlsx("data/ensemble/selected_models.xlsx")

## ----inspect-lag-selected-models, echo=FALSE----------------------------------
knitr::kable(
  selected_models,
  digits = 4,
  align = "c"
)

## ----eval = FALSE-------------------------------------------------------------
# lag_ensemble$model_weights

## ----echo=FALSE---------------------------------------------------------------
model_weights = readxl::read_xlsx("data/ensemble/model_weights.xlsx")

## ----inspect-lag-model-weights, echo=FALSE------------------------------------
knitr::kable(
  model_weights,
  digits = 4,
  align = "c"
)

## ----eval = FALSE-------------------------------------------------------------
# lag_ensemble$ensemble_by_lag

## ----echo=FALSE---------------------------------------------------------------
ensemble_by_lag = readxl::read_xlsx("data/ensemble/ensemble_by_lag.xlsx")

## ----inspect-lag-ensemble, echo=FALSE-----------------------------------------
knitr::kable(
  utils::head(
    ensemble_by_lag,
    12L
  ),
  digits = 4,
  align = "c"
)

## ----eval = FALSE-------------------------------------------------------------
# ensemble_by_lag = lag_ensemble$ensemble_by_lag

## ----plot-lag-ensemble-signed, fig.width=10, fig.height=7---------------------
lag_ensemble_signed_plot <- ggplot(
  ensemble_by_lag,
  aes(
    x = lag,
    y = ECI_weighted_ens,
    color = method
  )
) +
  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_reverse(
    breaks = seq(
      0,
      85,
      by = 15
    )
  ) +
  scale_color_manual(
    values = c(
      unweighted = "#4C72B0",
      weighted = "#C44E52"
    ),
    labels = c(
      unweighted = "Unweighted",
      weighted = "Weighted"
    )
  ) +
  theme_bw() +
  labs(
    x = "Retrospective lag",
    y = "Ensemble weighted ECI",
    color = "Ensemble method"
  ) +
  theme(
    text = element_text(
      size = 10,
      face = "bold"
    ),
    strip.background = element_rect(
      fill = "white",
      color = "black"
    ),
    strip.text = element_text(
      face = "bold"
    ),
    legend.position = "top"
  )

## -----------------------------------------------------------------------------
lag_ensemble_signed_plot

## ----plot-lag-ensemble-percent, fig.width=10, fig.height=7--------------------
lag_ensemble_percent_plot <- ggplot(
  ensemble_by_lag,
  aes(
    x = lag,
    y = ECI_percent_ens,
    color = method
  )
) +
  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_reverse(
    breaks = seq(
      0,
      85,
      by = 15
    )
  ) +
  scale_color_manual(
    values = c(
      unweighted = "#4C72B0",
      weighted = "#C44E52"
    ),
    labels = c(
      unweighted = "Unweighted",
      weighted = "Weighted"
    )
  ) +
  theme_bw() +
  labs(
    x = "Retrospective lag",
    y = "Absolute lag contribution (%)",
    color = "Ensemble method"
  ) +
  theme(
    text = element_text(
      size = 10,
      face = "bold"
    ),
    strip.background = element_rect(
      fill = "white",
      color = "black"
    ),
    strip.text = element_text(
      face = "bold"
    ),
    legend.position = "top"
  )

## -----------------------------------------------------------------------------
lag_ensemble_percent_plot

## ----strict-lag-ensemble, eval=FALSE------------------------------------------
# lag_ensemble_strict <- ensemble_bestfit(
#   bestfit = best_models,
#   data = epi_data1,
#   group = "epi_id",
#   var = c(
#     "tmean",
#     "rain",
#     "wetness"
#   ),
#   ensemble_scope = "lag",
#   method = c(
#     "unweighted",
#     "weighted"
#   ),
#   top_n = 3,
#   weight_transform = "softmax",
#   lag_group_cols = c(
#     "epi_id",
#     "var",
#     "lag"
#   ),
#   lg_strategy = "strict",
#   rl_weights = TRUE,
#   compute_ecilag_args = list(
#     uncertainty = FALSE
#   )
# )