Bug in seed infections for intervention scenarios

[innovative-simulator]

In UK.R simulations for intervention scenarios are being performed with the default seed_times (1 person on day 1), not the seed_times used for the base scenario (2 per day, for 28 days, starting at the randomly sampled seed_start, as per the Lancet paper).

To check, see UK.R:

• Line 372: For the Base scenario, the parameter set named "params" is set to default values by duplicating "parameters". (The default value for variable “seed_times” is set to 1 person on day 1, at line 237 of script “covidm_params.R”.)
• Lines 376, 377: seed_times and dist_seed_ages are generated and stored in params. For seed_times, the values equate to 2 people per day for 28 days, starting at the previously defined seed_start.
• Line 421: The simulation is run with this params for the “Base” scenario.
• Line 443: For each intervention scenario, params is re-set to default values by duplicating parameters (as in line 372).
• No new attempt is made in this or any dependent script to alter the values of seed_times and dist_seed_ages.
• Line 490: The simulation is run with params for the intervention scenario, with seed_times and dist_seed_ages still set to the default, not to the values used for the Base scenario.

Running with 1 seed infection, instead of 2*28, means many counties fail to have an epidemic. This is difficult to see from the aggregated output time series in the X-dynamics.qs file, but try:

    library(data.table)
    library(qs)
    D <- qread("1-dynamics.qs") # Assuming this file is in the working directory.
    dim(D[compartment=="E" & scenario !="Base" & t<=0 & value != 0]) # Returns 0 rows.
    dim(D[compartment=="E" & scenario !="Base" & t==1 & value == 0]) # Returns 0 rows.
    dim(D[compartment=="E" & scenario =="Base" & t<=0 & value != 0]) # Returns lots of rows.
    dim(D[compartment=="E" & scenario =="Base" & t==1 & value == 0]) # Returns lots of rows.

The effects of the bug vary between runs, but in general are that interventions produce much bigger reductions in numbers of cases, beds, and deaths than they would if they received the same seed infections as the Base scenario.

The authors have been notified, so a bug fix should follow.

In the meantime, if you are using UK.R try the following correction. (I have concatenated the lines of code to preserve the line numbering.):

Line 372:
params = duplicate(parameters); initial_seed_times <- list(); initial_dist_seed_ages <- list(); # CW's addition: Setup lists to store populations' seed_times and dist_seed_ages
Lines: 376, 377:

        initial_seed_times[[j]] = rep(seed_start[j] + 0:27, each = 2); params$pop[[j]]$seed_times = initial_seed_times[[j]]; # CW's addition: Store seed_times. Then use them.
        initial_dist_seed_ages[[j]] = cm_age_coefficients(25, 50, 5 * 0:16); params$pop[[j]]$dist_seed_ages = initial_dist_seed_ages[[j]]; # CW's addition: Store dist_seed_ages. Then use them.

Line 446:
params$pop[[j]]$y = covy; params$pop[[j]]$seed_times = initial_seed_times[[j]]; params$pop[[j]]$dist_seed_ages = initial_dist_seed_ages[[j]]; # CW's addition: Re-use Base scenario's stored seed_times and dist_seed_ages.

[nicholasdavies]

Thanks for pointing this out — we are preparing a fix and correction to the article.

We don't expect this issue to affect our qualitative conclusions, specifically in terms of the relative impact of different interventions, and it has no impact upon the calculation of R values under different interventions, meaning that our conclusions regarding which interventions may be effective for controlling the UK epidemic will not change. This error occurred while preparing the manuscript for Lancet ID, so did not affect the modelling evidence originally submitted to government advisory bodies.

Thanks for submitting the feedback!

Nick