nee_and_composite_scheme_intervals.qmd

---
title: Calculate composite intervals from selected schemes' final report params
author:
  name: Fran Barton
  email: francis.barton@nhs.net
date: 30 January 2026
lang: en-GB
theme: lumen
embed-resources: true
minimal: true
fig-dpi: 96
fig-responsive: true
lightbox: true
out-width: 92%
knitr:
  opts_chunk:
    dev: ragg_png
    message: false
execute:
  warning: true
df-print: kable
lightbox: true
mainfont: Fira Sans
---

## Introduction

This document generates a set of "average mitigation assumptions" for each of
96 Types of Potentially Mitigatable Activity (TPMA).

The results are shown in the final table of the document.
They are also made available in various formats as xlsx and rds files.

The document also shows some intermediate data objects for reference, and
documents the assumptions made in the code in order to arrive at the final
table, for checking and transparency.

## Assumptions and decisions made about the process

The assumptions made are, in summary:

1. Only model runs a 'final' run stage are used, and maximum of 1 run per scheme
2. Only `activity_avoidance` and `efficiencies` intervals are included
3. All TPMAs are included, by default, with the exception of a small number of
  TPMAs labelled `bads_*`, which are excluded due to not matching the list
  provided in a reference `mitigators.json` file.
4. Intervals for all types of activity (inpatient, outpatient, aae) are
  included.
5. `p10`, and `p90` values are reported, taken from the distributions created
  from all interval values available for each TPMA, across all selected runs.

Where these assumptions are not in line with what is required, it is simple to
make changes to the process.
Please contact Fran Barton with details of any amendments, or any questions.

## Data process

### Initial data

We retrieve the params data for all tagged runs, from a `pins` board, and a
csv file with a list of "mitigators" and their names.

We also retrieve a set of NEE (national expert elicitation) values for 78 TPMAs.

```{r get pinned params data}
#| cache: true
#| filename: "get pinned params data"
board <- pins::board_connect(auth = "envvar")

nhp_tagged_runs_meta <- board |>
  pins::pin_read("matt.dray/nhp_tagged_runs_meta")

params_data_init <- board |>
  pins::pin_read("matt.dray/nhp_tagged_runs_params")

tpma_lookup <- azkit::get_container("supporting-data") |>
  azkit::read_azure_csv("mitigator-lookup.csv") |> View()
  dplyr::rename(
    change_factor = "mitigator_type",
    strategy = "mitigator_variable",
    strategy_name = "mitigator_name"
  )
```

```{r get nee values}
#| filename: "get nee values"
#| cache: true

nee_intervals <- azkit::get_container("supporting-data") |>
  AzureStor::storage_load_rds("nee_table.rds", type = "none") |>
  # azkit::read_azure_rds("nee_table.rds") |>
  dplyr::rename_with(stringr::str_to_snake) |>
  dplyr::select(c(
    change_factor = "type",
    type = "strategy_type",
    strategy = "param_name",
    # strategy_name = "strategy_label",
    p10 = "lower_ci",
    p90 = "upper_ci"
  )) |>
  dplyr::mutate(
    dplyr::across("type", \(x) {
      dplyr::case_when(
        grepl("^inpatient", x) ~ "ip",
        grepl("^length of stay", x) ~ "ip",
        grepl("^outpatient", x) ~ "op",
        grepl("^A&E", x) ~ "aae",
        .default = x
      )
    }),
    dplyr::across(c("p10", "p90"), \(x) x / 100)
  )

```

We begin with data for `r length(params_data_init)` runs in `params_data_init`.

`tpma_lookup` provides a list of `r nrow(tpma_lookup)` "strategies".
Later we will use this to ensure that only intervals for strategies included
in this list are used in the final data output.

