What-if: calibrating exposure

library(ambre)
set.seed(2024)
library(dplyr)
library(purrr)

Two numbers decide how much of the water each person actually meets:

Together with the pathogen concentration, these drive the dose, and the dose drives the risk. This vignette is about the exposure side of the model: where those two numbers come from by default, how to read them, and how to override them to ask “what if?” — including the two most common management measures, personal protective equipment and reduced-contact irrigation.

If you are new to the pipeline, read vignette("a-get-started", package = "ambre") first. The barrier-side levers (treatment plants, on-field practices, their log-reductions) are a different story, told in vignette("b-initial-vs-new-scenario", package = "ambre").

How exposure is set by default

You never set volume, frequency or concentration by hand in a normal run. When you call inflow_concentration() (which run_qmra_intial_situation() and run_qmra_supplementary_process() call for you), it runs three updates in order:

inflow_concentration()
  ├─ update_pathogen()    # pathogen concentration
  ├─ update_frequency()   # number of events / year, per path
  └─ update_volume()      # litres per event, per path
scenario <- create_scenario(
  system.file("input_1culture_2pop.xlsx", package = "ambre")
)
scenario_conc <- inflow_concentration(
  scenario = scenario,
  pathogenName = "Campylobacter jejuni"
)

update_pathogen() look up inflow database. update_frequency() and update_volume() look up the exposure path of each scenario row and copy the path-specific numbers out of the shipped database config_ambre$path into that row’s embedded config$exposure table. So the model does not use one generic volume for every path — each crop × population × path combination gets its own exposure profile.

You can see the result. The config$exposure table has three rows — number_of_repeatings (the Monte-Carlo count), number_of_exposures (events per year) and volume_perEvent:

# Row 1 of this two-row scenario: irrigation-staff droplet ingestion
scenario_conc$config[[1]]$exposure
#> # A tibble: 3 × 10
#>   name                 type    value   min   max mode  mean  sd    meanlog sdlog
#>   <chr>                <chr>   <dbl> <dbl> <dbl> <lgl> <lgl> <lgl> <lgl>   <lgl>
#> 1 number_of_repeatings value    1000 NA    NA    NA    NA    NA    NA      NA   
#> 2 number_of_exposures  value      60 NA    NA    NA    NA    NA    NA      NA   
#> 3 volume_perEvent      triang…    NA  1e-3  1e-3 NA    NA    NA    NA      NA

Two quirks are worth knowing up front, and we return to them at the end:

Inspect the base numbers

Before changing anything, look at what the model is assuming.

Start from the PathogenName, you can see the default value of concentration for the simulated pathogen :

# Row 1 of this two-row scenario: Campylobacter jejuni concentration
dplyr::filter(scenario_conc$config[[1]]$inflow, PathogenName == "Campylobacter jejuni")
#> # A tibble: 1 × 13
#>   PathogenID PathogenName   PathogenGroup simulate type  value   min   max mode 
#>        <dbl> <chr>          <chr>            <dbl> <chr> <lgl> <dbl> <dbl> <lgl>
#> 1          1 Campylobacter… Bacteria             1 unif… NA      100  5000 NA   
#> # ℹ 4 more variables: mean <lgl>, sd <lgl>, meanlog <lgl>, sdlog <lgl>

Then start from a path description and get its PathID with query_exp_path():

pid <- query_exp_path(
  pathName = "Ingestion of water droplets during maintenance of the irrigation system"
)
pid
#> [1] 1

Then read the per-path volume and frequency straight from the database. Both helpers accept a vector of IDs, so you can inspect a whole scenario at once:

scenario$PathID          # the two paths in this input file
#> [1] 1 4

query_volume(scenario$PathID)     # litres per event: min / max
#> # A tibble: 2 × 2
#>     min   max
#>   <dbl> <dbl>
#> 1 0.001 0.001
#> 2 0.001 0.001
query_frequency(scenario$PathID)  # events per year: min / max
#> # A tibble: 2 × 2
#>     min   max
#>   <dbl> <dbl>
#> 1    40    60
#> 2    48    48

These are exactly the numbers update_volume() and update_frequency() inject. The volumes here are tiny — 0.001 litre is one millilitre of incidentally swallowed droplets — the frequency 4060 becomes 60 because only the maximum is used and the concentration is between 100 and 500 log/L.

Override for a what-if

To explore sensitivity you replace those numbers with your own. Two helpers do this, and both expect one value per scenario row.

update_volume_with_desired_value() takes a volume argument that is a data frame (or list) with a min, a max and a type column. Use min == max for a fixed volume, or min < max for a distribution spread according to the type specified:

scenario_ppe <- update_volume_with_desired_value(
  scenario = scenario,
  volume = data.frame(
    min = c(0.0005, 0.0005),   # 0.5 mL per event, one entry per row
    max = c(0.0005, 0.0005),
    type = c("triangle", "triangle")
  )
)

scenario_ppe$config[[1]]$exposure
#> # A tibble: 3 × 10
#>   name                 type    value   min   max mode  mean  sd    meanlog sdlog
#>   <chr>                <chr>   <dbl> <dbl> <dbl> <lgl> <lgl> <lgl> <lgl>   <lgl>
#> 1 number_of_repeatings value    1000 NA    NA    NA    NA    NA    NA      NA   
#> 2 number_of_exposures  value     365 NA    NA    NA    NA    NA    NA      NA   
#> 3 volume_perEvent      triang…    NA  5e-4  5e-4 NA    NA    NA    NA      NA

