README.Rmd

---
output:
  md_document:
    variant: markdown_github
---

<!-- README.md is generated from README.Rmd. Please edit that file -->

### Quantifying the impact of COVID-19 control measures using a Bayesian model of physical distancing

This repository contains code associated with a manuscript investigating the impact of COVID-19 control measures in British Columbia, Canada.

The main statistical model written in [Stan](https://mc-stan.org/) is available [here](analysis/seeiqr.stan) and the main R function that calls this model for a vector of daily case counts is available [here](analysis/fit_seeiqr.R). A function to make projection plots is available [here](analysis/make_projection_plot.R). This model may be released at a later date in a proper R package form.

A more fully featured and documented version of this model, which can accommodate multiple types of case data at once (e.g., reported cases, hospitalizations, ICU admissions) and estimate segments of positive-case sampling fractions for reported cases, is also available as an R package [covidseir](https://github.com/seananderson/covidseir).

Generally, any part of the analysis can be re-created by running one of the numbered R files starting with `01-...R` in the [`analysis`](analysis) folder. Alternatively, the file [`00-run-all.R`](analysis/00-run-all.R) can be sourced to run the entire analysis.

You will need the following packages installed:

```{r, eval=FALSE}
install.packages(c("tidyverse", "remotes", "rstan", "here", 
  "future", "deSolve", "furrr", "cowplot", "reshape2"))
remotes::install_github("seananderson/ggsidekick")
```

See the C++ compiler [installation instructions for covidseir](https://github.com/seananderson/covidseir) first, then:

```{r, eval=FALSE}
remotes::install_github("seananderson/covidseir", ref = "preprint")
```

An example of how to run the model:

```{r options, echo = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "README-figs/",
  cache.path = "README-cache/"
)
```

```{r load-data, message=FALSE, warning=FALSE}
library("rstan")
library("dplyr")
library("ggplot2")
rstan_options(auto_write = TRUE) # cache the compiled model
options(mc.cores = parallel::detectCores() / 2) # Stan parallel processing
seeiqr_model <- rstan::stan_model("analysis/seeiqr.stan")

source("analysis/fit_seeiqr.R")
source("analysis/make_projection_plot.R")

d <- readr::read_csv("data-generated/daily-cases.csv")
d
```

```{r fit-model, results='hide', message=FALSE, warning=FALSE}
# Using fewer iterations for a quick example:
fit <- fit_seeiqr(d$cases, seeiqr_model = seeiqr_model,
  iter = 300, chains = 4)
```

```{r summary}
print(fit$fit, pars = c("R0", "f2", "phi"))
```

```{r proj-plot}
make_projection_plot(list(fit)) + theme_light()
```