-
-
Save mok0/ae122af7006791af192b46ba3f561855 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
## Infrastructure notes ######################################################## | |
# I am running this in docker image rocker/tidyverse:3.6.1 | |
# Load packages | |
library(tidyverse) | |
library(EpiEstim) | |
## Fetch data from the covid19data API ######################################### | |
url <- "https://api.covid19data.dk/ssi_newly_hospitalized" | |
hosp_raw <- jsonlite::fromJSON(url) | |
hosp_raw <- hosp_raw %>% | |
as_tibble() %>% | |
mutate(date = date %>% lubridate::ymd_hms() %>% as.Date()) | |
## Plot the raw data ########################################################### | |
hosp_raw %>% | |
ggplot() + | |
geom_bar(aes(date, newly_hospitalized), stat = "identity", fill = "Purple") + | |
theme_minimal() | |
## Do data transformations ##################################################### | |
# The last four days are usually not fully up-to-date and are missing | |
# hospitalizations. So we will cut the last four entries and create a 7-day | |
# rolling mean | |
ads <- hosp_raw %>% | |
slice(1:(n()-4)) %>% # remove last four entries | |
mutate(rollingmean = c(rep(NA, 3), # Adding a rolling mean | |
zoo::rollmean(x = newly_hospitalized, k = 7), | |
rep(NA, 3))) | |
## Plot the new data ########################################################### | |
ads %>% | |
ggplot() + | |
geom_bar(aes(date, newly_hospitalized), stat = "identity", fill = "Purple") + | |
geom_line(aes(date, rollingmean)) + | |
theme_minimal() | |
## Choosing parameters for estimating R ######################################## | |
# We need to supply the estimate_R function with some parameters in order for it | |
# to estimate R. We will use two sets of parameters: | |
# 1) from the CDC where they have estimated mean and the standard deviation for | |
# the serial interval. | |
# 2) Numbers I have learned that the danish SSI uses. I have only mean informed | |
# of their mean, so the standard deviation I have chosen to be the same as | |
# the mean. This number can be changed as we get more info from SSI. | |
# CDC | |
# https://wwwnc.cdc.gov/eid/article/26/6/20-0357_article | |
cdc_mean_si = 3.96 | |
cdc_std_si = 4.75 | |
# SSI | |
ssi_mean_si = 4.7 | |
ssi_std_si = 4.7 | |
## Estimating R ################################################################ | |
confirmed_cases <- ads %>% drop_na() %>% select(I = rollingmean) | |
cdc_R <- estimate_R(confirmed_cases, | |
method = "parametric_si", | |
config = make_config(list(mean_si = cdc_mean_si, | |
std_si = cdc_std_si))) | |
ssi_R <- estimate_R(confirmed_cases, | |
method = "parametric_si", | |
config = make_config(list(mean_si = ssi_mean_si, | |
std_si = ssi_std_si))) | |
## Plotting the estimated R #################################################### | |
# First we create a data set that is ready for plotting | |
cdc_pd <- cdc_R$R %>% | |
as_tibble() %>% | |
select(t = t_end, | |
R = `Mean(R)`, | |
lower = `Quantile.0.05(R)`, | |
upper = `Quantile.0.95(R)`) %>% | |
mutate(source = "cdc") | |
ssi_pd <- ssi_R$R %>% | |
as_tibble() %>% | |
select(t = t_end, | |
R = `Mean(R)`, | |
lower = `Quantile.0.05(R)`, | |
upper = `Quantile.0.95(R)`) %>% | |
mutate(source = "ssi") | |
pd <- bind_rows(cdc_pd, ssi_pd) %>% | |
mutate(t = (t-1) + as.Date("2020-02-28")) | |
# Plot it | |
facetlabs <- c(cdc = "Center for Disease Control", | |
ssi = "Statens Serum Institut") | |
ggplot(pd) + | |
geom_ribbon(aes(x = t, ymin = lower, ymax = upper), fill = "lightgrey") + | |
geom_line(aes(x = t, y = R)) + | |
geom_hline(yintercept = 1, linetype = "dashed") + | |
facet_wrap(~ source, ncol = 1, | |
labeller = labeller(source = facetlabs)) + | |
theme_minimal() + | |
labs(x = "", y = "Effektiv reproduktionstal") + | |
scale_x_date(breaks = "weeks", labels = scales::date_format("%d-%m")) | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment