Create functions for generating forecast distribution

index.qmd

@@ -131,29 +131,29 @@
   contact_rate_weekend <- params[7]
 
   # Set recovery rate
   set_param(ef_model, "Recovery rate", recovery_rate)
 
   # Global event to change contact and transmission rates
   change_c_and_t_rates <- function(model) {
     # Get the current model day (step)
     current_model_day <- today(model)
 
     ## Update contact rate based on weekday/weekend
     if (any(c(6, 0) %in% (current_model_day %% 7L))) {
       set_param(model, "Contact rate", contact_rate_weekend)
     } else {
       set_param(model, "Contact rate", contact_rate_weekday)
     }
 
     ## Update transmission rate each season
     if (current_model_day == spring_start) {
       set_param(model, "Transmission rate", t_rate_spring)
     } else if (current_model_day == summer_start) {
       set_param(model, "Transmission rate", t_rate_summer)
     } else if (current_model_day == fall_start) {
       set_param(model, "Transmission rate", t_rate_fall)
     } else if (current_model_day == winter_start) {
       set_param(model, "Transmission rate", t_rate_winter)
     }
 
     invisible(model)
@@ -161,8 +161,8 @@
 
   # Add global event to the model
   change_c_and_t_event_name <- "Change Contact and Transmission Rates"
   globalevent_fun(change_c_and_t_rates, name = change_c_and_t_event_name) |>
     add_globalevent(model = ef_model)
 
   # Run the model
   verbose_off(ef_model)
@@ -173,10 +173,12 @@
 
   # Get infected cases
   hist_matrix <- get_hist_transition_matrix(ef_model)
-
-  # - Model returns 91 values so we drop the last value
+  # - Drop the last day of data because the model returns data for (model_ndays + 1) days
   infected_cases <- head(
-    hist_matrix[grep("Infected", hist_matrix$state_to), "counts"],
+    hist_matrix[
+      hist_matrix$state_to == "Infected" &
+        hist_matrix$state_from == "Susceptible",
+      c("counts")],
     -1
   )
 
@@ -334,16 +336,21 @@
 
 We can now run the forecast.
 Our model prevalence is set according to the reported case counts of the most recent day of the UDHHS data.
-For the other model parameters, we find the second quantile from the LFMCMC accepted parameters.
-We then run the SIR connected model with the new prevalence and the selected parameters.
-The model runs for two weeks, giving us a 14-day forecast of COVID-19 in Utah.
+We then take a sample of n = 50 from the LFMCMC accepted parameters and run the SIR connected model with the new prevalence for each set of parameters.
+Each simulation is for two weeks, giving us a 14-day forecast of COVID-19 in Utah.
+The forecast mean is shown below along with the 50% and 95% confidence intervals.
+
 ```{r}
+#| label: run-forecast
+library(tidyr)
+forecast_sample_n <- 50
+# Run model for 2 weeks
+model_ndays <- 14
 # Compute prevalance based on most recent day
 forecast_prevalence <- forecast_data_cases[90] / model_n
 # Select parameters from second quantile of accepted params
 accepted_params <- get_accepted_params(lfmcmc_model)
-params_median <- accepted_params[-c(1:(n_samples / 2)), ] |>
-  apply(2, quantile, probs = c(.5))
+params_sample <- accepted_params[sample(nrow(accepted_params), forecast_sample_n), ]
 # Create the new model
 ef_model <- ModelSIRCONN(
   name              = "COVID-19",
@@ -353,14 +360,71 @@
   transmission_rate = init_transmission_rate,
   recovery_rate     = init_lfmcmc_params[1]
 )
-# Run model for 2 weeks
-model_ndays <- 14
 
-forecasted_cases <- lfmcmc_simulation_fun(params_median)
+# Run the simulation for each set of params in the sample
+forecast_dist <- apply(params_sample, 1, lfmcmc_simulation_fun)
+# Compute mean and median
+fc_sample_mean <- apply(forecast_dist, 1, mean)
+fc_sample_median <- apply(forecast_dist, 1, median)
+
+# Compute confidence interval
+# - Standard error
+std_error <- function(x) sd(x) / sqrt(length(x))
+margin_error <- function(ci, se) {
+  qt(ci, df = model_ndays - 1) * se
+}
+fc_sample_se <- apply(forecast_dist, 1, std_error)
+# - Margin of error
+fc_sample_me95 <- margin_error(0.975, fc_sample_se)
+fc_sample_me50 <- margin_error(0.75, fc_sample_se)
+
+# Combine observed data with sample median for plotting forecast
+observed_df <- data.frame(
+  date = forecast_data$Date[60:90],
+  counts = forecast_data$Daily.Cases[60:90],
+  observed = TRUE,
+  lb_95 = NA,
+  ub_95 = NA,
+  lb_50 = NA,
+  ub_50 = NA
+)
+sample_df <- data.frame(
+  date = forecast_data$Date[90] + 0:13,
+  counts = fc_sample_mean,
+  observed = FALSE,
+  lb_95 = fc_sample_mean - fc_sample_me95,
+  ub_95 = fc_sample_mean + fc_sample_me95,
+  lb_50 = fc_sample_mean - fc_sample_me50,
+  ub_50 = fc_sample_mean + fc_sample_me50
+)
 
-# Print Model Summary and Plot Incidence
-summary(ef_model)
-plot_incidence(ef_model)
+forecast_df <- rbind(observed_df, sample_df)
+```
+```{r}
+#| label: plot-forecast
+# Use color-blind-friendly palette:
+cbb_light_blue <- "#56B4E9"
+cbb_palette <- c(cbb_light_blue, "black")
+
+ggplot(forecast_df,
+  aes(x = date)) +
+  geom_ribbon(aes(ymin = lb_95, ymax = ub_95),
+    fill = cbb_light_blue,
+    alpha = 0.4,
+    na.rm = TRUE) +
+  geom_ribbon(aes(ymin = lb_50, ymax = ub_50),
+    fill = cbb_light_blue,
+    alpha = 0.4,
+    na.rm = TRUE) +
+  geom_point(aes(y = counts,
+    color = observed)) +
+  geom_line(aes(y = counts,
+    color = observed)) +
+  labs(x = "Date", y = "Daily Cases") +
+  scale_colour_manual(values = cbb_palette,
+    labels = c("Forecasted Cases", "Observed Cases")) +
+  scale_y_continuous(n.breaks = 20) +
+  theme_bw()
 ```
 
 ## Acknowledgements