Episode suggestion: disease burden model

[adamkucharski]

It would be useful to have an episode illustrating how new infections in {epidemics} translate into hospitalisations and stock-and-flow issues with hospital capacity (i.e. people can remain in a bed for multiple days. Key learning objectives could include how models of transmission can be separated from models of burden if hospitalised cases contribute little to transmission (as well as mentioning where this assumption breaks down, e.g. Ebola).

Could also cover difference between age distribution of infections vs severe disease, and delayed outcomes (linking with this howto: https://epiverse-trace.github.io/howto/analyses/reconstruct_transmission/estimate_infections.html)

Could also be nice opportunity to show how to get parameters in/out of epiparameter.

Some rough illustrative code:

# Simulate some random new cases
n_days <- 200 # Number of days in the simulation
new_cases <- rpois(n_days, lambda = 10) # Random daily new cases (Poisson distributed)
new_cases[101:200] <- 0 # End outbreak after day 100

# Define delay parameters
onset_to_admission <- epiparameter(
  disease = "Disease X",
  epi_name = "onset to admission",
  prob_distribution = create_prob_distribution(
    prob_distribution = "gamma",
    prob_distribution_params = c(shape = 5, scale = 4),
    discretise = T
  )
)

admission_to_discharge <- epiparameter(
  disease = "Disease X",
  epi_name = "admission to discharge",
  prob_distribution = create_prob_distribution(
    prob_distribution = "gamma",
    prob_distribution_params = c(shape = 10, scale = 2),
    discretise = T
  )
)

ihr <- 0.1 # Infection-hospitalisation ratio

# Simulate hospitalizations and discharges
simulate_hospitalizations <- function(new_cases, onset_to_admission, admission_to_discharge, ihr) {
  hosp_day <- c()    # Track daily hospitalizations
  discharge_day <- c() # Track daily discharges
  
  for (day in seq_along(new_cases)) {
    # Hospitalizations from current day's new cases
    n_cases_today <- rbinom(1, new_cases[day], ihr)
    if (n_cases_today > 0) {
      hosp_delays <- generate(onset_to_admission, n_cases_today)
      hosp_day <- c(hosp_day, day + hosp_delays)
    }
  }
  
  # Calculate discharges based on hospitalizations
  for (hosp_date in hosp_day) {
    disch_delays <- generate(admission_to_discharge, 1)
    discharge_day <- c(discharge_day, hosp_date + disch_delays)
  }
  
  list(hosp_day = hosp_day, discharge_day = discharge_day)
}

# Run the simulation
sim_results <- simulate_hospitalizations(new_cases, onset_to_admission, admission_to_discharge, ihr)

# Calculate daily hospitalizations and discharges
hospitalizations <- table(factor(sim_results$hosp_day, levels = 1:n_days))
discharges <- table(factor(sim_results$discharge_day, levels = 1:n_days))

# Total in hospital calculation
in_hospital <- cumsum(hospitalizations) - cumsum(discharges)

# Plot results
plot(new_cases,type="l",ylim=c(0,30))
lines(as.numeric(hospitalizations),col="red")
lines(as.numeric(in_hospital),col="blue")

[amanda-minter]

Thanks for the code - very useful to get things started.

The output from new_infections() in {epidemics} is continuous so we can't use the binomial distribution to generate the hospitalizations as currently implemented in your code

n_cases_today <- rbinom(1, new_cases[day], ihr)

Is there a method you have used before to convert cases from continuous to discrete for this purpose?

[adamkucharski]

I suppose if the number of infections is large, we could use the normal approximation to the binomial distribution?

[amanda-minter]

The number of new infections early in the outbreak is small but the issue of continuous cases also arises with simulating the number of hospitalization delay times here:

hosp_delays <- generate(onset_to_admission, n_cases_today)

we need a discrete number of cases here otherwise we would try to simulate a delay time for a proportion of a case

[adamkucharski]

The alternative would be to do everything with expectations (e.g. multiply continuous value by the continuous distributional value). Actually, perhaps this is cleaner, because the above code has a stochastic element, so we'd ideally simulate more than once...

[amanda-minter]

Something like :

# Simulate some random new cases
n_days <- 200 # Number of days in the simulation
new_cases <- rpois(n_days, lambda = 10) # Random daily new cases (Poisson distributed)
new_cases[101:200] <- 0 # End outbreak after day 100

time <- seq(1,200, by = 1)
# hospitalised cases = new cases * ihr
hosp <- new_cases * ihr 
# time of hospitalization
time_hosp <- time + (5*4) # mean time to hosp from gamma dist
# time of discharge
time_disch <- time_hosp + (10*2) # mean time to discharge from gamma dist

results_hosp <- data.frame(time = time_hosp, number = hosp)
# assume same number are discharged
results_disch <- data.frame(time = time_disch, number = hosp) 

# fill empty times with 0s to calculate cumulative sum
results_hosp <- complete(results_hosp, time = seq(1, max(time_disch)), 
                         fill = list(number = 0))
results_disch <- complete(results_disch, time = seq(1, max(time_disch)), 
                         fill = list(number = 0))

# total in hospital calculation
in_hospital <- cumsum(results_hosp$number) - cumsum(results_disch$number) 
  
plot(time, new_cases, type = "l", ylim = c(0, 30))
lines(results_hosp$time, results_hosp$number, col = "red")
lines(in_hospital, col = "blue")

[adamkucharski]

Thanks. It would be good to represent the distributional spread rather than just the mean (e.g. if 100 people get ill today, it would be useful to plot the expected distribution of admissions over the coming weeks, according to onset_to_admission)

Below is quickest way I could think to code this up - would need to explain to learner what a convolution is (but perhaps not a bad thing, given it's so fundamental to modelling, and methods like this act by doing de-convolution: https://epiverse-trace.github.io/howto/analyses/reconstruct_transmission/estimate_infections.html):

# Simulate some random new cases
n_days <- 200 # Number of days in the simulation
new_cases <- rpois(n_days, lambda = 10) # Random daily new cases (Poisson distributed)
new_cases[101:200] <- 0 # End outbreak after day 100

# Define delay parameters
onset_to_admission <- epiparameter(
  disease = "Disease X",
  epi_name = "onset to admission",
  prob_distribution = create_prob_distribution(
    prob_distribution = "gamma",
    prob_distribution_params = c(shape = 5, scale = 4),
    discretise = T
  )
)

admission_to_discharge <- epiparameter(
  disease = "Disease X",
  epi_name = "admission to discharge",
  prob_distribution = create_prob_distribution(
    prob_distribution = "gamma",
    prob_distribution_params = c(shape = 10, scale = 2),
    discretise = T
  )
)

ihr <- 0.1 # Infection-hospitalisation ratio

# Calculate expected numbers with convolution:
hosp <- new_cases * ihr 
hospitalizations <- convolve(hosp, rev(density(onset_to_admission,0:50)), type = "open")[1:length(hosp)]
discharges <- convolve(hospitalizations, rev(density(admission_to_discharge,0:50)), type = "open")[1:length(hospitalizations)]
in_hospital <- cumsum(hospitalizations) - cumsum(discharges)


plot(new_cases,type="l",ylim=c(0,30))
lines(as.numeric(hospitalizations),col="red")
lines(as.numeric(in_hospital),col="blue")

[adamkucharski]

Also, would need to tweak the hard coded 50 - could use quantile(onset_to_admission,0.999) to define the tail of the delay (which is 59 in this case)?