From 96ab826326e5ccf4a5d3e0934132afcf2148b948 Mon Sep 17 00:00:00 2001
From: cazodi <cazodi@svi.edu.au>
Date: Fri, 3 Apr 2020 14:54:15 +1100
Subject: [PATCH] add code for setting up time sensitive experiments where
parameter values change over time
---
R/results-parse.R | 12 ++++---
R/time_experiments.R | 68 +++++++++++++++++++++++++++++++++++++
vignettes/sirplus_intro.Rmd | 45 +++++++++++++++---------
3 files changed, 105 insertions(+), 20 deletions(-)
create mode 100644 R/time_experiments.R
diff --git a/R/results-parse.R b/R/results-parse.R
index 843e6c9..80ac246 100644
--- a/R/results-parse.R
+++ b/R/results-parse.R
@@ -153,7 +153,8 @@ plot_models <- function(sims = baseline_sim,
scale_x_date(date_breaks = "1 week", date_labels = "%m-%d") +
scale_colour_manual(values = compcols, labels = complabels) +
labs(title = plot_title, x = x_axis, y = "Prevalence") +
- theme_bw()
+ theme_bw() +
+ theme(axis.text.x = element_text(angle = 90))
} else (
plot_df %>% ggplot(aes(x = Date, y = count, colour = compartment,
linetype = sim)) +
@@ -162,7 +163,8 @@ plot_models <- function(sims = baseline_sim,
geom_line(size = 1.5, alpha = 0.8) +
scale_colour_manual(values = compcols, labels = complabels) +
labs(title = plot_title, x = x_axis, y = "Prevalence") +
- theme_bw()
+ theme_bw() +
+ theme(axis.text.x = element_text(angle = 90))
)
}else{
if(length(sim_id) == 1){
@@ -173,7 +175,8 @@ plot_models <- function(sims = baseline_sim,
geom_line(size = 1.5, alpha = 0.8) +
scale_colour_manual(values = compcols, labels = complabels) +
labs(title = plot_title, x = x_axis, y = "Prevalence") +
- theme_bw()
+ theme_bw() +
+ theme(axis.text.x = element_text(angle = 90))
} else (
plot_df %>% ggplot(aes(x = Date, y = count+1, colour = compartment,
linetype = sim)) +
@@ -183,7 +186,8 @@ plot_models <- function(sims = baseline_sim,
geom_line(size = 1.5, alpha = 0.8) +
scale_colour_manual(values = compcols, labels = complabels) +
labs(title = plot_title, x = x_axis, y = "Prevalence") +
- theme_bw()
+ theme_bw() +
+ theme(axis.text.x = element_text(angle = 90))
)
}
}
diff --git a/R/time_experiments.R b/R/time_experiments.R
new file mode 100644
index 0000000..5b2afac
--- /dev/null
+++ b/R/time_experiments.R
@@ -0,0 +1,68 @@
+#' Generating time-dependent parameters
+#'
+#' Function to generate parameter values for the time range of the simulation
+#' that can change over time.
+#'
+#' @param nstep Number of time steps to generate parameter values for.
+#' @param vals List of parameter values to include over the nsteps.
+#' @param timing List of the step numbers at which to start changes, with the
+#' last number reflecting when to hit the last parameter value in vals.
+#'
+#' @return list of parameter values for length t
+#'
+vary_param <- function(nstep = nstep, vals = vals, timing = timing) {
+
+ stopifnot(length(vals) == length(timing))
+ y <- list()
+
+ for(t in seq(1:nstep)){
+ if(t <= timing[1]){ # If before first jump set to val[1]
+ y.t <- vals[1]
+ }else if(t > tail(timing, n=1)){ # If after last jump set to val[-1]
+ y.t <- tail(vals, n=1)
+ }else{ # If intermediate step, calculate...
+ for(j in (2:length(timing))){
+ if(t > timing[j - 1] & t <= timing[j]){
+ start <- vals[j - 1]
+ end <- vals[j]
+ start_t <- timing[j - 1]
+ end_t <- timing[j]
+ y.t <- start - (t-start_t)*(start - end)/(end_t - start_t)
+ }
+ }
+ }
+ y <- append(y, y.t)
+ }
+
+ return(unlist(y))
+
+}
+
+
+vary_fun <- function(t,start_t=1,start_val,
+ end_t = 9,
+ end_val,
+ end_t2=22, end_val2,
+ do_relax = FALSE,
+ relax_start = NULL,
+ relax_end = NULL,
+ relax_val = NULL) {
+ stopifnot(start_t < end_t)
+
+ if(do_relax){
+ stopifnot(relax_start <= relax_end)
+ stopifnot(end_t < relax_start)
+
+ y <- ifelse(t <= start_t, start_val,
+ ifelse(t <= end_t, start_val - (t-start_t)*(start_val - end_val)/(end_t - start_t),
+ ifelse(t <= end_t2,end_val - (t - end_t)*(end_val - end_val2)/(end_t2 - end_t),
+ ifelse(t<=relax_start, end_val2,
+ ifelse(t <= relax_end, end_val2 - (t-relax_start)*(end_val2 - relax_val)/(relax_end - relax_start), relax_val)))))
+ } else {
+ y <- ifelse(t <= start_t, start_val,
+ ifelse(t <= end_t, start_val - (t-start_t)*(start_val - end_val)/(end_t - start_t),
+ ifelse(t <= end_t2,end_val - (t - end_t)*(end_val - end_val2)/(end_t2 - end_t),end_val2)))
+ }
+ return(round(y,4))
+
+}
diff --git a/vignettes/sirplus_intro.Rmd b/vignettes/sirplus_intro.Rmd
index 150233f..6dec74f 100644
--- a/vignettes/sirplus_intro.Rmd
+++ b/vignettes/sirplus_intro.Rmd
@@ -56,11 +56,12 @@ default parameters for disease spread (i.e. no additional interventions).
```{r simulate baselines}
s.num <- 2000 # number susceptible
-i.num <- 50 # number infected
-q.num <- 10 # number in self-isolation
+i.num <- 15 # number infected
+q.num <- 5 # number in self-isolation
h.num <- 1 # number in the hospital
+nstep <- 90 # number of steps (e.g. days) to simulate
-baseline_sim <- simulate_seiqhrf(s.num = s.num, i.num = i.num,
+baseline_sim <- simulate_seiqhrf(nstep = nstep, s.num = s.num, i.num = i.num,
q.num = q.num, h.num = h.num)
head(baseline_sim$df)
@@ -86,7 +87,6 @@ over time in each compartment.
```{r viz prevalence}
plot_models(sims = baseline_sim,
- time_lim = 50,
start_date = lubridate::ymd("2020-01-01"),
comp_remove = c('s.num', 'r.num'),
plot_title = 'Baseline Model')
@@ -94,22 +94,35 @@ plot_models(sims = baseline_sim,
## Run an experiment
-With the sirplus package you can also set up experiments. For example, let's say that one week after the beginning of the epidemic, schools and non-essential businesses are closed to encourage social distancing. We can model the impacts that policies may have and compare to our baseline model. In this example, we model a policy-based increase in social distancing by ramping down the act.rate.e (# of exposure events, or acts, between individuals in the E and S compartment, per day) after 7 days and plot this model next to our baseline.
+With the sirplus package you can also set up experiments. We will set up two
+experiments here:
-```{r experiment 1}
-closures_RampOnday7 <- function(t) {
- ifelse(t < 7, 10, ifelse(t <= 14, 10 - (t-7)*(10 - 5)/7, 5))
-}
+- Experiment #1: One week after the beginning of the epidemic, schools and non-essential businesses are closed to encourage social distancing. This causes the act.rate to gradually drop from 10 to 6 over the course of the next week. In this experiment, we imagine these policies are never lifted, so act.rate remains at 6 for the duration of the simulation.
+- Experiment #2: Again, one week after the beginning of the epidemic, social distancing policies are put into place resulting in act.rate dropping from 10 to 6 over the next week. But after two weeks these policies are lifted and the act.rate returns to normal within the next week.
+```{r Experiment example with act.rate}
+# Experiment #1
+vals <- c(10, 7)
+timing <- c(7, 14)
+act_rate <- vary_param(nstep = nstep, vals = vals, timing = timing)
-closures_sim <- simulate_seiqhrf(s.num = s.num, i.num = i.num,
- q.num = q.num, h.num = h.num,
- act.rate.e = closures_RampOnday7(1:366),
- act.rate.i = closures_RampOnday7(1:366))
+sim_exp <- simulate_seiqhrf(nstep = nstep, s.num = s.num, i.num = i.num,
+ q.num = q.num, h.num = h.num, act.rate.e = act_rate,
+ act.rate.i = act_rate * 0.5)
-plot_models(sims = c(baseline_sim, closures_sim),
- sim_id = c('Baseline', 'Closures (d14)'),
- time_lim = 50,
+# Experiment #2
+vals <- c(10, 7, 7, 10)
+timing <- c(7, 14, 21, 28)
+act_rate_relax <- vary_param(nstep, vals = vals, timing = timing)
+
+sim_exp_relax <- simulate_seiqhrf(nstep = nstep, s.num = s.num, i.num = i.num,
+ q.num = q.num, h.num = h.num,
+ act.rate.e = act_rate_relax,
+ act.rate.i = 0.5 * act_rate_relax)
+
+# Compare experiments 1 and 2 to the baseline simulation
+plot_models(sims = c(baseline_sim, sim_exp, sim_exp_relax),
+ sim_id = c('Baseline', 'Closures', 'Closures (2 mo)'),
start_date = lubridate::ymd("2020-01-01"),
comp_remove = c('s.num', 'r.num'),
plot_title = 'Closures Experiment')
--
GitLab