Skip to content

Instantly share code, notes, and snippets.

@klpn
Last active Feb 28, 2021
Embed
What would you like to do?
Plot regional or national Rt for COVID-19 based on FHM data to reproduce FHM reports
# Plot regional or national Rt for COVID-19 based on FHM data, to reproduce the
# results published on https://www.folkhalsomyndigheten.se/smittskydd-beredskap/utbrott/aktuella-utbrott/covid-19/statistik-och-analyser/analys-och-prognoser/
# Adapted from:
# Churches (2020, Feb. 18). Tim Churches Health Data Science Blog: Analysing COVID-19 (2019-nCoV) outbreak data with R - part 1. Retrieved from https://timchurches.github.io/blog/posts/2020-02-18-analysing-covid-19-2019-ncov-outbreak-data-with-r-part-1/, CC BY-SA
library(readxl)
library(EpiEstim)
# path to latest FHM Excel file
secov <- read_excel("data/FHM/FHM_latest.xlsx")
secov$Statistikdatum <- as.Date(secov$Statistikdatum)
regframe <- function(reg, sdate, edate) {
df <- data.frame(secov$Statistikdatum, secov[[reg]])
colnames(df) <- c("dates", "I")
return(df[df$dates>=sdate & df$dates<=edate, ])
}
secov_par_si <- function(rf) {
return(estimate_R(rf, method = "parametric_si",
config = make_config(list(mean_si = 4.8, std_si = 2.3))))
}
secov_par_si_u <- function(rf) {
return(estimate_R(rf, method = "uncertain_si",
config = make_config(list(mean_si = 4.8, std_mean_si = 2, min_mean_si = 1.2,
max_mean_si = 8.4, std_si = 2.3, std_std_si = 1,
min_std_si = 0.6, max_std_si = 4, n1 = 100, n2 = 100))))
}
# plot_Ri(secov_par_si(regframe("Uppsala", "2020-09-01", "2021-02-25")))
plot_Ri <- function(estimate_R_obj) {
p_I <- plot(estimate_R_obj, "incid", add_imported_cases = TRUE) # plots the incidence
p_SI <- plot(estimate_R_obj, "SI") # plots the serial interval distribution
p_Ri <- plot(estimate_R_obj, "R")
return(gridExtra::grid.arrange(p_I, p_SI, p_Ri, ncol = 1))
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment