Improve correlation tests. (#287)

The tests for the correlation parameters were using a tolerance of 1e2
everywhere. Since the tolerance is relative, assertions would
effectively succeed for any value within a few orders of magnitude.
Setting obviously wrong expected results in the assertions did not
produce any errors.

This switches the tolerance to 0.1. This was chosen through
experimenting, as something that didn't cause false negative even after
many runs, while also being reasonably close. When using 0.01 I did get
1 failure in 100 runs of the test suite.

Some of the assertions were using incorrected expected values, which
flew under the radar because of the huge tolerance. I've also cleared up
the tests and bit and made them more consitent with one another.

Finally I changed the `CorrelationParameters` constructor to accept only
the values it needs (population size and intervention booleans), rather
than the full simulation parameters. This makes the test a bit more
concise, and will also help with upcoming tests that work at restoring
correlation state while adding interventions. The existing public
wrapper `get_correlation_parameters` still has the same interface as
before.

R/correlation.R

@@ -24,11 +24,11 @@ CorrelationParameters <- R6::R6Class(
   public = list(
 
     #' @description initialise correlation parameters
-    #' @param parameters model parameters
-    initialize = function(parameters) {
-      # Find a list of enabled interventions
-      enabled <- vlapply(INTS, function(name) parameters[[name]])
-      private$interventions <- INTS[enabled]
+    #' @param population popularion size
+    #' @param interventions character vector with the name of enabled interventions
+    initialize = function(population, interventions) {
+      private$population <- population
+      private$interventions <- interventions
 
       # Initialise a rho matrix for our interventions
       n_ints <- private$n_ints()
@@ -38,9 +38,6 @@ CorrelationParameters <- R6::R6Class(
         ncol = n_ints,
         dimnames = list(private$interventions, private$interventions)
       )
-
-      # Store population for mvnorm draws
-      private$population <- parameters$human_population
     },
 
     #' @description Add rho between rounds
@@ -183,7 +180,10 @@ CorrelationParameters <- R6::R6Class(
 #' 
 #' # You can now pass the correlation parameters to the run_simulation function
 get_correlation_parameters <- function(parameters) {
-  CorrelationParameters$new(parameters)
+  # Find a list of enabled interventions
+  enabled <- vlapply(INTS, function(name) parameters[[name]])
+
+  CorrelationParameters$new(parameters$human_population, INTS[enabled])
 }
 
 #' @title Sample a population to intervene in given the correlation parameters

tests/testthat/test-correlation.R

@@ -1,105 +1,136 @@
 test_that('1 correlation between rounds gives sensible samples', {
   pop <- 1e6
   target <- seq(pop)
-  vaccine_coverage <- .2
-  parameters <- get_parameters(list(
-    human_population = pop,
-    pev = TRUE
-  ))
-  correlations <- get_correlation_parameters(parameters)
+
+  coverage_1 <- .2
+  coverage_2 <- .4
+
+  correlations <- CorrelationParameters$new(pop, c('pev'))
   correlations$inter_round_rho('pev', 1)
-  round_1 <- sample_intervention(target, 'pev', vaccine_coverage, correlations)
-  round_2 <- sample_intervention(target, 'pev', vaccine_coverage, correlations)
-  expect_equal(sum(round_1), pop * .2, tolerance=1e2)
-  expect_equal(sum(round_2), pop * .2, tolerance=1e2)
-  expect_equal(sum(round_1 & round_2), pop * .2, tolerance=1e2)
+
+  round_1 <- sample_intervention(target, 'pev', coverage_1, correlations)
+  round_2 <- sample_intervention(target, 'pev', coverage_2, correlations)
+
+  expect_equal(sum(round_1), pop * coverage_1, tolerance=.1)
+  expect_equal(sum(round_2), pop * coverage_2, tolerance=.1)
+
+  expect_equal(
+    sum(round_1 & round_2),
+    pop * min(coverage_1, coverage_2),
+    tolerance=.1)
+
+  expect_equal(
+    sum(round_1 | round_2),
+    pop * max(coverage_1, coverage_2),
+    tolerance=.1)
 })
 
 test_that('0 correlation between rounds gives sensible samples', {
   pop <- 1e6
   target <- seq(pop)
-  vaccine_coverage <- .5
-  parameters <- get_parameters(list(
-    human_population = pop,
-    pev = TRUE
-  ))
-  correlations <- get_correlation_parameters(parameters)
+
+  coverage_1 <- .2
+  coverage_2 <- .4
+
+  correlations <- CorrelationParameters$new(pop, c('pev'))
   correlations$inter_round_rho('pev', 0)
-  round_1 <- sample_intervention(target, 'pev', vaccine_coverage, correlations)
-  round_2 <- sample_intervention(target, 'pev', vaccine_coverage, correlations)
+
+  round_1 <- sample_intervention(target, 'pev', coverage_1, correlations)
+  round_2 <- sample_intervention(target, 'pev', coverage_2, correlations)
+
+  expect_equal(sum(round_1), pop * coverage_1, tolerance=.1)
+  expect_equal(sum(round_2), pop * coverage_2, tolerance=.1)
+
   expect_equal(
-    length(intersect(which(round_1), which(round_2))),
-    pop * .5,
-    tolerance=1e2
-  )
-  expect_equal(sum(round_1), sum(round_2), tolerance=1e2)
-  expect_equal(sum(round_1), pop * .5, tolerance=1e2)
+    sum(round_1 & round_2),
+    pop * coverage_1 * coverage_2,
+    tolerance=.1)
+
+  expect_equal(
+    sum(round_1 | round_2),
+    pop * (coverage_1 + coverage_2 - (coverage_1 * coverage_2)),
+    tolerance=.1)
 })
 
 test_that('1 correlation between interventions gives sensible samples', {
   pop <- 1e6
   target <- seq(pop)
-  vaccine_coverage <- .2
-  mda_coverage <- .2
-  parameters <- get_parameters(list(
-    human_population = pop,
-    pev = TRUE,
-    mda = TRUE
-  ))
-  correlations <- get_correlation_parameters(parameters)
+
+  pev_coverage <- .2
+  mda_coverage <- .4
+
+  correlations <- CorrelationParameters$new(pop, c('pev', 'mda'))
   correlations$inter_round_rho('pev', 1)
   correlations$inter_round_rho('mda', 1)
   correlations$inter_intervention_rho('pev', 'mda', 1)
-  vaccine_sample <- sample_intervention(target, 'pev', vaccine_coverage, correlations)
+
+  pev_sample <- sample_intervention(target, 'pev', pev_coverage, correlations)
   mda_sample <- sample_intervention(target, 'mda', mda_coverage, correlations)
 
-  expect_equal(sum(vaccine_sample), pop * .2, tolerance=1e2)
-  expect_equal(sum(mda_sample), pop * .2, tolerance=1e2)
-  expect_equal(sum(vaccine_sample & mda_sample), pop * .2, tolerance=1e2)
+  expect_equal(sum(pev_sample), pop * pev_coverage, tolerance=.1)
+  expect_equal(sum(mda_sample), pop * mda_coverage, tolerance=.1)
+
+  expect_equal(
+    sum(pev_sample & mda_sample),
+    pop * min(pev_coverage, mda_coverage),
+    tolerance=.1)
+
+  expect_equal(
+    sum(pev_sample | mda_sample),
+    pop * max(pev_coverage, mda_coverage),
+    tolerance=.1)
 })
 
 test_that('0 correlation between interventions gives sensible samples', {
   pop <- 1e6
   target <- seq(pop)
-  vaccine_coverage <- .2
-  mda_coverage <- .2
-  parameters <- get_parameters(list(
-    human_population = pop,
-    pev = TRUE,
-    mda = TRUE
-  ))
-  correlations <- get_correlation_parameters(parameters)
+
+  pev_coverage <- .2
+  mda_coverage <- .4
+
+  correlations <- CorrelationParameters$new(pop, c('pev', 'mda'))
   correlations$inter_round_rho('pev', 1)
   correlations$inter_round_rho('mda', 1)
   correlations$inter_intervention_rho('pev', 'mda', 0)
-  vaccine_sample <- sample_intervention(target, 'pev', vaccine_coverage, correlations)
+
+  pev_sample <- sample_intervention(target, 'pev', pev_coverage, correlations)
   mda_sample <- sample_intervention(target, 'mda', mda_coverage, correlations)
+
+  expect_equal(sum(pev_sample), pop * pev_coverage, tolerance=.1)
+  expect_equal(sum(mda_sample), pop * mda_coverage, tolerance=.1)
+
   expect_equal(
-    length(intersect(which(vaccine_sample), which(mda_sample))),
-    pop * .5,
-    tolerance=1e2
-  )
-  expect_equal(sum(vaccine_sample), sum(mda_sample), tolerance=1e2)
-  expect_equal(sum(vaccine_sample), pop * .5, tolerance=1e2)
+    sum(pev_sample & mda_sample),
+    pop * pev_coverage * mda_coverage,
+    tolerance=.1)
+
+  expect_equal(
+    sum(pev_sample | mda_sample),
+    pop * (pev_coverage + mda_coverage - (pev_coverage * mda_coverage)),
+    tolerance=.1)
 })
 
 test_that('-1 correlation between interventions gives sensible samples', {
   pop <- 1e6
   target <- seq(pop)
-  vaccine_coverage <- .2
-  mda_coverage <- .2
-  parameters <- get_parameters(list(
-    human_population = pop,
-    pev = TRUE,
-    mda = TRUE
-  ))
-  correlations <- get_correlation_parameters(parameters)
+
+  pev_coverage <- .2
+  mda_coverage <- .4
+
+  correlations <- CorrelationParameters$new(pop, c('pev', 'mda'))
   correlations$inter_round_rho('pev', 1)
   correlations$inter_round_rho('mda', 1)
   correlations$inter_intervention_rho('pev', 'mda', -1)
-  vaccine_sample <- sample_intervention(target, 'pev', vaccine_coverage, correlations)
+
+  pev_sample <- sample_intervention(target, 'pev', pev_coverage, correlations)
   mda_sample <- sample_intervention(target, 'mda', mda_coverage, correlations)
-  expect_equal(length(intersect(which(vaccine_sample), which(mda_sample))), 0)
-  expect_equal(sum(vaccine_sample), .2 * pop, tolerance=1e2)
-  expect_equal(sum(mda_sample), .2 * pop, tolerance=1e2)
+
+  expect_equal(sum(pev_sample), pop * pev_coverage, tolerance=.1)
+  expect_equal(sum(mda_sample), pop * mda_coverage, tolerance=.1)
+
+  expect_equal(sum(pev_sample & mda_sample), 0, tolerance=.1)
+  expect_equal(
+    sum(pev_sample | mda_sample),
+    pop * (pev_coverage + mda_coverage),
+    tolerance=.1)
 })