Hacks for starting with a partially immune population

[pratikunterwegs]

This is a small reprex that shows how to start with a partially immune population.

In PR #118 I wasn't clear that since the vaccination parameters are taken from the value passed to the vaccine_investment argument, passing vaccine_investment = NULL (the default) results in the (internal) use of dummy_vaccination() which specifies a $\psi$ (waning rate) of 1.0 - so everybody is shifted out of the vaccinated compartment pretty quickly, leading to very small differences in the epidemic size.

I think that with these bespoke scenarios it's best to define a <daedalus_vaccination> manually, passing a zero vaccination rate $\nu$ with $\psi > 1$ (a waning period of 6 months was often assumed during Covid), to keep initially immune individuals in the initially immune compartment.

library("daedalus.data")
library("daedalus")
library(data.table)
library(ggplot2)
library(scales)

# use england = GB to hide data
england <- daedalus_country("GB")

# note: this is not actually being used?
# assign to england hosp capacity if you want state-dependent responses
# to trigger on this
response_threshold <- 1000

# small function to get initially immune
f = function(x) {
  d = get_data(x)
  d = d[d$vaccine_group == "vaccinated" & d$time == 0, ]
  sum(d$value)
}

# pass a custom vaccination object
# this is needed to specify the parameter `psi`, the waning rate
# the default value of `vaccine_investment` is `NULL`, which assumes a `psi`
# of 1.0, which moves everybody out of the vaccinated compartment
vax = daedalus_vaccination(
  "low",
  england,
  start_time = 600,
  0.0,
  0.0,
  50,
  waning_period = 180
)

# iterate over initial immunity levels
init_immune_levels = c(0.9, 0.6, 0.3, 0.1, 0.0)
data_list <- lapply(
  init_immune_levels,
  function(x) {
    daedalus(
      england,
      daedalus_infection("sars_cov_2_omicron", r0 = 3.7),
      vaccine_investment = vax,
      response_strategy = "none",
      initial_state_manual = list(
        p_immune = x
      )
    )
  }
)
names(data_list) = sprintf("immune %.2f", init_immune_levels)

# get the total number initially immune
# this looks okay, reduces as initially immune reduces
lapply(data_list, f)
#> $`immune 0.90`
#> [1] 61097404
#> 
#> $`immune 0.60`
#> [1] 40731602
#> 
#> $`immune 0.30`
#> [1] 20365801
#> 
#> $`immune 0.10`
#> [1] 6788600
#> 
#> $`immune 0.00`
#> [1] 0

# check epidemic sizes and total hospitalisations
lapply(data_list, get_epidemic_summary)
#> $`immune 0.90`
#>         value                measure
#> 1    553036.8           total_deaths
#> 2 131403644.5          epidemic_size
#> 3   1801826.2 total_hospitalisations
#> 
#> $`immune 0.60`
#>         value                measure
#> 1    565893.7           total_deaths
#> 2 133748252.1          epidemic_size
#> 3   1837790.9 total_hospitalisations
#> 
#> $`immune 0.30`
#>         value                measure
#> 1    646556.3           total_deaths
#> 2 135706211.3          epidemic_size
#> 3   1867425.4 total_hospitalisations
#> 
#> $`immune 0.10`
#>         value                measure
#> 1    653670.1           total_deaths
#> 2 136900509.9          epidemic_size
#> 3   1885831.3 total_hospitalisations
#> 
#> $`immune 0.00`
#>         value                measure
#> 1    657056.1           total_deaths
#> 2 137471379.5          epidemic_size
#> 3   1894598.0 total_hospitalisations

#### get hospital occupancy ####
data_timeseries = Map(
  data_list,
  names(data_list),
  f = function(x, y) {
    d = get_data(x)
    d$scenario = y

    d
  }
)

data <- rbindlist(data_timeseries)

# sum over age and econ strata as total is more relevant
data_hosp <- data[
  compartment == "hospitalised",
  .(value = sum(value)),
  by = c("time", "scenario")
]

# check actual outcomes of interest - these don't look as good
ggplot(data_hosp) +
  geom_line(aes(time, value, colour = scenario)) +
  geom_hline(yintercept = response_threshold, linetype = "dotted")

^{Created on 2025-09-05 with reprex v2.1.1}

[pratikunterwegs]

@pabloperguz and @patcatgit - here's a presumably mid-pandemic/pre-existing immunity modelling use case from UKHSA, with my recommended solution.

The main avenue for development is editing the model to have three (or a flexible number?) of immunisation strata, as passing a custom <daedalus_vaccination> prevents to specify pre-existing immunity disallows specifying a vaccination regime that kicks in later.