Contact tracing behaviour

[wfmackey]

Using the various *contact_tracing* arguments, I'm not getting the changes to the proportion of infected people isolated as I had imagined. @HughParsonage can you tell me if I am interpreting/using them incorrectly?

library(tidyverse)
library(purrr)
library(data.table)
#> 
#> Attaching package: 'data.table'
#> The following objects are masked from 'package:dplyr':
#> 
#>     between, first, last
#> The following object is masked from 'package:purrr':
#> 
#>     transpose
library(covid19.model.sa2)

repeat_simulation <- function(repeat_times = 1, 
                              .days = 20, 
                              .PolicyPars = set_policypars(),
                              .EpiPars = set_epipars(),
                              ...) {

  return_simulation <- function(x) {
    d <- 
    simulate_sa2(days_to_simulate = .days, 
                 returner = 3,
                 PolicyPars = .PolicyPars,
                 EpiPars = .EpiPars)

    Status12 <- d[["Status12"]]

    rep_each <- function(x, len) {
      rep(x, each = len / length(x), length.out = len)
    }

    DT <- data.table(N = Status12)
    DT[, day := rep_each(1:.days, .N)]
    DT[, state := rep_each(1:9, .N), by = .(day)]
    DT[, status := rep_each(-2:3, .N), by = .(day, state)]
    DT[, isolated := rep_len(0:1, .N)]
    DT[, runid := x]

    return(DT)

  }

  map_dfr(1:repeat_times, return_simulation)
}

contact_tracing_03_100 <- repeat_simulation(
  PolicyPars = set_policypars(do_contact_tracing = TRUE,
                              contact_tracing_days_before_test = 0,
                              contact_tracing_days_until_result = 3,
                              contact_tracing_success = 1))

contact_tracing_03_50 <- repeat_simulation(
  PolicyPars = set_policypars(do_contact_tracing = TRUE,
                              contact_tracing_days_before_test = 0,
                              contact_tracing_days_until_result = 3,
                              contact_tracing_success = 0.5))

contact_tracing_25_50 <- repeat_simulation(
  PolicyPars = set_policypars(do_contact_tracing = TRUE,
                              contact_tracing_days_before_test = 2,
                              contact_tracing_days_until_result = 5,
                              contact_tracing_success = 0.5))

plot_status <- function(data) {
data %>% 
  filter(status >= 1,
         state <= 6) %>% 
  group_by(state, isolated, day) %>% 
  summarise(N = sum(N)) %>% 
  ggplot(aes(day, N, colour = factor(isolated))) +
  geom_line() +
  facet_grid(vars(state), scales = "free_y") + 
  theme(legend.position = "top")
}

plot_status(contact_tracing_03_100)
#> `summarise()` regrouping by 'state', 'isolated' (override with `.groups` argument)

plot_status(contact_tracing_03_50)
#> `summarise()` regrouping by 'state', 'isolated' (override with `.groups` argument)

plot_status(contact_tracing_25_50)
#> `summarise()` regrouping by 'state', 'isolated' (override with `.groups` argument)

^{Created on 2020-05-25 by the reprex package (v0.3.0)}

[HughParsonage]

"Success" in terms of contact tracing means the number of households who are notified if one of their members is identified as infected. I'll admit I've tried to create some policy and epi pars that make this a bit more visible but I haven't been able to.

That said, I think schools are probably too eager to isolate everyone so I should probably change that.

[wfmackey]

The proportion of infected people who are isolated looks like (is almost zero) for all these scenarios though. Is that expected?

[HughParsonage]

Only 1% of asymptomatic people are tested by default.

[wfmackey]

Which is controlled by contact_tracing_only_sympto? I still don't see much of a difference:

library(tidyverse)
library(purrr)
library(data.table)
#> 
#> Attaching package: 'data.table'
#> The following objects are masked from 'package:dplyr':
#> 
#>     between, first, last
#> The following object is masked from 'package:purrr':
#> 
#>     transpose
library(covid19.model.sa2)

run_days <- 50

