Enable basic model calibration

.devcontainer/Dockerfile

@@ -3,7 +3,7 @@
 FROM rocker/tidyverse:4.4.1
 
 # Install epiworldR
-RUN install2.r epiworldR
+RUN installGithub.r UofUEpiBio/epiworldR
 
 # Run Bash
 CMD ["bash"]

.gitignore

@@ -47,3 +47,4 @@ Temporary Items
 .Rhistory
 .RData
 .Ruserdata
+*.Rproj

README.md

@@ -3,21 +3,20 @@
 This project is currently under construction and will have the following functionality:
 - Gets health data from public source (e.g., DHHS) for Salt Lake county or all of Utah
 - Calibrates a model using the data
-- Does a forecast using the calibrated model and the [`epiworld`](https://github.com/UofUEpiBio/epiworld/) simulation tool
+- Does a forecast using the calibrated model and the [`epiworldR`](https://github.com/UofUEpiBio/epiworldR/) simulation tool
 - Generates a report with figures and descriptions
 - Publishes the report via GitHub Pages
 
 The forecast will update automatically each week using GitHub Actions.
 
 ## Data Sources
-We are currently planning to use the [weekly reports](https://coronavirus.utah.gov/case-counts/) published by Utah DHHS on COVID-19.
+We are currently using the [weekly reports](https://coronavirus.utah.gov/case-counts/) published by Utah DHHS on COVID-19.
 
 Other possible data sources include:
 - Utah DHHS publishes [weekly reports](https://coronavirus.utah.gov/case-counts/) on COVID-19, Influenza, and Respiratory syncytial virus (RSV).
 - Germ Watch could be another great public source, but is currently unavailable (links, such as those mentioned [here](https://epi.utah.gov/influenza-reports/) redirect to Intermountain's home page)
 - DELPHI maintains a frequently updated COVID data API [here](https://cmu-delphi.github.io/delphi-epidata/api/covidcast.html) and additional endpoints (less frequently updated) for influenza, dengue, and norovirus [here](https://cmu-delphi.github.io/delphi-epidata/api/README.html)
 - CDC's [public datasets](https://data.cdc.gov), some are updated infrequently, others are weekly estimates (e.g., [weekly flu vaccine estimates](https://data.cdc.gov/Vaccinations/Weekly-Cumulative-Estimated-Number-of-Influenza-Va/ysd3-txwj/about_data)
-- Germ watch
 
 ### Other Relevant Resources
 It's worth also taking a look at:
@@ -27,4 +26,4 @@ It's worth also taking a look at:
 
 ## Code of Conduct
 
-Please note that the epiworld-forecasts project is released with a [Contributor Code of Conduct](./CODE_OF_CONDUCT.md). By contributing to this project, you agree to abide by its terms. More information about how to contribute to the project can be found under [`DEVELOPMENT.md`](DEVELOPMENT.md).
+Please note that the Epiworld Forecasts project is released with a [Contributor Code of Conduct](./CODE_OF_CONDUCT.md). By contributing to this project, you agree to abide by its terms. More information about how to contribute to the project can be found under [`DEVELOPMENT.md`](DEVELOPMENT.md).

index.qmd

@@ -22,6 +22,7 @@
 data_url <- "https://coronavirus-dashboard.utah.gov/Utah_COVID19_data.zip"
 target_file_regex <- "Trends_Epidemic+"
 forecast_data <- get_forecast_data(data_url, target_file_regex)
+forecast_data$Date <- as.Date(forecast_data$Date)
 # Check for errors
 if (length(forecast_data) > 1) {
   # Plot the observed data
@@ -35,22 +36,225 @@
 }
 ```
 
-## Epiworld Forecast
-We calibrate a SIR Connected model using the above data and run the model in epiworldR.
-Here are the results of a single model run:
+Extracting only the past 90 days:
+
+```{r}
+library(ggplot2)
+
+# Getting the last 3 months of data
+last_date <- max(forecast_data$Date)
+forecast_data <- forecast_data[forecast_data$Date > (last_date - 90), ]
+ggplot(forecast_data, aes(x = Date, y = Daily.Cases)) +
+  geom_line(aes(group = 1)) +
+  labs(x = "Date", y = "Daily Cases", title = "Daily COVID-19 Cases in Utah")
+
+# Selecting start date of each season from data
+get_date_season <- function(date) {
+  date_month <- as.integer(format(as.Date(date, format = "%d/%m/%Y"), "%m"))
+
+  if (date_month >= 3 && date_month <= 5) {
+    return("spring")
+  } else if (date_month >= 6 && date_month <= 8) {
+    return("summer")
+  } else if (date_month >= 9 && date_month <= 11) {
+    return("fall")
+  } else {
+    return("winter")
+  }
+}
+forecast_seasons <- mapply(get_date_season, forecast_data$Date)
+
+spring_start <- match("spring", forecast_seasons, nomatch = -1)
+summer_start <- match("summer", forecast_seasons, nomatch = -1)
+fall_start <- match("fall", forecast_seasons, nomatch = -1)
+winter_start <- match("winter", forecast_seasons, nomatch = -1)
+```
+
+Creating an LFMCMC model that adjusts contact and transmission rates based on the date:
+
 ```{r}
 library(epiworldR)
-model_sir <- ModelSIRCONN(
+
+# Set parameters
+model_seed <- 112
+model_ndays <- 90
+
+init_lfmcmc_params <- c(
+  1 / 7,  # recovery_rate
+  0.05,   # t_rate_spring
+  0.04,   # t_rate_summer
+  0.06,   # t_rate_fall
+  0.07,   # t_rate_winter
+  10,     # contact_rate_weekday
+  2       # contact_rate_weekend
+)
+
+# Create the base SIRCONN model and run once to ensure it is fully initialized
+ef_model <- ModelSIRCONN(
   name              = "COVID-19",
-  n                 = 50000,
+  n                 = 100000,
   prevalence        = 0.0001,
-  contact_rate      = 2,
-  transmission_rate = 0.5,
-  recovery_rate     = 1 / 3
+  contact_rate      = 10,
+  transmission_rate = 0.05,
+  recovery_rate     = init_lfmcmc_params[1]
+)
+verbose_off(ef_model)
+run(ef_model, ndays = model_ndays, seed = model_seed)
+```
+
+```{r}
+# Define the LFMCMC simulation function
+lfmcmc_simulation_fun <- function(params) {
+  # Extract parameters
+  recovery_rate <- params[1]
+  t_rate_spring <- params[2]
+  t_rate_summer <- params[3]
+  t_rate_fall <- params[4]
+  t_rate_winter <- params[5]
+  contact_rate_weekday <- params[6]
+  contact_rate_weekend <- params[7]
+  # contact_rate_holiday <- params[8]
+
+  # Set recovery rate
+  set_param(ef_model, "Recovery rate", recovery_rate)
+
+  # Add a global event that changes the contact rate parameter
+  change_c_and_t_rates <- function(model) {
+    current_model_day <- today(model)
+
+    ## Update contact rate
+    # Check whether today is a weekday or 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)
+    }
+
+    # TODO: Add cases to update contact rate for specific holidays (extract from the 90 days of data)
+
+    ## Update transmission rate each season
+    # TODO: do we need current_model_day + 1?
+    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)
+  }
+
+  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)
+  run(ef_model, ndays = model_ndays, seed = model_seed)
+
+  # Clean up the global event so we can set it the next time we run the simulation
+  rm_globalevent(ef_model, change_c_and_t_event_name)
+
+  # Printing Model Summary (TODO: delete this)
+  # summary(ef_model)
+  # plot_incidence(ef_model)
+
+  # Return infected cases
+  hist_matrix <- get_hist_transition_matrix(ef_model)
+  # - model returns 91 values, drop last value
+  infected_cases <- head(hist_matrix[grep("Infected", hist_matrix$state_to), "counts"], -1)
+  return(infected_cases)
+}
+
+# lfmcmc_simulation_fun(init_lfmcmc_params)
+
+# Expects dat to include case counts
+lfmcmc_summary_fun <- function(dat) {
+  # Extract summary statistics from the data
+  time_to_peak <- which.max(dat)
+  # TODO: compute as % of pop
+  size_of_peak <- dat[time_to_peak]
+  sd_cases <- sd(dat)
+  mean_cases <- mean(dat)
+  # TODO: Repeat for each season in data
+  return(as.integer(c(
+    time_to_peak,
+    size_of_peak,
+    sd_cases,
+    mean_cases
+  )))
+}
+
+lfmcmc_proposal_fun <- function(params_prev) {
+  res_1 <- abs(params_prev[1] + rnorm(1, mean = 0, sd = 0.5))
+  res_2_to_5 <- abs(params_prev[2:5] + rnorm(4, mean = 0, sd = 0.1))
+  res_6_to_7 <- abs(params_prev[6:7] + rnorm(2, mean = 0, sd = 1))
+
+  res <- c(
+    res_1,
+    res_2_to_5,
+    res_6_to_7
+  )
+
+  # res <- params_prev + rnorm(length(params_prev), )
+  return(res)
+}
+
+# lfmcmc_proposal_fun(init_lfmcmc_params)
+
+lfmcmc_kernel_fun <- function(stats_now, stats_obs, epsilon) {
+
+  ans <- sum(mapply(function(v1, v2) (v1 - v2)^2,
+    stats_obs,
+    stats_now))
+
+  return(ifelse(sqrt(ans) < epsilon, 1.0, 0.0))
+}
+
+# Create the LFMCMC model
+lfmcmc_model <- LFMCMC(ef_model) |>
+  set_simulation_fun(lfmcmc_simulation_fun) |>
+  set_summary_fun(lfmcmc_summary_fun) |>
+  set_proposal_fun(lfmcmc_proposal_fun) |>
+  set_kernel_fun(lfmcmc_kernel_fun) |>
+  set_observed_data(forecast_data$Daily.Cases)
+```
+
+```{r}
+# Run the LFMCMC simulation
+run_lfmcmc(
+  lfmcmc = lfmcmc_model,
+  params_init_ = init_lfmcmc_params,
+  n_samples_ = 100,
+  epsilon_ = 1.0,
+  seed = model_seed
+)
+
+# Print the results
+par_names <- c(
+  "recovery_rate",
+  "t_rate_spring",
+  "t_rate_summer",
+  "t_rate_fall",
+  "t_rate_winter",
+  "contact_rate_weekday",
+  "contact_rate_weekend"
 )
+set_par_names(lfmcmc_model, par_names)
+
+print(lfmcmc_model)
+# TODO: Process results as percentage of population, relative to the population size in the model
+```
 
+## Epiworld Forecast
+We calibrate a SIR Connected model using the above data and run the model in epiworldR.
+Here are the results of a single model run:
+```{r}
 # Printing Model Summary
-summary(model_sir)
+# run(ef_model, ndays = model_ndays, seed = model_seed)
 ```
 
 ## Methodology