```{r filter to preferred params per dataset}
#| filename: "filter to preferred params per dataset"
run_stages <- c("final_report_ndg2", "final_report_ndg3", "final_report_ndg1")

preferred_runs <- nhp_tagged_runs_meta |>
  dplyr::filter(dplyr::if_any("run_stage", \(x) x %in% {{ run_stages }})) |>
  dplyr::slice_min(
    match(.data[["run_stage"]], {{ run_stages }}),
    by = "dataset"
  ) |>
  dplyr::select(!"file")

preferred_runs
```

The initial set of `r nrow(nhp_tagged_runs_meta)` is now filtered down to
`r nrow(preferred_runs)` "preferred" runs.

### Summary of selected runs to be included

```{r filter params data to preferred runs per dataset only}
#| filename: "filter params data to preferred runs per dataset only"
preferred_run_names <- preferred_runs |>
  glue::glue_data("{dataset}_{scenario}_{create_datetime}")

params_data <- params_data_init |>
  purrr::keep_at(preferred_run_names)

# check that all names were found
length(params_data) == length(preferred_run_names)

core_elements <- c(
  "dataset",
  "scenario",
  "start_year",
  "end_year",
  "app_version",
  "create_datetime"
)
params_data_tbl_summary <- params_data |>
  purrr::map(\(x) purrr::keep_at(x, core_elements)) |>
  purrr::map(tibble::as_tibble_row) |>
  purrr::list_rbind() |>
  dplyr::mutate(dplyr::across(tidyselect::ends_with("year"), as.integer))

```

This table is not used in the further calculation below, but is shown here for
information.

```{r display summary table}
#| filename: "display summary table"
params_data_tbl_summary |>
  dplyr::arrange(.data[["create_datetime"]])
```

### Extract intervals from selected params

Here we take the `activity_avoidance` and `efficiencies` elements of the params
data only.
For each of the `r nrow(preferred_runs)` runs, and for each "strategy", we
extract the provided `p10` and `p90` values and pivot these into columns in a
data frame.

We then filter `combined_intervals_tbl` against the `tpma_lookup` table that was
read in at the top of the document.
Intervals are only retained if both their `change_factor` (activity_avoidance or
efficiencies) and `strategy` fields are matched.

```{r extract intervals}
#| filename: "extract intervals"
create_intervals_tbl <- function(lst) {
  lst |>
    purrr::map_depth(2, \(x) list(x[["interval"]])) |>
    purrr::imap(\(x, y) dplyr::bind_cols(type = y, tibble::as_tibble_row(x))) |>
    purrr::list_rbind() |>
    tidyr::pivot_longer(
      !"type",
      names_to = "strategy",
      values_drop_na = TRUE
    ) |>
    tidyr::unnest_wider("value", names_sep = "_") |>
    dplyr::rename(p10 = "value_1", p90 = "value_2")
}


combined_intervals_tbl <- params_data |>
  purrr::map(\(x) purrr::keep_at(x, c("activity_avoidance", "efficiencies"))) |>
  purrr::map_depth(2, create_intervals_tbl) |>
  purrr::map(\(x) purrr::list_rbind(x, names_to = "change_factor")) |>
  purrr::list_rbind(names_to = "run_label") |>
  # only include "approved" list of strategies (excludes "bads_*"))
  dplyr::inner_join(tpma_lookup, c("change_factor", "strategy")) |>
  dplyr::relocate("strategy_name", .after = "strategy")

```

The `combined_intervals_tbl` data frame contains
`r nrow(combined_intervals_tbl)` rows.


### Visualising scheme-supplied intervals

This will help us notice supplied values that don't look right, and that may
cause problems when creating the distributions (later) from which we will create
our composite values.