repeat_simulation <- function(repeat_times = 1, 
                              .days = run_days, 
                              .PolicyPars = set_policypars(),
                              .EpiPars = set_epipars(),
                              ...) {

  return_simulation <- function(x, to_return = .to_return) {

     rep_each <- function(x, len) {
        rep(x, each = len / length(x), length.out = len)
      }

    d <- 
    simulate_sa2(days_to_simulate = .days, 
                 returner = 3,
                 PolicyPars = .PolicyPars,
                 EpiPars = .EpiPars)

    # Status
      Status12 <- d[["Status12"]]

      DT <- data.table(N = Status12)
      DT[, day := rep_each(1:.days, .N)]
      DT[, state := rep_each(1:9, .N), by = .(day)]
      DT[, status := rep_each(-2:3, .N), by = .(day, state)]
      DT[, isolated := rep_len(0:1, .N)]
      DT[, runid := x]

    # New infections
      NewInfectionsByState <- d[["NewInfectionsByState"]]

      DT2 <- data.table(new_infections = NewInfectionsByState)
      DT2[, day := rep_each(1:.days, .N)]
      DT2[, state := rep_each(1:9, .N), by = .(day)]

      return(left_join(DT, DT2))

  }

  map_dfr(1:repeat_times, return_simulation)

}

plot_status <- function(data) {

data %>% 
  filter(status >= 1,
         state <= 6 & state != 5) %>% 
  group_by(state, isolated, day) %>% 
  summarise(N = sum(N)) %>% 
  ggplot(aes(day, N,
             colour = factor(isolated))) +
  geom_line() +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  scale_y_continuous(limits = c(0, NA)) + 
  facet_grid(vars(state), scales = "free_y") + 
  theme(legend.position = "top")
}

testing90 <- repeat_simulation(
  .PolicyPars = set_policypars(contact_tracing_only_sympto = .9,
                               do_contact_tracing = TRUE,
                               contact_tracing_days_before_test = 1,
                               contact_tracing_days_until_result = 3))
#> Joining, by = c("day", "state")
plot_status(testing90)
#> `summarise()` regrouping by 'state', 'isolated' (override with `.groups` argument)

testing10 <- repeat_simulation(
    .PolicyPars = set_policypars(contact_tracing_only_sympto = .1,
                               do_contact_tracing = TRUE,
                               contact_tracing_days_before_test = 1,
                               contact_tracing_days_until_result = 3))
#> Joining, by = c("day", "state")

plot_status(testing10)
#> `summarise()` regrouping by 'state', 'isolated' (override with `.groups` argument)

^{Created on 2020-05-25 by the reprex package (v0.3.0)}

[HughParsonage]

No it's not implemented yet.

https://github.com/grattan/covid19.model.sa2/blob/d3ec23d92879b0f0df1ac887f03d7d8cefc4c310/src/do_aus_simulate.cpp#L1612-L1614

[wfmackey]

Ah okay cool. The final thing I noticed was the difficulty in getting to zero active cases. Numbers seem to linger (see base and low below). What's the reason for this?

library(tidyverse)
library(purrr)
library(data.table)
#> 
#> Attaching package: 'data.table'
#> The following objects are masked from 'package:dplyr':
#> 
#>     between, first, last
#> The following object is masked from 'package:purrr':
#> 
#>     transpose
library(patchwork)
library(covid19.model.sa2)

run_days <- 50

repeat_simulation <- function(repeat_times = 1, 
                              .days = run_days, 
                              PolicyPars = set_policypars(),
                              EpiPars = set_epipars(),
                              ...) {

  return_simulation <- function(x, 
                                to_return = .to_return,
                                .PolicyPars = PolicyPars,
                                .EpiPars = EpiPars) {

     rep_each <- function(x, len) {
        rep(x, each = len / length(x), length.out = len)
     }

    d <- 
    simulate_sa2(days_to_simulate = .days, 
                 returner = 3,
                 PolicyPars = .PolicyPars,
                 EpiPars = .EpiPars)

    # Status
      Status12 <- d[["Status12"]]

      DT <- data.table(N = Status12)
      DT[, day := rep_each(1:.days, .N)]
      DT[, state := rep_each(1:9, .N), by = .(day)]
      DT[, status := rep_each(-2:3, .N), by = .(day, state)]
      DT[, isolated := rep_len(0:1, .N)]
      DT[, runid := x]

    # New infections
      NewInfectionsByState <- d[["NewInfectionsByState"]]

      DT2 <- data.table(new_infections = NewInfectionsByState)
      DT2[, day := rep_each(1:.days, .N)]
      DT2[, state := rep_each(1:9, .N), by = .(day)]

      return(left_join(DT, DT2))

  }

  map_dfr(1:repeat_times, return_simulation)

}

