rm(list = ls())
library(patchwork)
library(dplyr)
library(ggplot2)
remotes::install_github("Bisaloo/sir_age", upgrade = TRUE)
library(sirage)
remotes::install_github("Bisaloo/Npieurope", upgrade = TRUE)
library(NpiEurope)
temp <- read.csv(system.file("extdata", "COVID_time_series_v4_2020-09-16.csv", package = "NpiEurope"),
stringsAsFactors = FALSE) %>%
filter(Date <= "2020-09-08")
countryVec <- sort(unique(temp$Country))
temp <- temp %>%
group_by(Country) %>%
mutate(cum_cases = cumsum(NewCases),
simulated = NA_real_)
folder <- "MCMCOK5"
pGraph <- list()
for (i in 1:length(countryVec)) {
country <- countryVec[i]
print(country)
country_data <- NpiEurope::load_country_data(country)
country_data$NewCases[which(country_data$NewCases < 0)] <- -country_data$NewCases[which(country_data$NewCases < 0)]
epi_data <- country_data[, c("Date", "NewCases", "NewDeaths")]
npi_data <- country_data[, !colnames(country_data) %in% c("NewCases", "NewDeaths", "Country", "Population")]
epi_data <- estimate_asympto(epi_data, smooth = 7)
contact_data <- NpiEurope::load_contact_data(country)
age_data <- load_age_data(country)
plot_asympto(epi_data)
ggsave(sprintf("%s/Report_%s.pdf", folder, country))
res <- read.csv(sprintf("%s/%s.csv", folder, country))
plot_MCMC(sprintf("%s/%s.csv", folder, country))
ggsave(sprintf("%s/mcmc_%s.pdf", folder, country))
write.csv(
summarise_estimation(sprintf("%s/%s.csv",folder,country), npi_data, contact_data, age_data),
sprintf("%s/estim_%s.csv",folder,country)
)
# knitr::kable(summarise_estimation(sprintf("MCMCNew/MCMC_%s.csv",country), country_data))
sims <- simulate_trajectory(sprintf("%s/%s.csv", folder, country),
epi_data, npi_data, contact_data, age_data,
Npost = 5)
sims <- t(do.call(rbind, sims))
temp$simulated[temp$Country==country & !is.na(temp$cum_cases) & temp$cum_cases > 0] <- apply(sims, MARGIN = 2, median)
tempG <- plot_simulations(sims, country_data) + ggtitle(country)
pGraph[[length(pGraph) + 1]] <- tempG
ggsave(sprintf("%s/sim_%s.pdf", folder, country))
}
p <- ggplot(temp, aes(x = NewCases, y = simulated, color = Country)) +
geom_point() +
scale_x_continuous(trans = "log10") +
scale_y_continuous(trans = "log10") +
geom_abline(slope = 1, col = "black", linetype = "dashed") +
xlab("Observed cases") + ylab("Simulated cases") +
theme_minimal()
# pGraph[[length(pGraph) + 1]] <- p
#
# layout <- mapply(area, t = c(rep(1:4, each = 6), rep(5:6, each = 4)), l = c(rep(1:6, times = 4), rep(1:4, times = 2)), SIMPLIFY = FALSE)
# layout <- do.call(patchwork:::c.patch_area, layout)
#
# layout <- c(layout, area(5,5,6,6))
#
# wrap_plots(pGraph, design = layout)
wrap_plots(pGraph)
ggsave("figure2.pdf", width = 11.69, height = 8.27)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.