```{r lollipop plot functions}
#| filename: "lollipop plot functions"
plot_lollipop <- function(dat, strat) {
  c1 <- "violetred1"
  c2 <- "cadetblue"
  dat |>
    dplyr::arrange(.data[["p10"]]) |>
    dplyr::mutate(dplyr::across("scheme", forcats::fct_inorder)) |>
    ggplot2::ggplot(ggplot2::aes(y = .data[["scheme"]])) +
    ggplot2::geom_point(ggplot2::aes(x = .data[["p10"]]), colour = c1) +
    ggplot2::geom_point(ggplot2::aes(x = .data[["p90"]]), colour = c2) +
    ggplot2::scale_x_continuous(limits = c(0, 1)) +
    ggplot2::geom_segment(
      ggplot2::aes(x = .data[["p10"]], xend = .data[["p90"]]),
      colour = "grey67"
    ) +
    ggplot2::labs(
      title = strat,
      x = "scheme-supplied p10 (red) and p90 (blue) values"
    ) +
    ggplot2::theme(title = ggplot2::element_text(size = 8))
}
create_patchwork <- function(dat, type) {
  type_full <- dplyr::case_match(
    type,
    "ip" ~ "Inpatient",
    "op" ~ "Outpatient",
    "aae" ~ "A&E"
  )
  n_cols <- dplyr::if_else(type == "ip", 12, 8)
  dat |>
    tidyr::nest(.by = c("strategy_label", "n_schemes")) |>
    dplyr::arrange(.data[["n_schemes"]]) |>
    dplyr::select(!"n_schemes") |>
    tibble::deframe() |>
    purrr::imap(plot_lollipop) |>
    patchwork::wrap_plots(
      ncol = n_cols,
      axes = "collect_x",
      axis_titles = "collect"
    ) +
    patchwork::plot_annotation(type_full)
}

patchworks <- combined_intervals_tbl |>
  dplyr::filter(dplyr::if_all(c("p10", "p90"), \(x) dplyr::between(x, 0, 1))) |>
  dplyr::mutate(
    scheme = stringr::str_extract(.data[["run_label"]], "^[[:alnum:]]{3}"),
    cf = dplyr::if_else(grepl("^eff", .data[["change_factor"]]), "eff", "aa"),
    .keep = "unused"
  ) |>
  dplyr::add_count(.data[["strategy"]], name = "n_schemes") |>
  dplyr::mutate(
    strategy_label = stringr::str_wrap(
      gsub("_", " ", glue::glue("{cf}: {strategy} ({n_schemes} schemes)")),
      32
    )
  ) |>
  tidyr::nest(.by = "type") |>
  tibble::deframe() |>
  purrr::imap(create_patchwork)

```

```{r inpatient validation plot}
#| fig-height: 15
#| fig-width: 27
#| echo: false

patchworks[["ip"]]

```


```{r outpatient validation plot}
#| fig-height: 6
#| fig-width: 18
#| echo: false

patchworks[["op"]]

```

```{r aae validation plot}
#| fig-height: 6
#| fig-width: 18
#| echo: false

patchworks[["aae"]]

```


```{r tidy intervals summary}
#| eval: false
combined_intervals_tbl |>
  dplyr::arrange(.data[["p10"]]) |>
  openxlsx2::write_xlsx("scheme_supplied_intervals.xlsx")

```

### Remove rows with inapplicable values

```{r remove rows with inapplicable values}
#| filename: "remove rows with inapplicable values"

# 12 rows were identified through inspection to contain inapplicable values
remove_rows <- combined_intervals_tbl |>
  dplyr::select(c("run_label", "strategy")) |>
  dplyr::filter(
    (grepl("^RH5", .data[["run_label"]]) &
      .data[["strategy"]] == "virtual_wards_activity_avoidance_ari") |
      (grepl("^RAS", .data[["run_label"]]) &
        (grepl("adult_non-surgical$", .data[["strategy"]]) |
          .data[["strategy"]] == "pre-op_los_1-day")) |
      (grepl("^RGN", .data[["run_label"]]) &
        grepl("^discharged_no_treatment_.*_ambulance$", .data[["strategy"]])) |
      (grepl("^RTX", .data[["run_label"]]) &
        .data[["strategy"]] %in%
          c(
            "stroke_early_supported_discharge",
            "virtual_wards_efficiencies_ari",
            "raid_ip",
            "virtual_wards_efficiencies_heart_failure"
          )
      )
  )

nrow(remove_rows) == 12 # check

# We also remove any rows where p10 or p90 are not between 0 and 1
tidy_intervals_tbl <- combined_intervals_tbl |>
  dplyr::anti_join(remove_rows, c("run_label", "strategy")) |>
  # values greater than 1 should have already been removed in the above step,
  # but just to make sure:
  dplyr::filter(dplyr::if_all(c("p10", "p90"), \(x) dplyr::between(x, 0, 1)))

nrow(combined_intervals_tbl) - nrow(tidy_intervals_tbl) == 12 # check

```


