add code for setting up time sensitive experiments where parameter values change over time

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))
         )
     }
 }

R/time_experiments.R

+#' 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))
+
+}

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')