plot_status <- function(data) {

data %>% 
  filter(status >= 1,
         state <= 6 & state != 5) %>% 
  group_by(state, day) %>% 
  summarise(N = sum(N)) %>% 
  ggplot(aes(day, N)) +
  geom_line() +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  scale_y_continuous(limits = c(0, NA)) + 
  facet_grid(vars(state), scales = "free_y") + 
  theme(legend.position = "top")
}

plot_new_cases <- function(data) {
  data %>% 
    filter(state <= 6 & state != 5) %>% 
    group_by(state, day) %>% 
    filter(row_number() == 1) %>% 
    ggplot(aes(day, new_infections)) +
    geom_col() + 
    facet_grid(vars(state))

}

base <- repeat_simulation()
#> Joining, by = c("day", "state")
plot_new_cases(base) + plot_status(base)
#> `summarise()` regrouping by 'state' (override with `.groups` argument)

high <- repeat_simulation(
  EpiPars = set_epipars(q_workplace = .3,
                        a_workplace_rate = 1,
                        q_supermarket = .3),
  PolicyPars = set_policypars(supermarkets_open = TRUE,
                              workplaces_open = 1, 
                              workplace_size_max = 1000))
#> Joining, by = c("day", "state")

plot_new_cases(high) + plot_status(high)
#> `summarise()` regrouping by 'state' (override with `.groups` argument)

low <- repeat_simulation(
  EpiPars = set_epipars(q_workplace = .01,
                        a_workplace_rate = .1,
                        q_supermarket = .001),
  PolicyPars = set_policypars(supermarkets_open = TRUE,
                              workplaces_open = .2, 
                              workplace_size_max = 10))
#> Joining, by = c("day", "state")

plot_new_cases(low) + plot_status(low)
#> `summarise()` regrouping by 'state' (override with `.groups` argument)

^{Created on 2020-05-25 by the reprex package (v0.3.0)}

[HughParsonage]

(This is excellent feedback Will. Thank you!)

Bear in mind that 'most people' visit supermarkets and cafes at least once in during an incubation period, so once persons are distributed throughout SA2s you get that effect that sustains the epidemic. Setting supermarkets_open or cafes_open to FALSE generally causes the model to reach zero within 50 days.

[wfmackey]

In the case of low, there are no new cases in SA for ~30 days and there is still an active case; and the same in Tas for 50 days. Is this because of absurdly long outliers in incubation/illness periods?

[wfmackey]

Even with everything closed so there are no new cases anywhere (according to NewInfectionsByState), active cases still linger permanently. 100 days in the example below:

lower <- repeat_simulation(.days = 100,
  EpiPars = set_epipars(q_workplace = .01,
                        a_workplace_rate = .1,
                        q_supermarket = .001),
  PolicyPars = set_policypars(supermarkets_open = FALSE,
                              cafes_open = FALSE,
                              workplaces_open = .2, 
                              workplace_size_max = 10))

plot_new_cases(lower) + plot_status(lower)

[wfmackey]

Breaking them out, it is both symptomatic and asymptomatic cases that linger:

[HughParsonage]

Hmm I can't reproduce that lingering effect. Mine all hit zero by day 77 at the latest.

[wfmackey]

Weird! Can you share your run of the following?

library(tidyverse)
  options(dplyr.summarise.inform = FALSE)
library(purrr)
library(data.table)
#> 
#> Attaching package: 'data.table'
#> The following objects are masked from 'package:dplyr':
#> 
#>     between, first, last
#> The following object is masked from 'package:purrr':
#> 
#>     transpose
library(patchwork)
library(covid19.model.sa2)

run_days <- 50