### Standardise intervals

Because the provided intervals (`p10` and `p90` values) were supplied as
projections of activity avoidance / efficiencies over different time spans (for
example, some were 2019-2036, others 2023-2041), we need to standardise each set
of intervals to an annualised rate, before creating the mixture distributions,
and then multiply up our results to a standard time period.


```{r annualise scheme intervals}
#| filename: "annualise scheme intervals"
run_lookup <- params_data_tbl_summary |>
  dplyr::select(!"app_version") |>
  tidyr::unite(run_label, c("dataset", "scenario", "create_datetime"))

annualise <- \(x, p) `^`(x, (1 / p))

annualised_intervals_tbl <- tidy_intervals_tbl |>
  dplyr::left_join(run_lookup, "run_label") |>
  dplyr::mutate(
    period = .data[["end_year"]] - .data[["start_year"]],
    .keep = "unused"
  ) |>
  dplyr::mutate(
    dplyr::across(c("p10", "p90"), \(x) annualise(x, .data[["period"]])),
    .keep = "unused"
  )

```


```{r annualise nee intervals}
#| filename: "annualise nee intervals"

nee_period <- 20 # years between baseline and horizon

annualised_nee_intervals <- nee_intervals |>
  dplyr::mutate(dplyr::across(c("p10", "p90"), \(x) annualise(x, nee_period)))

```

Now we use a function to generate a distribution across all `p10` and `p90`
values for each strategy.
We then apply a "reinflation" function that converts the annualised p10 and p90
values back to quantities that would apply over, in this case, 5 and 7 years.
(This raises the annualised values to the power of 5 and 7 respectively).

We then report back the `p10`, `p50` and `p90` values for each of those
reinflated distributions.

```{r create_mixtures function}
#| filename: "create_mixtures function"
create_mixtures <- function(intervals_tbl, colnames_vec) {
  intervals_tbl |>
    dplyr::mutate(
      mu = (.data[["p10"]] + .data[["p90"]]) / 2,
      sigma = (.data[["p90"]] - .data[["mu"]]) /
        qnorm(p = 0.90, mean = 0, sd = 1),
      dplyr::across("sigma", \(x) dplyr::if_else(x == 0, 0.0001, x)),
      dist = purrr::map2(.data[["mu"]], .data[["sigma"]], \(m, s) {
        distr::Truncate(distr::Norm(mean = m, sd = s), lower = 0, upper = 1)
      })
    ) |>
    dplyr::summarise(
      mixture = purrr::map(list(.data[["dist"]]), \(l) {
        distr::UnivarMixingDistribution(Dlist = l)
      }),
      .by = tidyselect::all_of(colnames_vec)
    )
}

```

```{r create mixtures across all scheme runs}
#| cache: true
#| filename: "create mixtures across all scheme runs"

scheme_mixtures_tbl <- annualised_intervals_tbl |>
  create_mixtures(c("type", "change_factor", "strategy"))

```

