estimate_R0.R

library(dplyr)
library(tidyr)
library(ggplot2)
library(R0)
library(EpiEstim)

# read data
url_confirmed <- "https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv"
corona <- data.table::fread(url_confirmed) %>% as_tibble()
corona[c("Province/State", "Lat", "Long")] <- NULL

corona <-
  pivot_longer(corona, cols = -c(`Country/Region`), names_to = "date", values_to = "cases") %>%
  rename(country = `Country/Region`) %>%
  filter(country %in% c("Spain")) %>%
  mutate(date = as.Date(date, format = "%m/%d/%y"))


# set cases manually (non-updated data)
corona[corona$date == "2020-03-12", ]$cases <- 3033
corona[corona$date == "2020-03-23", ]$cases <- 35212

corona <-
  corona %>%
  mutate(daily_cases = cases - lag(cases))


espana <-
  corona %>%
  filter(!is.na(daily_cases), daily_cases > 10) %>%
  dplyr::select(daily_cases) %>%
  mutate(index = row_number())


# Model
mGT <- generation.time("gamma", c(3, 1.5))
estR0espana <- estimate.R(
  epid = as.integer(espana$daily_cases), GT = mGT, begin = 1, end = 31, methods = c("TD"),
  pop.size = 47000000, nsim = 500
)



# Estimates results tibble
tst_esp <- tibble::tibble(
  R = estR0espana$estimates$TD$R,
  inferior = estR0espana$estimates$TD$conf.int$lower,
  superior = estR0espana$estimates$TD$conf.int$upper,
  fecha = corona %>%
    filter(!is.na(daily_cases), daily_cases > 10) %>%
    dplyr::pull(date)
)


tst_esp <- tst_esp[1:nrow(tst_esp) - 1, ]

p1 <-
  ggplot(tst_esp, aes(fecha)) +
  geom_ribbon(aes(ymin = inferior, ymax = superior), fill = "grey", alpha = .5) +
  geom_line(aes(y = R)) +
  geom_point(aes(y = R)) +
  theme_minimal() +
  ylim(0, 5) +
  geom_hline(
    yintercept = 1, linetype = "dashed",
    color = "red", size = .5
  ) +
  labs(title = "Basic reproduction number in Spain") +
  theme(
    plot.title = element_text(size = 20),
    plot.subtitle = element_text(size = 17),
    axis.text = element_text(size = 14),
    axis.title = element_text(size = 14)
  )


ggsave(plot = p1, filename = "~/Descargas/basic_reproduction_number.png", width = 14, height = 10)


###

# Based on EpiEstim package
# serial interval for COVID-19: https://www.medrxiv.org/content/10.1101/2020.02.19.20025452v4
# https://cran.r-project.org/web/packages/EpiEstim/vignettes/demo.html
# http://spatial-r.com/en/2015/07/R0/

res <- estimate_R(espana$daily_cases, 
           method="parametric_si",
           config = make_config(list(
             mean_si = 3.96, 
             std_si = 4.75,
             t_start = c(2:29),
             t_end = c(2:29)+2)
))

fechas <- corona %>%
  filter(!is.na(daily_cases), daily_cases > 10) %>%
  dplyr::pull(date)

fechas <- fechas[-c(1:3)]

tst_esp <- tibble::tibble(
  R = res$R$`Mean(R)`,
  inferior = res$R$`Quantile.0.025(R)`,
  superior = res$R$`Quantile.0.975(R)`,
  fecha = fechas
)


p1 <-
  ggplot(tst_esp, aes(fecha)) +
  geom_ribbon(aes(ymin = inferior, ymax = superior), fill = "grey", alpha = .5) +
  geom_line(aes(y = R)) +
  geom_point(aes(y = R)) +
  theme_minimal() +
  ylim(0, 5) +
  geom_hline(
    yintercept = 1, linetype = "dashed",
    color = "red", size = .5
  ) +
  labs(title = "Basic reproduction number in Spain") +
  theme(
    plot.title = element_text(size = 20),
    plot.subtitle = element_text(size = 17),
    axis.text = element_text(size = 14),
    axis.title = element_text(size = 14)
  )


ggsave(plot = p1, filename = "~/Descargas/basic_reproduction_number.png", width = 14, height = 10)