update_frequency_with_desired_value() takes a plain numeric vector of integer, one number of events per year per row:

scenario_night <- update_frequency_with_desired_value(
  scenario = scenario,
  frequency = c(30L, 30L)   # cap both paths at 30 events / year
)

scenario_night$config[[1]]$exposure
#> # A tibble: 3 × 10
#>   name                 type    value   min   max mode  mean  sd    meanlog sdlog
#>   <chr>                <chr>   <dbl> <dbl> <dbl> <lgl> <lgl> <lgl> <lgl>   <lgl>
#> 1 number_of_repeatings value    1000  NA      NA NA    NA    NA    NA      NA   
#> 2 number_of_exposures  value      30  NA      NA NA    NA    NA    NA      NA   
#> 3 volume_perEvent      triang…    NA   0.5     3 NA    NA    NA    NA      NA

The volume_perEvent row now reads 0.0005, and number_of_exposures reads 30: your what-if values have replaced the ones the database would have injected.

update_concentration() takes a data.frame with PathogenName, min, max and type of ditribution law values.

concentration_custom <- data.frame(PathogenName = c("Campylobacter jejuni"),
                                   min = c(1),
                                   max = c(2),
                                   type = c("uniform"))

scenario_pathogen <- update_pathogen(scenario = scenario, pathoName = concentration_custom$PathogenName)
scenario_low_pathogen <- update_concentration(scenario = scenario_pathogen ,
                     concentration = concentration_custom)

dplyr::filter(scenario_low_pathogen$config[[1]]$inflow, PathogenName == "Campylobacter jejuni")
#> # A tibble: 1 × 13
#>   PathogenID PathogenName PathogenGroup simulate value mode  mean  sd    meanlog
#>        <dbl> <chr>        <chr>            <dbl> <lgl> <lgl> <lgl> <lgl> <lgl>  
#> 1          1 Campylobact… Bacteria             1 NA    NA    NA    NA    NA     
#> # ℹ 4 more variables: sdlog <lgl>, min <dbl>, max <dbl>, type <chr>

Modelling management measures

Read those two overrides as interventions. This is the exposure-side complement to the barrier log-reductions in vignette("b-treatment-vs-multibarrier", package = "ambre"): instead of removing pathogens from the water, you change how much of the water each person meets.

Both are assumptions you choose and should defend — they are not credited from a barrier database. (The barriere_path / barriere_specific tables in config_ambre that would encode route-specific exposure reductions are shipped data, not yet wired into the engine.)

A mini sensitivity study

Let us quantify how much the volume assumption matters. One catch first: run_qmra_initial_situation() and run_qmra_supplementary_process() re-derives the per-path volumes from the database on every call (via inflow_concentration()), so it would overwrite any override you set. To inject a what-if value you should use run_qmra_custom function:

Run it twice on the same file and pathogen — once with a pessimistic bare-hand volume, once with a PPE volume ten times smaller:

sc <- create_scenario(system.file("input_1culture_2pop.xlsx", package = "ambre"))

regulation_reduction <- config_ambre$regulation$regulation_value |> 
  dplyr::filter(Country == "France") |>
  dplyr::select(-c(Concentration, Country, RegulationID))
regulation_concentration <- config_ambre$regulation$regulation_value |> 
  dplyr::filter(Country == "France") |>
  dplyr::select(-c(Country, RegulationID, Reduction))

concentration_custom <- data.frame(PathogenName = c("Rotavirus"),
                                   min = c(1000),
                                   max = c(2000),
                                   type = c("uniform"))

set.seed(2024)
bare <- run_qmra_custom(scenario = sc,
                  concentration = concentration_custom,
                  volume = data.frame(min=c(0.005,0.005), max = c(0.01, 0.01), type = c("triangle", "triangle")),
                  frequency = c(48L, 60L),
                  regulationLog = regulation_reduction,
                  regulationConcentration = regulation_concentration,
                  initialSituation = TRUE) # min 5 mL, max 10 mL

set.seed(2024)
ppe  <- run_qmra_custom(scenario = sc,
                  concentratio = concentration_custom,
                  volume = data.frame(min=c(0.0005,0.0005), max = c(0.001, 0.001), type = c("triangle", "triangle")),
                  frequency = c(48L, 60L),
                  regulationLog = regulation_reduction,
                  regulationConcentration = regulation_concentration,
                  initialSituation = TRUE) # min 0.5 mL min 0.1 mL

As the the other functions run_qmra_* this custom function return 2 plots and to formattable. See d-interpreting-riskfor more detail on the output.

library(ggplot2)
cowplot::plot_grid(
      bare$dalys$Rotavirus +
        labs(subtitle = "min 5 mL, max 10 mL") +
        theme(plot.subtitle = element_text(hjust = 0.5)),
      ppe$dalys$Rotavirus +
        labs(subtitle = "min 0.5 mL min 0.1 mL") +
        theme(plot.subtitle = element_text(hjust = 0.5)),
      ncol = 2,
      align = "h"
    )

The risk tracks the volume almost proportionally: a tenfold cut in swallowed volume gives roughly a tenfold cut in DALYs. More tellingly, the 95th percentile crosses the line — the bare-hand upper tail sits above 1e-6 while the PPE tail drops below it. Deciding whether a scenario meets the target can come down to this single exposure assumption. For how to read these ranges against the regulatory line, see vignette("d-interpreting-risk", package = "ambre"); for why every number is a distribution rather than a point, see vignette("g-monte-carlo-engine", package = "ambre").

Technical notes and caveats