repeat_simulation <- function(repeat_times = 1, 
                              .days = run_days, 
                              PolicyPars = set_policypars(),
                              EpiPars = set_epipars(),
                              ...) {

  return_simulation <- function(x, 
                                to_return = .to_return,
                                .PolicyPars = PolicyPars,
                                .EpiPars = EpiPars) {

     rep_each <- function(x, len) {
        rep(x, each = len / length(x), length.out = len)
     }

    d <- 
    simulate_sa2(days_to_simulate = .days, 
                 returner = 3,
                 PolicyPars = .PolicyPars,
                 EpiPars = .EpiPars)

    # Status
      Status12 <- d[["Status12"]]

      DT <- data.table(N = Status12)
      DT[, day := rep_each(1:.days, .N)]
      DT[, state := rep_each(1:9, .N), by = .(day)]
      DT[, status := rep_each(-2:3, .N), by = .(day, state)]
      DT[, isolated := rep_len(0:1, .N)]
      DT[, runid := x]

    # New infections
      NewInfectionsByState <- d[["NewInfectionsByState"]]

      DT2 <- data.table(new_infections = NewInfectionsByState)
      DT2[, day := rep_each(1:.days, .N)]
      DT2[, state := rep_each(1:9, .N), by = .(day)]

      return(left_join(DT, DT2))

  }

  map_dfr(1:repeat_times, return_simulation)

}

plot_status <- function(data) {

data %>% 
  filter(status >= 1,
         state <= 6 & state != 5) %>% 
  group_by(status, state, day) %>% 
  summarise(N = sum(N)) %>% 
  ggplot(aes(day, N, colour = factor(status))) +
  geom_line() +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  scale_y_continuous(limits = c(0, NA)) + 
  facet_grid(vars(state), scales = "free_y") + 
  theme(legend.position = "top")
}

plot_new_cases <- function(data) {
  data %>% 
    filter(state <= 6 & state != 5) %>% 
    group_by(state, day) %>% 
    filter(row_number() == 1) %>% 
    ggplot(aes(day, new_infections)) +
    geom_col() + 
    facet_grid(vars(state))

}

lower <- repeat_simulation(.days = 100,
  EpiPars = set_epipars(q_workplace = .01,
                        a_workplace_rate = .1,
                        q_supermarket = .001),
  PolicyPars = set_policypars(supermarkets_open = FALSE,
                              cafes_open = FALSE,
                              workplaces_open = .2, 
                              workplace_size_max = 10))
#> Joining, by = c("day", "state")

plot_status(lower)

^{Created on 2020-05-25 by the reprex package (v0.3.0)}

[HughParsonage]

Ah I see now. I was looking at non-isolated cases. Overcorrected the previous issue.

[wfmackey]

Cool cool. Once that's fixed, I'll have another look and then I think we're ready to run a PolicyGrid through the servers again. What do you think?

[HughParsonage]

Yep sounds great. Thanks so much! I've pushed that change.

[wfmackey]

Looking good! The default epi and policy pars now send cases increasing rapidly. That's not necessarily an issue -- but is it what you would have expected?

library(tidyverse)
  options(dplyr.summarise.inform = FALSE)
library(purrr)
library(data.table)
#> 
#> Attaching package: 'data.table'
#> The following objects are masked from 'package:dplyr':
#> 
#>     between, first, last
#> The following object is masked from 'package:purrr':
#> 
#>     transpose
library(patchwork)
library(covid19.model.sa2)

run_days <- 50

repeat_simulation <- function(repeat_times = 1, 
                              .days = run_days, 
                              PolicyPars = set_policypars(),
                              EpiPars = set_epipars(),
                              ...) {

  return_simulation <- function(x, 
                                to_return = .to_return,
                                .PolicyPars = PolicyPars,
                                .EpiPars = EpiPars) {

     rep_each <- function(x, len) {
        rep(x, each = len / length(x), length.out = len)
     }

    d <- 
    simulate_sa2(days_to_simulate = .days, 
                 returner = 3,
                 PolicyPars = .PolicyPars,
                 EpiPars = .EpiPars)

    # Status
      Status12 <- d[["Status12"]]

      DT <- data.table(N = Status12)
      DT[, day := rep_each(1:.days, .N)]
      DT[, state := rep_each(1:9, .N), by = .(day)]
      DT[, status := rep_each(-2:3, .N), by = .(day, state)]
      DT[, isolated := rep_len(0:1, .N)]
      DT[, runid := x]

    # New infections
      NewInfectionsByState <- d[["NewInfectionsByState"]]

      DT2 <- data.table(new_infections = NewInfectionsByState)
      DT2[, day := rep_each(1:.days, .N)]
      DT2[, state := rep_each(1:9, .N), by = .(day)]

      return(left_join(DT, DT2))

  }

  map_dfr(1:repeat_times, return_simulation)

}

