## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE,
  warning = FALSE,
  fig.width = 7,
  fig.height = 4.5,
  out.width = "100%"
)

## ----setup--------------------------------------------------------------------
library(ambre)
set.seed(2024)

## ----libs---------------------------------------------------------------------
library(dplyr)
library(purrr)

## ----chain, results = "hide"--------------------------------------------------
scenario <- create_scenario(
  system.file("input_1culture_2pop.xlsx", package = "ambre")
)
scenario <- inflow_concentration(scenario, pathogenName = "Campylobacter jejuni")
scenario <- scenario |> mutate(volume = map(config, ~ simulate_exposure(config = .x)))
scenario <- initial_dose_calculation(scenario)
scenario <- update_treatment_scheme(scenario, initial_situation = TRUE)
scenario <- scenario |> mutate(log_reduction = map(config, ~ simulate_treatment(.x)))
scenario <- final_dose_calculation(scenario)
scenario <- infection_probability_calculation(scenario)
scenario <- illness_probability_calculation(scenario)
scenario <- dalys_calculation(scenario)
scenario <- get_risk_total(scenario)

## ----unnest-------------------------------------------------------------------
rt <- bind_rows(scenario$risk_total)
nrow(rt)         # 1000 runs x 2 scenario rows

## ----quantiles----------------------------------------------------------------
quantile(rt$dalys_sum, c(0.05, 0.5, 0.95))
mean(rt$dalys_sum > 1e-6)     # share of Monte-Carlo runs above the WHO target

## ----dalys-plot, results = "hide"---------------------------------------------
library(dplyr)
scenario_example <- create_scenario(filepath = system.file("input_1culture_2pop.xlsx", package = "ambre"))
regulation_reduction <- config_ambre$regulation$regulation_value |> filter(Country == "France") |>
                          select(-c(Concentration, Country, RegulationID))
regulation_concentration <- config_ambre$regulation$regulation_value |> filter(Country == "France") |>
                        select(-c(Country, RegulationID, Reduction))

plots <- plot_comparison_qmra_initial_vs_supplementary_processes(
  scenario = scenario_example,
  pathogen = c("Campylobacter jejuni"),
  regulationLog = regulation_reduction,
  regulationConcentration = regulation_concentration
)

## ----dalys-show---------------------------------------------------------------
plots$dalys