COVID_Analysis.Rmd

---
title: "EID Class - Fall 2021"
author: "Grant Brown"
date: "9/28/2021"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


This document contains the basic structure of a stochastic, Bayesian epidemic model for COVID-19 in Iowa. We'll go through the steps in class. 



```{r}
# Load packages:
suppressPackageStartupMessages({
  suppressWarnings({
    library(ABSEIR)
    library(ggplot2)
    library(knitr)
    library(dplyr)
    library(openxlsx)
    library(splines)
    library(lubridate)
  })
})

#  https://github.com/nytimes/covid-19-data
NYT_state_data <- read.csv("https://raw.githubusercontent.com/nytimes/covid-19-data/master/us-states.csv")


# Iowa Population
state_population <- 3155000
stateDataFiltered <- filter(NYT_state_data, state == "Iowa") %>% arrange(date) %>%
    mutate(date = as.Date(date,format = "%Y-%m-%d")) %>% 
    filter(weekdays(date) == "Friday") %>% 
    mutate(idxDate = as.numeric(date - min(date))) %>%
    mutate(deathsNonCum = c(deaths[1], diff(deaths))) %>% 
    mutate(casesNonCum = c(cases[1], diff(cases))) %>%
    select(date, cases, deaths, idxDate, deathsNonCum, casesNonCum)


minIdx <- max(1, min(which(stateDataFiltered$casesNonCum > 0))-7)
stateDataFiltered <- stateDataFiltered[minIdx:nrow(stateDataFiltered),]

stateDataFiltered$idxDate <- as.numeric(stateDataFiltered$date -
                                        min(stateDataFiltered$date))

# Focus just on early period
stateDataFiltered <- filter(stateDataFiltered, 
                            date >= as.Date("2020-01-01"))
par(xaxt="n")
plot(stateDataFiltered$casesNonCum~stateDataFiltered$date)
par(xaxt="s")
abline(v = as.Date("2020-01-01") + months(1:24))
axis(side = 1, at = as.Date("2020-01-01") + months(1:24), 
     labels = months(as.Date("2020-01-01") + months(1:24)))




         
```

Model Based on Deaths
=============================

```{r}


# Build I and R data models
data_model_1 = DataModel(c(stateDataFiltered$deathsNonCum, rep(NA, 90)),
                         type = "fractional",
                         compartment="I_star",
                         cumulative=FALSE, 
                         params = list(report_fraction = 0.75,   # Approximate CFR
                                       report_fraction_ess = 200))
# Build Exposure Model

#interventionDate <- as.Date("2020-04-04")
#reopenDate <- as.Date("2020-05-01")

X <- cbind(1, ns(1:(nrow(stateDataFiltered)+90), df=6))
exposure_model = ExposureModel(X,
                               nTpt = (nrow(stateDataFiltered)+90),
                               nLoc = 1,
                               betaPriorPrecision = 0.5,
                               betaPriorMean = 0)
# There's no reinfection in this case, so we just use a "SEIR" model. 
reinfection_model = ReinfectionModel("SEIR")

# we have no distance model, because it's a non-spatial analysis
distance_model = DistanceModel(list(matrix(0)))

# Set initial population sizes
initial_value_container = InitialValueContainer(S0=state_population,
                                                E0=10,
                                                I0=10,
                                                R0=0, type = "uniform", 
                                                params = list(max_S0 = state_population + 100,
                                                              max_E0 = 1000,
                                                              max_I0 = 1000,
                                                              max_R0 = 2))


# Model to describe E to I and I to R transition probabilities.

# Latent period: 2-14 days with median 5

pickWeibullPars <- function(qdf){
  rslt <- optim(par = c(1,1), fn = function(par){
    sum((qweibull(p = qdf$q, shape = par[1], scale = par[2]) - qdf$x)^2)
  })
  rslt$par
}

pickGammaHyperPars <- function(mean, ESS){
  b <- ESS/(mean+1)
  a <- ESS - b
  c(a,b)
}

latent_par_means <- pickWeibullPars(qdf=data.frame(q=c(0.025,0.5,0.975),
                           x=c(2,5,14)))
infectious_par_means <- pickWeibullPars(qdf = data.frame(q=c(0.025,0.5,0.975),
                                                          x = c(10,14,32)))

Transition_priors = WeibullTransitionPriors(latent_shape_prior_alpha = pickGammaHyperPars(latent_par_means[1], 1000)[1],
                                           latent_shape_prior_beta = pickGammaHyperPars(latent_par_means[1], 1000)[2],
                                           latent_scale_prior_alpha = pickGammaHyperPars(latent_par_means[2], 1000)[1],
                                           latent_scale_prior_beta = pickGammaHyperPars(latent_par_means[2], 1000)[2],
                                           infectious_shape_prior_alpha = pickGammaHyperPars(infectious_par_means[1], 100)[1],
                                           infectious_shape_prior_beta = pickGammaHyperPars(infectious_par_means[1], 100)[2],
                                           infectious_scale_prior_alpha = pickGammaHyperPars(infectious_par_means[2], 100)[1],
                                           infectious_scale_prior_beta = pickGammaHyperPars(infectious_par_means[2], 100)[2]) 
                                           
sampling_control = SamplingControl(seed = 123124, 
                                   n_cores = 6,
                                   algorithm="Beaumont2009",
                                   list(init_batch_size = 100000,
                                        batch_size = 2000,
                                        epochs = 1000,
                                        max_batches = 20,
                                        shrinkage = 0.95,
                                        keep_compartments =TRUE
                                   ))

if (!file.exists(results.rda)){
  result = SpatialSEIRModel(data_model_1,
                            exposure_model,
                            reinfection_model,
                            distance_model,
                            transition_priors = Transition_priors,
                            initial_value_container,
                            sampling_control,
                            samples = 100,
                            verbose = 2)
  
  save("result", file = "results.rda", compress="bzip2")
} else{
  load("results.rda")
}
summary(result)

```