plot_status <- function(data) {

data %>% 
  filter(status >= 1,
         state <= 6 & state != 5) %>% 
  group_by(status, state, day) %>% 
  summarise(N = sum(N)) %>% 
  ggplot(aes(day, N, colour = factor(status))) +
  geom_line() +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  scale_y_continuous(limits = c(0, NA)) + 
  facet_wrap(vars(state), scales = "free_y", ncol = 1) + 
  theme(legend.position = "top")
}

plot_new_cases <- function(data) {
  data %>% 
    filter(state <= 6 & state != 5) %>% 
    group_by(state, day) %>% 
    filter(row_number() == 1) %>% 
    ggplot(aes(day, new_infections)) +
    geom_col() + 
    facet_wrap(vars(state), scales = "free_y", ncol = 1)

}

base <- repeat_simulation()
#> Joining, by = c("day", "state")
plot_new_cases(base) + plot_status(base)

high <- repeat_simulation(
  EpiPars = set_epipars(q_workplace = .3,
                        a_workplace_rate = 1,
                        q_supermarket = .3),
  PolicyPars = set_policypars(supermarkets_open = TRUE,
                              cafes_open = TRUE,
                              workplaces_open = 1, 
                              workplace_size_max = 1000))
#> Joining, by = c("day", "state")

plot_new_cases(high) + plot_status(high)

low <- repeat_simulation(.days = 100,
  EpiPars = set_epipars(q_workplace = .01,
                        a_workplace_rate = .1,
                        q_supermarket = .001),
  PolicyPars = set_policypars(supermarkets_open = FALSE,
                              cafes_open = FALSE,
                              workplaces_open = .2, 
                              workplace_size_max = 10))
#> Joining, by = c("day", "state")

plot_new_cases(low) + plot_status(low)

short_incubation <- repeat_simulation(
  EpiPars = set_epipars(incubation_distribution = "pois",
                        incubation_mean = 4,
                        incubation_sigma = .3))
#> Joining, by = c("day", "state")

long_incubation <- repeat_simulation(
  EpiPars = set_epipars(incubation_distribution = "pois",
                        incubation_mean = 10,
                        incubation_sigma = 1))
#> Joining, by = c("day", "state")

plot_status(short_incubation) + plot_status(long_incubation)

^{Created on 2020-05-25 by the reprex package (v0.3.0)}

[wfmackey]

If you're happy with how that all looks, I'm happy to start the run again. Same grid as covid19-results-0522?

[wfmackey]

Or are there different elements to the ones we had in covid19-results-0522 that you think would be interesting for the user to see, or are particularly policy relevant?

[HughParsonage]

I think the sensitivity of some of the contact probabilities (q_') would be interesting to view, but I think we should focus on policy parameters. To that end let's create a CJ but CJ( ..., policy ) since this will cover the policy grids first then if they don't take too long we can test out the same policy grid with different epi pars.

There may be some calibration that is needed too.

[wfmackey]

That's a good idea. But Stephen is v keen to get a version of the shiny app to people tomorrow, so can we first run through the original covid19-results-0522 grid (or something similar), and I can get the app up-and-running with that tomorrow morning.

Then we can keep adding to it?

[HughParsonage]

Sure.

[wfmackey]

Thanks! Let me know when the new runs start feeding through so I can check it again tonight.

[HughParsonage]

Wait which grid did you mean?

Grid <-
  CJ(schools_open = c(TRUE, FALSE),
     only_Year12 = c(TRUE, FALSE),
     school_days_per_wk = c(1L, 3L, 5L),
     do_contact_tracing = c(TRUE, FALSE),
     contact_tracing_days_before_test = c(1:2),
     contact_tracing_days_until_result = c(3L, 5L, 7L),
     workplaces_open = c(0, 0.5, 1.0),
     workplace_size_max = c(10L, 75L),
     incubation_distribution = c("pois", "cauchy"),
     incubation_mean = c(4L, 7L, 10L),
     a_workplace_rate = c(0.50, 0.75),
     q_school = c(1/3000, 0.01, 0.05),
     q_household = c(0.04, 0.09, 0.15)) %>%
  .[only_Year12 %implies% schools_open] %>%
  .[(school_days_per_wk > 1L) %implies% schools_open] %>%
  .[(workplace_size_max > 10) %implies% (workplaces_open > 0)]

[wfmackey]

Yep