```{r prepare distribution visualisations for scheme mixtures}
#| filename: "prepare distribution visualisations for scheme mixtures"
#| cache: true
add_pdf_values <- function(mixture_dist, x_length = 1001) {
  tibble::tibble(
    x_values = seq(0, 1, length.out = x_length),
    pdf_values = distr::d(mixture_dist)(x_values)
  )
}

mixture_dists_plotting_data <- scheme_mixtures_tbl |>
  dplyr::mutate(
    pdf_data = purrr::map(.data[["mixture"]], add_pdf_values),
    .keep = "unused"
  ) |>
  tidyr::unnest_wider("pdf_data") |>
  tidyr::unnest_longer(c("x_values", "pdf_values"))
```

```{r generate distribution plots}
#| filename: "generate distribution plots"
#| fig-height: 36
mixture_dists_plotting_data |>
  dplyr::mutate(
    strategy_label = stringr::str_wrap(gsub("_", " ", .data[["strategy"]]), 20)
  ) |>
  ggplot2::ggplot(ggplot2::aes(.data[["x_values"]], .data[["pdf_values"]])) +
  ggplot2::geom_line() +
  ggplot2::facet_wrap(~strategy_label, ncol = 4, scales = "free") +
  ggplot2::theme(strip.text = ggplot2::element_text(size = 8))

```

A table to show the `p10`, `p50` and `p90` values taken from the scheme mixture
distribution for each strategy:

```{r pull out intervals from scheme mixtures}
#| filename: "pull out intervals from scheme mixtures"

scheme_intervals_full <- scheme_mixtures_tbl |>
  dplyr::mutate(
    p10 = purrr::map_dbl(.data[["mixture"]], \(m) m@q(p = 0.1)),
    p50 = purrr::map_dbl(.data[["mixture"]], \(m) m@q(p = 0.5)),
    p90 = purrr::map_dbl(.data[["mixture"]], \(m) m@q(p = 0.9)),
    .keep = "unused"
  )

scheme_intervals_full

```

```{r reinflate values}
#| filename: "reinflate values"

scheme_intervals <- scheme_intervals_full
# dplyr::select(!"p50")

init_joined_table <- list(scheme_intervals, annualised_nee_intervals) |>
  rlang::set_names(c("schemes", "nee")) |>
  dplyr::bind_rows(.id = "source")

reinflate_intervals <- function(exp, tbl) {
  dplyr::mutate(tbl, dplyr::across(c("p10", "p50", "p90"), \(x) x^exp))
}

inflation_periods <- c(5, 7) # years to horizon

final_joined_table <- inflation_periods |>
  purrr::map(\(x) reinflate_intervals(x, init_joined_table)) |>
  rlang::set_names(paste0(inflation_periods, "yr")) |>
  purrr::list_rbind(names_to = "period") |>
  # pivot_wider creates NA values in the table, where we are missing a NEE value
  tidyr::pivot_wider(
    names_from = c("source", "period"),
    names_glue = "{source}_{period}_{.value}",
    values_from = c("p10", "p50", "p90")
  ) |>
  tidyr::pivot_longer(
    tidyselect::matches("p[159]0$"),
    names_to = c("source", "period", "name"),
    names_sep = "_",
    cols_vary = "slowest"
  ) |>
  dplyr::mutate(dplyr::across("type", \(x) {
    forcats::fct(x, c("ip", "op", "aae"))
  })) |>
  dplyr::arrange(dplyr::pick(c(
    "type",
    "change_factor",
    "strategy",
    "period",
    "source"
  ))) |>
  # fill missing NEE values with scheme values as a fallback
  tidyr::fill(
    "value",
    .by = c("type", "change_factor", "strategy", "period", "name"),
    .direction = "up" # "nee" rows will be above "schemes" rows per group
  ) |>
  tidyr::pivot_wider()

```


```{r display final joined table}
#| filename: "display final joined table"
final_joined_table |>
  dplyr::mutate(dplyr::across(c("p10", "p50", "p90"), \(x) round(x, 3)))
```

```{r export data}
#| echo: false
#| eval: false

final_joined_table |>
  dplyr::select(!"p50") |>
  readr::write_rds("nee_and_composite_scheme_intervals.rds")

```