Summary
=======================================

```{r}

# Ugly function for plotting
plotPosteriorPredictive = function(sim_obj, dat, parName, useRF=FALSE, 
                                   finish = function(){},maxT = NA,
                                   heightQtile= 0.8,
                                   ...)
{
  
  allSimulated <-  sapply(sim_obj$simulationResults,
                          function(x){(x[[parName]])})
  if (useRF){
    rf <- getResultSummary(sim_obj,parname = "report_fraction", f = c)
    allSimulated <- allSimulated*matrix(rf, nrow = nrow(allSimulated),
                                        ncol = ncol(allSimulated), byrow = TRUE)
  }
  
  
  if (class(maxT) != "logical"){
    allSimulated <- allSimulated[1:maxT,]
  } else{
    maxT <- nrow(allSimulated)
  }
  
  Nt <- dim(allSimulated)[1]
  
  upperMQuantile = apply(allSimulated, 1, quantile, probs = c(heightQtile))

  
  if (length(dat$date[1:maxT]) == Nt){
    dateObj <- dat$date[1:maxT]
  } else {
    dateObj <- 1:Nt
  }
  if(any(is.na(dateObj))){
    dateObj[is.na(dateObj)] <- max(dateObj, na.rm=TRUE) + cumsum(rep(7, sum(is.na(dateObj))))
  }
  
  layout(mat = matrix(c(1,2), nrow = 1), widths = c(8,2), heights = c(6,6))
  par(xaxt="n")
  plot(dateObj, 1:Nt, 
       ylim = c(0, max(upperMQuantile)), type = "n",
       ...)
  par(xaxt="s")

  abline(v = as.Date("2020-01-01") + months(1:24))
  axis(side = 1, at = as.Date("2020-01-01") + months(1:24), 
       labels = months(as.Date("2020-01-01") + months(1:24)))

  
  #gridV <- (ceiling(seq(0, par("usr")[4], length = 20)/10)*10)[2]
  ymax <- (round(par("usr")[4]))
  gridV <- as.numeric(
    paste0(c("1", rep("0", nchar(as.character(ymax)) - 1)), 
           collapse = ""))
  while (ymax/gridV < 5){
    gridV <- gridV/2  
  }
  
  
  #qtls <- seq(0.025, 0.975, length = 10)
  qtls <- c(0.025,0.1,0.25,0.5,0.75,0.9, 0.975)
  qtlwd <- 3-2*sqrt(abs(qtls - 0.5))
  qtlty <- c(4,3,2,1,2,3,4)
  qtalpha <- 1-sqrt(abs(qtls - 0.5))
  qtcol <- rgb(0,0,1,qtalpha)
  
  for (i in 1:length(qtls)){
    qtileCurve = apply(allSimulated, 1, quantile, probs = qtls[i])
    lines(dateObj, qtileCurve, lty = qtlty[i], lwd = qtlwd[i], 
          col = qtcol[i])
  }
  qtileMCurve = apply(allSimulated, 1, quantile, probs = c(0.5))
  lines(dateObj, qtileMCurve, lty = 1, lwd = 3, col = rgb(0,0,1,1))
    
  
  
  abline(v = Sys.Date(), lty = 2, col = "green")

  abline(h = seq(0, max(allSimulated), gridV), col = "lightgrey", lty = 2)
  finish()
  
  ## Add Legend
  
  par(xpd = TRUE)
  par(bty = "n")
  par(xaxt = "n")
  par(yaxt = "n")
  
  
  plot(0,0,xlim = c(0,2), ylim = c(0,6), type = "n",
       xlab = "", ylab = "", main = "")
  legend(x = 0, y = 3, legend = rev(as.character(qtls)), 
         lty = qtlty,lwd = qtlwd, col = qtcol,
         xjust = 0.5,
         yjust = 0.5)
  
  par(yaxt = "s")
  par(xaxt = "s")
  par(bty = "o")
  par(xpd = FALSE)
  return(allSimulated)
}

r1 <- plotPosteriorPredictive(sim_obj = result, 
                              dat = stateDataFiltered,
                          parName = "I", main = "COVID Model Based on Cases, I", 
                          maxT = nrow(stateDataFiltered) + 9)

r2 <- plotPosteriorPredictive(sim_obj = result, 
                              dat = stateDataFiltered,
                          parName = "S", main = "COVID Model Based on Cases, S",
                          maxT = nrow(stateDataFiltered) + 9)
r3 <- plotPosteriorPredictive(sim_obj = result, 
                              dat = stateDataFiltered,
                          parName = "I", main = "COVID Model Based on Cases, E_star",
                          maxT = nrow(stateDataFiltered) + 9)

```