Allow multiple concurrent interventions

[pratikunterwegs]

This issue is to request that epidemics::epi_demic() should be able to include multiple, stackable interventions. This would allow examining and comparing the effects of particular sequences of interventions (e.g. school closures before wider measures, or vice versa).

A template for fixes to this issue could be found in the implementation of the multiple (concurrent) vaccination regimes, including the implementation of a c() method for vaccination objects. See #53 and #52 for details.

As an example:

i1 = intervention(...)
i2 = intervention(...)
all_interventions = c(i1, i2)

epidemic(population, infection, intervention = all_interventions)

One conceptual question is whether overlapping contact reduction for a particular demographic group should be additive (in terms of percentage points), or multiplicative.

[bahadzie]

Given the following code

npi_1 <- intervention(
  time_begin = 30,
  time_end = 60,
  contact_reduction = matrix(0.15, nrow = 3)
)

print(npi_1)
# <intervention>
# Intervention name: NA
# Time begin:
#      [,1]
# [1,]   30
# Time end:
#      [,1]
# [1,]   60
# Contact reduction:
#      [,1]
# [1,] 0.15
# [2,] 0.15
# [3,] 0.15

# second dose regime
npi_2 <- intervention(
  time_begin = 45,
  time_end = 75,
  contact_reduction = matrix(0.1, nrow = 3)
)

print(npi_2)
# <intervention>
# Intervention name: NA
# Time begin:
#      [,1]
# [1,]   45
# Time end:
#      [,1]
# [1,]   75
# Contact reduction:
#      [,1]
# [1,]  0.1
# [2,]  0.1
# [3,]  0.1

multi_npi <- c(npi_1, npi_2)
print(multi_npi)

# <intervention>
# Intervention name: NA
# Time begin:
#      npi_1 npi_2
# [1,]    30    45
# Time end:
#      npi_1 npi_2
# [1,]    60    75
# Contact reduction:
#      npi_1 npi_2
# [1,]  0.15   0.1
# [2,]  0.15   0.1
# [3,]  0.15   0.1

Should multi_npi not be combined only where time_begin and time_end have overlapping days. So using the code above the output should be
Day 30-44: reduction 0.15
Day 45-60: reduction 0.25 (0.15 + 0.1)
Day 61-75: reduction 0.10

and the output of print(multi_npi) be

# <intervention>
# Intervention name: NA
# Time begin:
#      npi_1 npi_2 npi_3
# [1,]    30    46    61
# Time end:
#      npi_1 npi_2 npi_3
# [1,]    45    60    75
# Contact reduction:
#      npi_1 npi_2 npi_3
# [1,]  0.15  0.25   0.1
# [2,]  0.15  0.25   0.1
# [3,]  0.15  0.25   0.1

[pratikunterwegs]

Hi @bahadzie, thanks for taking a look - you're right that the NPIs overlap in the way you describe, and this is taken care of in the background in the code added in PR #65.

epidemics/inst/include/intervention.h

Lines 17 to 40 in c33d6f5

	inline Eigen::ArrayXd cumulative_intervention(const double &t,
	const Rcpp::List &intervention) {
	// copy the list elements
	// TODO(all): reduce copying
	Rcpp::NumericVector time_begin = intervention["time_begin"];
	Rcpp::NumericVector time_end = intervention["time_end"];
	Rcpp::NumericMatrix cr = intervention["contact_reduction"];
	Rcpp::NumericVector eff_con_red(cr.size());
	// iterate over the Rcpp vector time_begin, a member of the list-like
	// intervention-class
	for (size_t i = 0; i < time_begin.size(); i++) {
	if (t >= time_begin(i) && t <= time_end(i)) {
	eff_con_red += cr(Rcpp::_, i);
	}
	}
	Eigen::ArrayXd effective_contact_reduction(
	Rcpp::as<Eigen::ArrayXd>(eff_con_red));
	return effective_contact_reduction;
	}

What I thought would be more useful and easier for users to understand is to treat each NPI as a single unit, and handle the concatenation internally. I thought it might be confusing for users to find a third NPI if they have not really specified one. Happy to hear how policy folks are used to thinking about this from @BlackEdder and @TimTaylor et al.

[TimTaylor]

@pratikunterwegs No actual use to comment on but as a visual way of presenting, the current c() method seems nice and readable. With a disclaimer that I've not looked in to the intervention class the only thing I'd wonder about is whether the c() method is needed over just accepting a list of interventions.

You also have an option to provide a + operator method which could act in the way @bahadzie describes should you wish (again I'm not sure if that is needed).

[pratikunterwegs]

Thanks @TimTaylor - my question was perhaps more on the lines of, if I have two NPIs overlapping, is the overlap considered a distinct intervention of its own? Perhaps this taxonomy hasn't been worked out yet in public health?

Representing the start and end times in a single intervention object makes it really easy to iterate over the rows of the time and contact reduction matrices in the C++ code, I thought - so went for this rather than a list of interventions. The + operator idea is interesting, but the result would be the same as c() wouldn't it?

[TimTaylor]

I see. Yes I was getting confused. I feel like I could be persuaded in multiple directions here so will dwell on it a little more.