June 20, 2023
You will need to install the package, “pfilterCOVID”, before using the pfilter function in the document
# devtools::load_all()
# devtools::install()
devtools::install_github("kimfinale/pfilterCOVID", force = TRUE)library(pfilterCOVID)
library(deSolve)
library(tidyverse)
library(cowplot)
library(EpiEstim)
library(foreach)
library(doParallel)
source('R/deconvolve.R')
theme_set(theme_bw())# beta is updated by R0
THETA <- list(beta = 0.3, R0 = 1.2, epsilon = 1/2.5, delta = 1/5,
gamma = 1/2.5, bp = 1, ba = 1, fa = 0.3, fd = 1, I0 = 1e2,
time_dep_Rt = FALSE)
# This changes the default y0 in the pfilterCOVID package
THETA$Y0 <- c(S = 1e7 - THETA$I0, E = 0, P = 0, A = 0,
I = THETA$I0, R = 0, CE = 0, CI = 10, CR = 0)
THETA$tend <- 200
usethis::use_data(THETA, overwrite = TRUE)
params <- THETAparams <- THETA
tend <- params$tend
params$I0 <- 100
params$Y0 <- c(S = 1e7 - params$I0, E = 0, P = 0, A = 0,
I = params$I0, R = 0, CE = 0, CI = 0, CR = 0)
tstamp <- format(Sys.time(), "%Y%m%d")
times <- seq(0, tend, by = 1)
params$time_dep_Rt <- TRUE
ode(y = params$Y0, times = times, func = ode_sepair, parms = params) %>%
as.data.frame() -> out
out$daily_infected <- c(0, diff(out$CE))
out$daily_symptom <- c(0, diff(out$CI))
out$daily_confirmed <- c(0, diff(out$CR))
out$daily_Rt <- sapply(0:tend, function(x) get_Rt(x))
# saveRDS(out, paste0("outputs/ode_I0_100_", tstamp, ".rds"))tstamp <- "20230619"
out <- readRDS(paste0("outputs/ode_I0_100_", tstamp, ".rds"))
out_long <- out %>% pivot_longer(-time)
out_long$name <- factor(out_long$name,
levels = c("S", "E", "P", "A", "I", "R",
"CE", "CI", "CR", "daily_infected",
"daily_symptom", "daily_confirmed", "daily_Rt"))
fig2a <- out_long %>%
filter(name == "daily_infected" |
name == "daily_symptom" | name == "daily_confirmed") %>%
ggplot(aes(x = time, y = value, color = name)) +
geom_line(linewidth=1.2) +
labs(x='Time (day)', y = 'Number of individuals', color = "") +
theme(legend.position = c(.2, .9),
legend.background = element_rect(fill='transparent')) +
scale_color_hue(labels=c("daily_infected" = "Infection",
"daily_symptom"="Symptom",
"daily_confirmed"="Confirmation"))+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))Plots of Rt and infection/confirmation time series of stoch outputs NB. Figure 1 is a compartment flow diagram which was not generated by code
# Rt data set; Rt_dataset_I0=100_20220105.rds, Rt_dataset_I0=10_20220105.rds
tstamp <- "20220105"
dat <- readRDS(paste0("outputs/Rt_dataset_I0=100_", tstamp, ".rds"))
# Generate Figure 2B.
# plot for the pre-defined Rt
fig2b <- ggplot() +
geom_line(aes(x=dat$time, y = dat$Rt), linewidth=1.4) +
labs(x="Time (day)", y = expression(paste(italic(R)[italic(t)])))+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text = element_text(size=13))
figure2 <- plot_grid(fig2a, fig2b, nrow = 2, labels = "AUTO")
figure2
# Generate plots for a comparison of simulated data (deterministic vs. stochastic, both with perfect observation) for daily_infected or confirmed.
oderes <- dat$ode[, c("time", "daily_infected", "daily_confirmed")]
oderes %>% pivot_longer(cols = -time) -> oderes
stochres <- dat$stoch_perfect_obs$daily_confirmed[,1:length(dat$stoch_perfect_obs$daily_confirmed)]
stochres$time <- dat$time
stochres %>% pivot_longer(cols = -time) -> stochres
# plt_case <- ggplot() +
# geom_line(data = stochres, aes(time, value, group=name), color="grey50") +
# geom_line(data = oderes, aes(time, value, color=name), size=1.2) +
# labs(x="time", y="", color="")This is to create a stochastic data set based on the Gillespie’s direct method
set.seed(42)
# tstamp <- format(Sys.time(), "%Y%m%d")
nrun <- 20 # increase the number to obtain large samples
res <- gillespie_run(func = gillespie_sepair,
tend = tend,
nrun = nrun,
y0 = params$Y0,
params = params,
report_dt = 1)
# saveRDS(res, "outputs/res_20230619.rds")res <- readRDS("outputs/res_20230619.rds")
nrun <- 20
tend <- params$tend
daily_infected <- data.frame(matrix(0, nrow = tend + 1, ncol = nrun))
daily_symptom <- data.frame(matrix(0, nrow = tend + 1, ncol = nrun))
daily_confirmed <- data.frame(matrix(0, nrow = tend + 1, ncol = nrun))
for (i in seq_len(length(res))) {
inf <- diff(res[[i]]$CE)
onset <- diff(res[[i]]$CI)
conf <- diff(res[[i]]$CR)
daily_infected[1:(length(inf)+1), i] <- c(0, inf)
daily_symptom[1:(length(onset)+1), i] <- c(0, onset)
daily_confirmed[1:(length(conf)+1), i] <- c(0, conf)
}
daily_infected <- data.frame(matrix(0, nrow = tend + 1, ncol = nrow(dat$stoch_perfect_obs$daily_infected)))
for (i in 1:nrow(dat$stoch_perfect_obs$daily_infected)) {
daily_infected[,i] <- dat$stoch_perfect_obs$daily_infected[[i]]
}
stoch <- list()
stoch$daily_confirmed <- daily_confirmed
stoch$daily_symptom <- daily_symptom
stoch$daily_infected <- daily_infected
stoch$daily_Rt <- sapply(0:tend, function(x) get_Rt(x))
# "stoch" contains daily time series for infection, symptom, and confirmation (corresponds to the dataset, D2, in the manuscript)
# saveRDS(stoch, paste0("outputs/stoch_I0_100_", tstamp, ".rds"))
df <- daily_infected
df$time <- 0:tend
df %>% pivot_longer(cols = -time) -> df
# Generate Figure 3.
locy = c(2000,1920,1840)
fig3 <- ggplot() +
geom_line(data = df, aes(time, value, group = name),
color = "grey80") +
geom_line(data = out, aes(time, daily_infected), # ODE model output
color = "darkred",
size = 1.5, inherit.aes = FALSE) +
labs(x="Time (day)", y="Daily infected cases", color="") +
geom_segment(aes(x=c(0,0,0), xend=c(5,5,5), y=locy,
yend=locy), size=c(1,1,1),
color = c("grey80", "darkred","black"),
linetype = c("solid","solid","dotted")) +
geom_line(aes(dat$time, rowMeans(dat$stoch_perfect_obs$daily_infected)), # mean of stochastic simulations ## simply to load the data we are using (2,000 simulations) - save time
linetype = "dotted",
size = 1.5, inherit.aes = FALSE) +
geom_text(aes(x=c(7,7,7), y=locy), size=c(4.6,4.6,4.6),
hjust = c(0,0,0), vjust = c(0.5,0.5,0.5),
label = c("Stochastic simulations",
"Determinstic simulation",
"Mean of stochastic simulations"))+
theme(text = element_text(size=16),
axis.text = element_text(size=15),
legend.text=element_text(size=15))
fig3
# ggsave("plots/daily_case_20230619.png", plot=figure2,
# width=3.4*3, height=2.7*3, units="in")The data set that contains all the necessary information: Rt trajectory, ODE and stochastic model result
Rt_dataset <- list()
Rt_dataset$README <- "To see if potential bias of Rt estimation based on the confirmation data set, especially in the most recent time, is reduced or eliminated when the number of cases are large"
Rt_dataset$ode <- out
Rt_dataset$stoch_perfect_obs <- stoch
Rt_dataset$Rt <- sapply(0:tend, function(x) get_Rt(x))
Rt_dataset$time <- out$time
Rt_dataset$params <- params
Rt_dataset$y0 <- params$Y0
# saveRDS(Rt_dataset, paste0("outputs/Rt_dataset_I0=100_", tstamp, ".rds"))Particle filtering - particle_filter Expect errors since backward sampling is not in place
set.seed(1)
tstamp <- "20220105"
I0 <- 100
dat <- readRDS(paste0("outputs/Rt_dataset_I0=", I0, "_", tstamp, ".rds"))
# determine the data type to model (infection vs. confirmation)
dtype <- "confirmation"
# dtype <- "infection"
if (dtype == "confirmation") {
d <- dat$ode[, c("time", "daily_confirmed")]
} else {
d <- dat$ode[, c("time", "daily_infected")]
}
dat$params$time_dep_Rt <- FALSE # this must be set FALSE to produce daily estimates
params <- dat$params
params[["betavol"]] <- 0.2 #standard deviation of the prior distribution
npart = 1e4 # increase the number of particles to obtain more accurate estimates
pf <- pfilter(
params = params,
y = params$Y0,
data = d,
data_type = dtype,
npart = npart,
tend = nrow(d),
dt = 0.1,
error_pdf = "pois",
systematic_resampling = FALSE,
backward_sampling = TRUE,
stoch = FALSE)
# Rt <- readRDS("outputs/Rt_default.rds")
epsilon <- params[["epsilon"]]
delta <- params[["delta"]]
gamma <- params[["gamma"]]
bp <- params[["bp"]] # relative infectiousness of pre-symptomatic state
ba <- params[["ba"]] # relative infectiousness of asymptomatic state
fa <- params[["fa"]] # fraction of asymptomatic state
durP <- (1 / delta - 1 / epsilon)
durI <- (1 / gamma)
R0_dur <- ((1 - fa) + fa * ba) * durI + bp * durP
med <- apply(pf$beta_filtered, 1, quantile, probs = c(0.5))
# plot(dat$time + 2, dat$daily_Rt, type="l", ylim=c(0, max(upper * R0_dur)))
# Comarison between the pre-determined Rt and PF (without backward sampling) inferred Rt based on ODE daily_confirmed
plot(dat$ode$time + 2, dat$ode$daily_Rt, type="l")
lines(med * R0_dur, lwd=2, col=3)This uses the extract_trace function that calls pfilter function and also executes backward sampling
# apply on infection time series
tstamp <- "20220105"
dat <- readRDS(paste0("outputs/Rt_dataset_I0=100_", tstamp, ".rds"))
params <- dat$params
params$time_dep_Rt <- FALSE
params[["betavol"]] <- 0.2
# increase the number to improve accuracy
nrep <- 1e3
npart <- 1e4
dt <- 0.2
library(parallel)
library(doParallel)
ncores <- detectCores() - 1
d_inf <- data.frame(date = dat$time, daily_infected = round(dat$ode$daily_infected))
set.seed(42)
cl <- makeCluster(getOption("cl.cores", ncores))
doParallel::registerDoParallel(cl)pf_inf <- foreach(i = 1:nrep, .packages = c("pfilterCOVID"), .inorder = F) %dopar% {
extract_trace(params = params,
y = params$Y0,
data = d_inf,
data_type = "infection",
npart = npart,
tend = nrow(d_inf),
dt = dt,
error_pdf = "pois",
negbin_size = 15,
backward_sampling = TRUE,
stoch = FALSE)
}
parallel::stopCluster(cl)
# saveRDS(pf_inf, "outputs/pf_inf_20230619.rds")pf_inf <- readRDS("outputs/pf_inf_20230619.rds")
nrep <- 1e3
npart <- 1e4
dt <- 0.2
dtype <- "infection"
tstamp <- "20220105"
parset_chr <- paste0(dtype, "_I0=", params$I0, "_npart=", npart,
"_nrep=", nrep, "_dt=", dt, "_", tstamp)
pr <- c(0.025, 0.25, 0.5, 0.75, 0.975)
Rt_est <- as.data.frame(sapply(pf_inf, function(x) x[, "Rt"]))
Rt_quantile <- as.data.frame(t(apply(Rt_est, 1, function(x) quantile(x, pr))))
df_inf <- cbind(Rt_quantile, d_inf)
col_fill <- "deepskyblue4"
col_dat <- "grey70"
df_inf$Rt <- dat$ode$daily_Rt
fig4a <- ggplot(df_inf[81:201,], aes(x = date)) +
geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = col_fill, alpha = 0.5) +
geom_ribbon(aes(ymin = `25%`, ymax = `75%`), fill = col_fill, alpha = 0.8) +
geom_line(aes(y = `50%`), color = "darkblue", size = 1.2) +
geom_line(aes(x = date, y = dat$ode$daily_Rt[80:200]), size = 1, linetype = "dashed") +
geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+
labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") +
coord_cartesian(ylim = c(0.25, 2.25))+scale_y_continuous(breaks=seq(0,2.25,by=0.5))+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
# fig4a
# ggsave(paste0("plots/", parset_chr, ".png"), plt)# apply on confirmation time series
d_conf <- data.frame(date = dat$time,
daily_confirmed = round(dat$ode$daily_confirmed))set.seed(42)
cl <- makeCluster(getOption("cl.cores", ncores))
doParallel::registerDoParallel(cl)
pf_conf <- foreach(i = 1:nrep, .packages = c("pfilterCOVID"), .inorder = F) %dopar% {
extract_trace(params = params,
y = params$Y0,
data = d_conf,
data_type = "confirmation",
npart = npart,
tend = nrow(d_conf),
dt = dt,
error_pdf = "pois",
negbin_size = 15,
backward_sampling = TRUE,
stoch = FALSE)
}
parallel::stopCluster(cl)
# saveRDS(pf_conf, "outputs/pf_conf_20230619.rds") pf_conf <- readRDS("outputs/pf_conf_20230619.rds")
parset_chr <- paste0(dtype, "_I0=", params$I0, "_npart=", npart,
"_nrep=", nrep, "_dt=", dt, "_", tstamp)
Rt_est <- as.data.frame(sapply(pf_conf, function(x) x[, "Rt"]))
Rt_quantile <- as.data.frame(t(apply(Rt_est, 1, function(x) quantile(x, pr))))
df_conf <- cbind(Rt_quantile, d_conf)
df_conf$Rt <- dat$ode$daily_Rt
fig4b <- ggplot(df_conf[81:201,], aes(x = date)) +
geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = "orange", alpha = 0.5) +
geom_ribbon(aes(ymin = `25%`, ymax = `75%`), fill = "orange", alpha = 0.8) +
geom_line(aes(y = `50%`), color = "brown", size = 1.2) +
geom_line(aes(x = date, y = dat$ode$daily_Rt[80:200]), size = 1, linetype = "dashed") +
geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+
labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") +
coord_cartesian(ylim = c(0.25, 2.25))+scale_y_continuous(breaks=seq(0,2.25,by=0.5))+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
# fig4b
# put all together for figure 3
figure4 <- plot_grid(fig4a, fig4b, nrow = 2, labels = "AUTO")
figure4
# ggsave("plots/Rt_recovered_20230619.png", plot=figure3,
# width=3.4*2, height=2.7*2, units="in")tstamp <- "20220105"
dat <- readRDS(paste0("outputs/Rt_dataset_I0=100_", tstamp, ".rds"))
params <- dat$params
params$time_dep_Rt <- FALSE
params[["betavol"]] <- 0.2
nrep <- 1e3
npart <- 1e4
dt <- 0.2
d_stoc_19 <- data.frame(date = dat$time, daily_confirmed = round(dat$stoch_perfect_obs$daily_confirmed[,19]))library(parallel)
library(doParallel)
set.seed(42)
ncores <- detectCores()
cl <- makeCluster(getOption("cl.cores", 2))
doParallel::registerDoParallel(cl)
# apply on the stochastic seed #19 (fig 4b)
pf_stoc_19 <- foreach(i = 1:nrep, .packages = c("pfilterCOVID"), .inorder = F) %dopar% {
extract_trace(params = params,
y = params$Y0,
data = d_stoc_19,
data_type = "confirmation",
npart = npart,
tend = nrow(d_stoc_19),
dt = dt,
error_pdf = "pois",
negbin_size = 15,
backward_sampling = TRUE,
stoch = FALSE)
}
parallel::stopCluster(cl)
# saveRDS(pf_stoc_19, "outputs/pf_stoc_19_20230619.rds") pf_stoc_19 <- readRDS("outputs/pf_stoc_19_20230619.rds")
parset_chr <- paste0(dtype, "_I0=", params$I0, "_npart=", npart,
"_nrep=", nrep, "_dt=", dt, "_", tstamp)
pr <- c(0.025, 0.25, 0.5, 0.75, 0.975)
Rt_est <- as.data.frame(sapply(pf_stoc_19, function(x) x[, "Rt"]))
Rt_quantile <- as.data.frame(t(apply(Rt_est, 1, function(x) quantile(x, pr))))
df_stoc_19 <- cbind(Rt_quantile, d_stoc_19)
col_fill <- "red"
col_dat <- "grey70"
df_stoc_19$Rt <- dat$ode$daily_Rt
fig5b <- ggplot(df_stoc_19[81:201,], aes(x = date)) +
geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = col_fill, alpha = 0.5) +
geom_ribbon(aes(ymin = `25%`, ymax = `75%`), fill = col_fill, alpha = 0.8) +
geom_line(aes(y = `50%`), color = "darkred", size = 1.2) +
geom_line(aes(x = date + 1, y = Rt), size = 1, linetype = "dashed") +
geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+
labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") +
coord_cartesian(ylim = c(0.25, 2.75))+
scale_y_continuous(breaks=seq(0,2.75,by=0.5))+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
# fig5b# apply on the stochastic seed #04 (fig 4c)
d_stoc_04 <- data.frame(date = dat$time,
daily_confirmed = round(dat$stoch_perfect_obs$daily_confirmed[,4]))set.seed(42)
cl <- makeCluster(getOption("cl.cores", 2))
doParallel::registerDoParallel(cl)
pf_stoc_04 <- foreach(i = 1:nrep, .packages = c("pfilterCOVID"), .inorder = F) %dopar% {
extract_trace(params = params,
y = params$Y0,
data = d_stoc_04,
data_type = "confirmation",
npart = npart,
tend = nrow(d_stoc_04),
dt = dt,
error_pdf = "pois",
negbin_size = 15,
backward_sampling = TRUE,
stoch = FALSE)
}
parallel::stopCluster(cl)
# saveRDS(pf_stoc_04, "outputs/pf_stoc_04_20230619.rds")pf_stoc_04 <- readRDS("outputs/pf_stoc_04_20230619.rds")
parset_chr <- paste0(dtype, "_I0=", params$I0, "_npart=", npart,
"_nrep=", nrep, "_dt=", dt, "_", tstamp)
pr <- c(0.025, 0.25, 0.5, 0.75, 0.975)
Rt_est <- as.data.frame(sapply(pf_stoc_04, function(x) x[, "Rt"]))
Rt_quantile <- as.data.frame(t(apply(Rt_est, 1, function(x) quantile(x, pr))))
df_stoc_04 <- cbind(Rt_quantile, d_stoc_04)
col_fill <- "green"
col_dat <- "grey70"
df_stoc_04$Rt <- dat$ode$daily_Rt
fig5c <- ggplot(df_stoc_04[81:201,], aes(x = date)) +
geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = "#76ff7a", alpha = 0.5) +
geom_ribbon(aes(ymin = `25%`, ymax = `75%`), fill = "green", alpha = 0.8) +
geom_line(aes(y = `50%`), color = "darkgreen", size = 1.2) +
geom_line(aes(x = date + 1, y = Rt), size = 1, linetype = "dashed") +
geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+
labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") +
coord_cartesian(ylim = c(0.25, 2.5))+
scale_y_continuous(breaks=seq(0,2.5,by=0.5))+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
# fig5cIncidence level is lower than the deterministic model
# apply on the stochastic seed #02 (fig 4d)
d_stoc_02 <- data.frame(date = dat$time,
daily_confirmed = round(dat$stoch_perfect_obs$daily_confirmed[,2]))set.seed(42)
cl <- makeCluster(getOption("cl.cores", 2))
doParallel::registerDoParallel(cl)
pf_stoc_02 <- foreach(i = 1:nrep, .packages = c("pfilterCOVID"), .inorder = F) %dopar% {
extract_trace(params = params,
y = params$Y0,
data = d_stoc_02,
data_type = "confirmation",
npart = npart,
tend = nrow(d_stoc_02),
dt = dt,
error_pdf = "pois",
negbin_size = 15,
backward_sampling = TRUE,
stoch = FALSE)
}
parallel::stopCluster(cl)
# saveRDS(pf_stoc_02, "outputs/pf_stoc_02_20230619.rds") pf_stoc_02 <- readRDS("outputs/pf_stoc_02_20230619.rds")
parset_chr <- paste0(dtype, "_I0=", params$I0, "_npart=", npart,
"_nrep=", nrep, "_dt=", dt, "_", tstamp)
pr <- c(0.025, 0.25, 0.5, 0.75, 0.975)
Rt_est <- as.data.frame(sapply(pf_stoc_02, function(x) x[, "Rt"]))
Rt_quantile <- as.data.frame(t(apply(Rt_est, 1, function(x) quantile(x, pr))))
df_stoc_02 <- cbind(Rt_quantile, d_stoc_02)
col_fill <- "blue"
col_dat <- "grey70"
df_stoc_02$Rt <- dat$ode$daily_Rt
fig5d <- ggplot(df_stoc_02[81:201,], aes(x = date)) +
geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = "orange", alpha = 0.5) +
geom_ribbon(aes(ymin = `25%`, ymax = `75%`), fill = "orange", alpha = 0.8) +
geom_line(aes(y = `50%`), color = "brown", size = 1.2) +
geom_line(aes(x = date + 1, y = Rt), size = 1, linetype = "dashed") +
geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+
labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") +
coord_cartesian(ylim = c(0.25, 2.5))+
scale_y_continuous(breaks=seq(0, 2.5, by=0.5))+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
# fig5d
# Generate figure 5
df_stoc <- data.frame(date=dat$ode$time, ode=dat$ode$daily_confirmed, c2=dat$stoch_perfect_obs$daily_confirmed[,2], c4=dat$stoch_perfect_obs$daily_confirmed[,4], c19=dat$stoch_perfect_obs$daily_confirmed[,19])
fig5a <- ggplot(df_stoc) + geom_line(aes(x=date, y=ode), linewidth=1, color="black")+geom_line(aes(x=date,y=c2),linewidth=1,color="orange")+geom_line(aes(x=date,y=c4),color="green")+geom_line(aes(x=date,y=c19),color="red") +
labs(y = "Number of confirmation", x = "Time (day)") +
scale_color_hue(labels=c("ode" = "ODE", "c2"="Low", "c4"="Medium", "c19"="High")) +
geom_segment(aes(x=0, xend=5, y=1000, yend=1000),
size=1.1, color = "red", linetype = "solid") +
geom_segment(aes(x=0, xend=5, y=950, yend=950),
size=1.1, color = "green", linetype = "solid") +
geom_segment(aes(x=0, xend=5, y=900, yend=900),
size=1.1, color = "orange", linetype = "solid") +
geom_text(aes(x=7, y=1000), size=4.6, hjust=0, vjust=0.5, label="High intensity") +
geom_text(aes(x=7, y=950), size=4.6, hjust=0, vjust=0.5, label="Medium intensity") +
geom_text(aes(x=7, y=900), size=4.6, hjust=0, vjust=0.5, label="Low intensity")+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
# fig5a
left_col <- plot_grid(fig5a, labels = c("A"))
right_col <- plot_grid(fig5b, fig5c, fig5d, nrow = 3, labels = c("B", "C", "D"))
figure5 <- plot_grid(left_col, right_col, ncol = 2)
figure5
# ggsave("plots/Rt_stoch_recovered_20230619.png", plot=figure4,
# width=3.4*2, height=2.7*2, units="in")Computing mean and std of generation time
# infection (disease) age for generation time
norm_factor <- integrate(function(x) {
1 - pgamma(x, shape = 2, rate = 1/2.5)}, 0, Inf)
f <- function(x) {
(1 - pgamma(x, shape = 2, rate = 1/2.5))/norm_factor$value}
# latent time
g <- function(x) { dexp(x, rate = 1/2.5) }
z <- function(z) {
integrate(function(x, z) g(z-x)*f(x), 0, Inf, z)$value}
fz <- Vectorize(z)
## upper is set to arbitrarily large number
mean_int <- integrate(function(y) y*fz(y), lower = 0,
upper = 100, subdivisions = 1e4)
var_int <- integrate(function(y) y^2*fz(y), lower = 0,
upper = 100, subdivisions = 1e4)
var <- var_int$value - mean_int$value^2 # variance = E[x^2] - E[X]^2
mean_int$value
sqrt(var)
# from Mathematica
mean <- 6.25 # mean_si
std <- 4.14578 # std_siExpect discrepancy (i.e., delay) between the pre-defined Rt and the estimated Rt) #### Figure 5
# Moving average
MA <- function(input_ts, lookback) {
sum <- 0
look_back1 <- lookback + 1
ma_input <- rep(0, length(input_ts))
for (i in 1:length(input_ts)) {
sum = sum + input_ts[i]
if (i < look_back1) {
ma_input[i] <- sum / i
} else {
sum = sum - input_ts[i-7]
ma_input[i] <- sum / 7
}
}
return(ma_input)
}
# apply EpiEstim on confirmation time series (ODE version)
ma_input <- MA(dat$ode$daily_confirmed, 7)
T <- length(ma_input)
time_window <- 1
t_start <- seq(3, T-(time_window-1))
t_end <- t_start + time_window-1
res_daily <- estimate_R(ma_input, method = "parametric_si",
config = make_config((list(t_start = t_start, t_end = t_end, mean_si = 6.25, std_si = 4.14578))))
df_epi_ode_conf <- data.frame(sapply(res_daily$R,c))
df_epi_ode_conf$Rt <- dat$ode$daily_Rt[3:201]
df_epi_ode_conf$time <- t_start-1
col_fill <- "deepskyblue4"
col_dat <- "grey70"
figure6 <-
ggplot(df_epi_ode_conf[80:199,], aes(x = time)) +
geom_ribbon(aes(ymin = Quantile.0.025.R., ymax = Quantile.0.975.R.), fill = "orange", alpha = 0.5) +
geom_ribbon(aes(ymin = Quantile.0.25.R., ymax = Quantile.0.75.R.), fill = "orange", alpha = 0.8) +
geom_line(aes(y = Median.R.), color = "darkblue", size = 1.2) +
geom_line(aes(x = time, y = Rt), colour = "black", size = 1, linetype = "dashed") +
geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+
labs(y = expression(italic(R)[italic(t)]), x = "Time (day)")+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
figure6
# ggsave("plots/EpiEstim_Rt_delay_20230619.png", plot=figure5,
# width=3.4*2, height=2.7, units="in")delay_dist <- {
t <- seq(0, 19)
p_unnorm <- pgamma(t + 1, shape = 3, rate = 1/2.5) -
pgamma(t, shape = 3, rate = 1/2.5)
data.frame(delay = t, pmf = p_unnorm / sum(p_unnorm))
}Estimate Rt by applying deconvolution on confimration series, followed by EpiEstim
# apply on the stochastic seed #19 (fig 6a)
input_ts <- dat$stoch_perfect_obs$daily_confirmed[,19]
# obtain moving average before deconvolution
ma <- seq(0, length(input_ts)-1)
lookback <- 7
sum <- 0
for (i in 1:length(input_ts)){
sum = sum + input_ts[i]
if (i < lookback+1) {
ma[i] = sum / i
} else {
sum = sum - input_ts[i-7]
ma[i] = sum / 7
}
}
result <- deconvolve_rltype_goldsteinetal(
t_obs_min = dat$time[1],
y_obs = round(ma),
delay_min = delay_dist$delay[1],
pmf_delay = delay_dist$pmf,
t_unobs_min = dat$time[1],
n_unobs = length(dat$time),
n_iterations_max = 100,
n_iterations = NULL
)
# output to the above opearation is defined as "result$x_unobs"
unobs <- result$x_unobs
# simple plots for the output of deconvolution (i.e., inferred infection), and the infection from the ODE model
# determine an appropriate lag value by comparing the peaks
# plot(unobs)
# lines(dat$ode$daily_infected)
lag <- which.max(dat$ode$daily_infected) - which.max(unobs)
#lag
# adjust the inferred infection (i.e., unobs)
# exclude first 4 elements of unobs (when, e.g., lag=4) as out of estimation period of our interest
# this causes 4 less estimates (i.e., upto t=196 instead of t=200 for lag=4)
unobs <- unobs[5:length(unobs)]
# now apply EpiEstim
T <- length(unobs)
time_window <- 1
t_start <- seq(3, T-(time_window-1))
t_end <- t_start + time_window-1
res_daily <- estimate_R(unobs, method = "parametric_si",
config = make_config((list(t_start = t_start, t_end = t_end, mean_si = 6.25, std_si = 4.14578))))
# used the mean and std of generation time above for mean_si and std_si
df_epi_19 <- data.frame(sapply(res_daily$R,c))
df_epi_19 <- rbind(df_epi_19, c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA))
df_epi_19$time <- dat$time[3:201]
df_epi_19$Rt <- dat$ode$daily_Rt[3:201]
# df_deconv$sto_per_high <- df_epi_19$Median.R.
plt_epi_19 <-
ggplot(df_epi_19[80:199,], aes(x = time)) +
geom_ribbon(aes(ymin = Quantile.0.025.R., ymax = Quantile.0.975.R.), fill = "#ff6961", alpha = 0.5) +
geom_ribbon(aes(ymin = Quantile.0.25.R., ymax = Quantile.0.75.R.), fill = "red", alpha = 0.8) +
geom_line(aes(y = Median.R.), color = "darkred", size = 1.2) +
geom_line(aes(x = time, y = Rt), colour = "black", size = 1, linetype = "dashed") +
geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+
labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") +
coord_cartesian(ylim = c(0.5, 2))+
scale_y_continuous(breaks=seq(0,2,by=0.5))+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
# apply on the stochastic seed #4 (fig 6b)
input_ts <- dat$stoch_perfect_obs$daily_confirmed[,4]
ma <- seq(0, length(input_ts)-1)
lookback <- 7
sum <- 0
for (i in 1:length(input_ts)){
sum = sum + input_ts[i]
if (i < lookback+1) {
ma[i] = sum / i
} else {
sum = sum - input_ts[i-7]
ma[i] = sum / 7
}
}
result <- deconvolve_rltype_goldsteinetal(
t_obs_min = dat$time[1],
y_obs = round(ma),
delay_min = delay_dist$delay[1],
pmf_delay = delay_dist$pmf,
t_unobs_min = dat$time[1],
n_unobs = length(dat$time),
n_iterations_max = 100,
n_iterations = NULL
)
unobs <- result$x_unobs
unobs <- unobs[5:length(unobs)]
T <- length(unobs)
time_window <- 1
t_start <- seq(3, T-(time_window-1))
t_end <- t_start + time_window-1
res_daily <- estimate_R(unobs, method = "parametric_si",
config = make_config((list(t_start = t_start, t_end = t_end, mean_si = 6.25, std_si = 4.14578))))
# graph
df_epi_4 <- data.frame(sapply(res_daily$R,c))
df_epi_4 <- rbind(df_epi_4, c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA))
df_epi_4$time <- dat$time[3:201]
df_epi_4$Rt <- dat$ode$daily_Rt[3:201]
# df_deconv$sto_per_high <- df_epi_4$Median.R.
plt_epi_4 <-
ggplot(df_epi_4[80:199,], aes(x = time)) +
geom_ribbon(aes(ymin = Quantile.0.025.R., ymax = Quantile.0.975.R.), fill = "#76ff7a", alpha = 0.5) +
geom_ribbon(aes(ymin = Quantile.0.25.R., ymax = Quantile.0.75.R.), fill = "green", alpha = 0.8) +
geom_line(aes(y = Median.R.), color = "darkgreen", size = 1.2) +
geom_line(aes(x = time, y = Rt), colour = "black", size = 1, linetype = "dashed") +
geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+
labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") +
coord_cartesian(ylim = c(0.5, 2))+
scale_y_continuous(breaks=seq(0,2,by=0.5))+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
# apply on the stochastic seed #2 (fig 6c)
input_ts <- dat$stoch_perfect_obs$daily_confirmed[,2]
ma <- seq(0, length(input_ts)-1)
lookback <- 7
sum <- 0
for (i in 1:length(input_ts)){
sum = sum + input_ts[i]
if (i < lookback+1) {
ma[i] = sum / i
} else {
sum = sum - input_ts[i-7]
ma[i] = sum / 7
}
}
result <- deconvolve_rltype_goldsteinetal(
t_obs_min = dat$time[1],
y_obs = round(ma),
delay_min = delay_dist$delay[1],
pmf_delay = delay_dist$pmf,
t_unobs_min = dat$time[1],
n_unobs = length(dat$time),
n_iterations_max = 100,
n_iterations = NULL
)
unobs <- result$x_unobs
unobs <- unobs[5:length(unobs)]
T <- length(unobs)
time_window <- 1
t_start <- seq(3, T-(time_window-1))
t_end <- t_start + time_window-1
res_daily <- estimate_R(unobs, method = "parametric_si",
config = make_config((list(t_start = t_start, t_end = t_end, mean_si = 6.25, std_si = 4.14578))))
df_epi_2 <- data.frame(sapply(res_daily$R,c))
df_epi_2 <- rbind(df_epi_2, c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA))
df_epi_2$time <- dat$time[3:201]
df_epi_2$Rt <- dat$ode$daily_Rt[3:201]
# df_deconv$sto_per_high <- df_epi_2$Median.R.
plt_epi_2 <-
ggplot(df_epi_2[80:199,], aes(x = time)) +
geom_ribbon(aes(ymin = Quantile.0.025.R., ymax = Quantile.0.975.R.), fill = "orange", alpha = 0.5) +
geom_ribbon(aes(ymin = Quantile.0.25.R., ymax = Quantile.0.75.R.), fill = "orange", alpha = 0.8) +
geom_line(aes(y = Median.R.), color = "brown", size = 1) +
geom_line(aes(x = time, y = Rt), colour = "black", size = 1, linetype = "dashed") +
geom_hline(yintercept = 1, color = "darkblue", size = 1, linetype = "dotted")+
labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") +
coord_cartesian(ylim = c(0.5, 2))+
scale_y_continuous(breaks=seq(0,2,by=0.5))+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
figure7 <- plot_grid(plt_epi_19, plt_epi_4, plt_epi_2, nrow = 3, labels = c("A", "B", "C"))
figure7
# ggsave("plots/EpiEstim_Rt_deconvolution_20230619.png", plot=figure6,
# width=3.4*2, height=2.7*2, units="in")Estimate Rt by applying deconvolution on confimration series, followed by EpiEstim
# apply on the stochastic seed #19 (fig 7a)
input_ts <- dat$stoch_perfect_obs$daily_confirmed[,19]
ma <- seq(0, length(input_ts)-1)
lookback <- 7
sum <- 0
for (i in 1:length(input_ts)){
sum = sum + input_ts[i]
if (i < lookback+1) {
ma[i] = sum / i
} else {
sum = sum - input_ts[i-7]
ma[i] = sum / 7
}
}
result <- deconvolve_rltype_goldsteinetal(
t_obs_min = dat$time[1],
y_obs = round(ma),
delay_min = delay_dist$delay[1],
pmf_delay = delay_dist$pmf,
t_unobs_min = dat$time[1],
n_unobs = length(dat$time),
n_iterations_max = 100,
n_iterations = NULL
)
unobs <- result$x_unobs
unobs <- unobs[5:length(unobs)]
T <- length(unobs)
time_window <- 1
t_start <- seq(3, T-(time_window-1))
t_end <- t_start + time_window-1
res_daily <- estimate_R(unobs, method = "parametric_si",
config = make_config((list(t_start = t_start, t_end = t_end, mean_si = 5, std_si = 4.14578))))
# misspecification effect: used mean_si=5, istead of 6.25
df_epi_19m <- data.frame(sapply(res_daily$R,c))
df_epi_19m <- rbind(df_epi_19m, c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA))
df_epi_19m$time <- dat$time[3:201]
df_epi_19m$Rt <- dat$ode$daily_Rt[3:201]
# df_deconv$mis_sto_per_low <- df_epi_2m$Median.R.
plt_epi_19m <-
ggplot(df_epi_19m[80:199,], aes(x = time)) +
geom_ribbon(aes(ymin = Quantile.0.025.R., ymax = Quantile.0.975.R.), fill = "#ff6961", alpha = 0.5) +
geom_ribbon(aes(ymin = Quantile.0.25.R., ymax = Quantile.0.75.R.), fill = "red", alpha = 0.8) +
geom_line(aes(y = Median.R.), color = "darkred", size = 1.2) +
geom_line(aes(x = time, y = Rt), colour = "black", size = 1, linetype = "dashed") +
geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+
labs(y = expression(italic(R)[italic(t)]), x = "Time (day)")+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
# apply on the stochastic seed #4 (fig 7b)
input_ts <- dat$stoch_perfect_obs$daily_confirmed[,4]
ma <- seq(0, length(input_ts)-1)
lookback <- 7
sum <- 0
for (i in 1:length(input_ts)){
sum = sum + input_ts[i]
if (i < lookback+1) {
ma[i] = sum / i
} else {
sum = sum - input_ts[i-7]
ma[i] = sum / 7
}
}
result <- deconvolve_rltype_goldsteinetal(
t_obs_min = dat$time[1],
y_obs = round(ma),
delay_min = delay_dist$delay[1],
pmf_delay = delay_dist$pmf,
t_unobs_min = dat$time[1],
n_unobs = length(dat$time),
n_iterations_max = 100,
n_iterations = NULL
)
# output to the above opearation is defined as "result$x_unobs"
unobs <- result$x_unobs
unobs <- unobs[5:length(unobs)]
# epiestim_4m
T <- length(unobs)
time_window <- 1
t_start <- seq(3, T-(time_window-1))
t_end <- t_start + time_window-1
res_daily <-
estimate_R(unobs, method = "parametric_si",
config = make_config((list(t_start = t_start, t_end = t_end, mean_si = 5, std_si = 4.14578))))
df_epi_4m <- data.frame(sapply(res_daily$R,c))
df_epi_4m <- rbind(df_epi_4m, c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA))
df_epi_4m$time <- dat$time[3:201]
df_epi_4m$Rt <- dat$ode$daily_Rt[3:201]
# df_deconv$mis_sto_per_low <- df_epi_2m$Median.R.
plt_epi_4m <-
ggplot(df_epi_4m[80:199,], aes(x = time)) +
geom_ribbon(aes(ymin = Quantile.0.025.R., ymax = Quantile.0.975.R.), fill = "#76ff7a", alpha = 0.5) +
geom_ribbon(aes(ymin = Quantile.0.25.R., ymax = Quantile.0.75.R.), fill = "green", alpha = 0.8) +
geom_line(aes(y = Median.R.), color = "darkgreen", size = 1.2) +
geom_line(aes(x = time, y = Rt), colour = "black", size = 1, linetype = "dashed") +
geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+
labs(y = expression(italic(R)[italic(t)]), x = "Time (day)")+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
# apply on the stochastic seed #2 (fig 7c)
input_ts <- dat$stoch_perfect_obs$daily_confirmed[,2]
ma <- seq(0, length(input_ts)-1)
sum <- 0
for (i in 1:length(input_ts)){
sum = sum + input_ts[i]
if (i < lookback+1) {
ma[i] = sum / i
} else {
sum = sum - input_ts[i-7]
ma[i] = sum / 7
}
}
result <- deconvolve_rltype_goldsteinetal(
t_obs_min = dat$time[1],
y_obs = round(ma),
delay_min = delay_dist$delay[1],
pmf_delay = delay_dist$pmf,
t_unobs_min = dat$time[1],
n_unobs = length(dat$time),
n_iterations_max = 100,
n_iterations = NULL
)
unobs <- result$x_unobs
unobs <- unobs[5:length(unobs)]
T <- length(unobs)
time_window <- 1
t_start <- seq(3, T-(time_window-1))
t_end <- t_start + time_window-1
res_daily <- estimate_R(unobs, method = "parametric_si",
config = make_config((list(t_start = t_start, t_end = t_end, mean_si = 5, std_si = 4.14578))))
df_epi_2m <- data.frame(sapply(res_daily$R,c))
df_epi_2m <- rbind(df_epi_2m, c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA))
df_epi_2m$time <- dat$time[3:201]
df_epi_2m$Rt <- dat$ode$daily_Rt[3:201]
# df_deconv$mis_sto_per_low <- df_epi_2m$Median.R.
plt_epi_2m <-
ggplot(df_epi_2m[80:199,], aes(x = time)) +
geom_ribbon(aes(ymin = Quantile.0.025.R., ymax = Quantile.0.975.R.), fill = "orange", alpha = 0.5) +
geom_ribbon(aes(ymin = Quantile.0.25.R., ymax = Quantile.0.75.R.), fill = "orange", alpha = 0.8) +
geom_line(aes(y = Median.R.), color = "brown", size = 1) +
geom_line(aes(x = time, y = Rt), colour = "black", size = 1, linetype = "dashed") +
geom_hline(yintercept = 1, color = "darkblue", size = 1, linetype = "dotted")+
labs(y = expression(italic(R)[italic(t)]), x = "Time (day)")+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
figure8 <- plot_grid(plt_epi_19m, plt_epi_4m, plt_epi_2m, nrow = 3, labels = c("a", "b", "c"))
figure8
# ggsave("plots/EpiEstim_stoch_Rt_deconvolution_20230619.png", plot=figure7,
# width=3.4*2, height=2.7*2, units="in")Plot for the curve of pre-defined Rt by periods of different epidemic characteristics
R_flat <- c(dat$Rt[1:101], rep(NA,32), dat$Rt[134:195], rep(NA, 6))
R_curve <- c(rep(NA, 100), dat$Rt[101:127], rep(NA, 74))
R_curve_last <- c(rep(NA, 126), dat$Rt[127:134], rep(NA, 67))
R_last <- c(rep(NA, 194), dat$Rt[195:201])
figure9 <- ggplot()+
geom_line(aes(x=dat$time, y=R_flat), color="black", linewidth=1.2)+
geom_line(aes(x=dat$time, y=R_curve), color="red", linewidth=1.2)+
geom_line(aes(x=dat$time, y=R_curve_last), linetype="dotted", color="red", linewidth=1.2)+
geom_line(aes(x=dat$time, y=R_last), linetype="dotted", linewidth=1.2)+
labs(y=expression(italic(R)[italic(t)]), x="Time (day)")+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
figure9
# ggsave("plots/Rt_sections_20230619.png", plot=figure8,
# width=3.4*2, height=2.7*2, units="in")input_ts <- dat$stoch_perfect_obs$daily_confirmed[,4][1:133]
# PF
d_stoc_04_mid <- data.frame(date = dat$time[1:133],
daily_confirmed = round(dat$stoch_perfect_obs$daily_confirmed[,4][1:133]))set.seed(42)
cl <- makeCluster(getOption("cl.cores", 2))
doParallel::registerDoParallel(cl)
pf_stoc_04_mid <- foreach(i = 1:nrep, .packages = c("pfilterCOVID"), .inorder = F) %dopar% {
extract_trace(params = params,
y = params$Y0,
data = d_stoc_04_mid,
data_type = "confirmation",
npart = npart,
tend = nrow(d_stoc_04_mid),
dt = dt,
error_pdf = "pois",
negbin_size = 15,
backward_sampling = TRUE,
stoch = FALSE)
}
parallel::stopCluster(cl)
# saveRDS(pf_stoc_04_mid, "outputs/pf_stoc_04_mid_20230619.rds") pf_stoc_04_mid <- readRDS("outputs/pf_stoc_04_mid_20230619.rds")
parset_chr <- paste0("confirmation", "_I0=", params$I0, "_npart=", npart,
"_nrep=", nrep, "_dt=", dt, "_", tstamp)
pr <- c(0.025, 0.25, 0.5, 0.75, 0.975)
Rt_est <- as.data.frame(sapply(pf_stoc_04_mid, function(x) x[, "Rt"]))
Rt_quantile <- as.data.frame(t(apply(Rt_est, 1, function(x) quantile(x, pr))))
df_stoc_04_mid <- cbind(Rt_quantile, d_stoc_04_mid)
col_fill <- "deepskyblue4"
col_dat <- "grey70"
df_stoc_04_mid$Rt <- dat$ode$daily_Rt[1:133]
fig10a <- ggplot(df_stoc_04_mid[81:133,], aes(x = date)) +
geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill="orange", alpha = 0.5) +
geom_ribbon(aes(ymin = `25%`, ymax = `75%`), fill="orange", alpha = 0.8) +
geom_line(aes(y = `50%`), color="brown", size = 1.2) +
geom_line(aes(x = date, y = df_stoc_04_mid$Rt[80:132]), size = 1, linetype = "dashed") +
geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+
labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") +
coord_cartesian(ylim = c(0.25, 2.75))+scale_y_continuous(breaks=seq(0,2.75,by=0.5))+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
# EpiEstim
ma <- seq(0, length(input_ts)-1)
lookback <- 7
sum <- 0
for (i in 1:length(input_ts)){
sum = sum + input_ts[i]
if (i < lookback+1) {
ma[i] = sum / i
} else {
sum = sum - input_ts[i-7]
ma[i] = sum / 7
}
}
result <- deconvolve_rltype_goldsteinetal(
t_obs_min = dat$time[1],
y_obs = round(ma),
delay_min = delay_dist$delay[1],
pmf_delay = delay_dist$pmf,
t_unobs_min = dat$time[1],
n_unobs = length(dat$time[1:133]),
n_iterations_max = 100,
n_iterations = NULL
)
unobs <- result$x_unobs
# adjustment of unobs
unobs <- unobs[5:133]
T <- length(unobs)
time_window <- 1
t_start <- seq(3, T-(time_window-1))
t_end <- t_start + time_window-1
# res_daily: without misspecification of SI
# res_daily_m: with misspecification of SI
res_daily <- estimate_R(unobs, method = "parametric_si",
config = make_config((list(t_start = t_start, t_end = t_end, mean_si = 6.25, std_si = 4.14578))))
res_daily_m <- estimate_R(unobs, method = "parametric_si",
config = make_config((list(t_start = t_start, t_end = t_end, mean_si = 5, std_si = 4.14578))))
df_epi_4_mid <- data.frame(sapply(res_daily$R, c))
df_epi_4_mid_m <- data.frame(sapply(res_daily_m$R, c))
df_epi_4_mid <- df_epi_4_mid[,c(5,7,8,9,11)]
df_epi_4_mid_m <- df_epi_4_mid_m[,c(5,7,8,9,11)]
df_epi_4_mid <- rbind(df_epi_4_mid, c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA))
df_epi_4_mid_m <- rbind(df_epi_4_mid_m, c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA), c(NA,NA,NA,NA,NA))
df_epi_4_mid$time <- dat$time[3:133]
df_epi_4_mid$Rt <- dat$ode$daily_Rt[3:133]
df_epi_4_mid_m$time <- dat$time[3:133]
df_epi_4_mid_m$Rt <- dat$ode$daily_Rt[3:133]
# fig10b
fig10b <-
ggplot(df_epi_4_mid[79:131,], aes(x = time)) +
geom_ribbon(aes(ymin = Quantile.0.025.R., ymax = Quantile.0.975.R.), fill = col_fill, alpha = 0.5) +
geom_ribbon(aes(ymin = Quantile.0.25.R., ymax = Quantile.0.75.R.), fill = col_fill, alpha = 0.8) +
geom_line(aes(y = Median.R.), color = "darkblue", size = 1.2) +
geom_line(aes(y = Rt), colour = "black", size = 1, linetype = "dashed") +
geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+
labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") +
coord_cartesian(ylim = c(0.25, 2.75))+
scale_y_continuous(breaks=seq(0,2.75,by=0.5))+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
# fig10c
fig10c <-
ggplot(df_epi_4_mid_m[79:131,], aes(x = time)) +
geom_ribbon(aes(ymin = Quantile.0.025.R., ymax = Quantile.0.975.R.),
fill="orange", alpha = 0.5) +
geom_ribbon(aes(ymin = Quantile.0.25.R., ymax=Quantile.0.75.R.),
fill="orange", alpha = 0.8) +
geom_line(aes(y = Median.R.), color = "brown", size = 1.2) +
geom_line(aes(y = Rt), colour = "black", size = 1, linetype = "dashed") +
geom_hline(yintercept = 1, color = "black", size = 1, linetype = "dotted")+
labs(y = expression(italic(R)[italic(t)]), x = "Time (day)") +
coord_cartesian(ylim = c(0.25, 2.75))+
scale_y_continuous(breaks=seq(0,2.75,by=0.5))+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13))
figure10 <- plot_grid(fig10a, fig10b, fig10c, nrow=3, labels=c("A", "B", "C"))
figure10
# ggsave("plots/Rt_PF_EpiEstim_20230619.png", plot=figure9,
# width=3.4*2, height=2.7*2, units="in")dg <- df_stoc_04[82:201,]
dg1 <- df_epi_4[84:199,]
dg2 <- df_epi_4m[84:199,]
# #
flat_Rt <- c(dg$Rt[1:20], dg$Rt[53:120])
flat_PF_Rt <- c(dg$`50%`[1:20], dg$`50%`[53:120])
flat_epi_Rt <- c(dg1$Median.R.[1:20], dg1$Median.R.[53:116], NA, NA, NA, NA)
flat_epim_Rt <- c(dg2$Median.R.[1:20], dg2$Median.R.[53:116], NA, NA, NA, NA)
# #
general_period <- list("PF" = sqrt(mean((df_stoc_04[82:201,]$Rt - df_stoc_04[82:201,]$`50%`)^2, na.rm=TRUE)),
"EpiEstim" = sqrt(mean((df_stoc_04[82:197,]$Rt - df_epi_4[84:199,]$Median.R.)^2, na.rm=TRUE)),
"EpiEstim_m" = sqrt(mean((df_stoc_04[82:197,]$Rt - df_epi_4m[84:199,]$Median.R.)^2, na.rm=TRUE)))
flat_period <- list("PF" = sqrt(mean((flat_Rt - flat_PF_Rt)^2, na.rm=TRUE)),
"EpiEstim" = sqrt(mean((flat_Rt - flat_epi_Rt)^2, na.rm=TRUE)),
"EpiEstim_m" = sqrt(mean((flat_Rt - flat_epim_Rt)^2, na.rm=TRUE)))
last_flat <- list("PF" = sqrt(mean((flat_Rt[82:88] - flat_PF_Rt[82:88])^2, na.rm=TRUE)),
"EpiEstim" = sqrt(mean((flat_Rt[82:88] - flat_epi_Rt[82:88])^2, na.rm=TRUE)),
"EpiEstim_m" = sqrt(mean((flat_Rt[82:88] - flat_epim_Rt[82:88])^2, na.rm=TRUE)))
fluctuation <- list("PF" = sqrt(mean((dat$ode$daily_Rt[99:133] - df_stoc_04_mid[99:133,]$`50%`)^2, na.rm=TRUE)),
"EpiEstim" = sqrt(mean((dat$ode$daily_Rt[99:133] - df_epi_4_mid[99:133,]$Median.R.)^2, na.rm=TRUE)),
"EpiEstim_m" = sqrt(mean((dat$ode$daily_Rt[99:133] - df_epi_4_mid_m[99:133,]$Median.R.)^2, na.rm=TRUE)))
last_fluctuaion <- list("PF" = sqrt(mean((dat$ode$daily_Rt[127:133] - df_stoc_04_mid[127:133,]$`50%`)^2, na.rm=TRUE)),
"EpiEstim" = sqrt(mean((dat$ode$daily_Rt[127:133] - df_epi_4_mid[127:133,]$Median.R.)^2, na.rm=TRUE)),
"EpiEstim_m" = sqrt(mean((dat$ode$daily_Rt[127:133] - df_epi_4_mid_m[127:133,]$Median.R.)^2, na.rm=TRUE)))
# RMSE table (in the order of PF, EpiEstim,and EpiEstim with misspecification)
table2 <- list("general period" = general_period, "flat period" = flat_period, "last 1 week flat" = last_flat, "fluctuation period" = fluctuation, "last 1 week fluctuation" = last_fluctuaion)set.seed(42)
cl <- makeCluster(getOption("cl.cores", detectCores()))
registerDoParallel(cl)
nrep <- 1000
npart <- 10000
params[["betavol"]] <- 0.20
y0 <- params$Y0
y0["I"] <- 100
y0["E"] <- 100
data <- read_csv("data/daily_confirmed_202009_202110.csv")
# data <- data[order(data$date),]
# write_csv(data, "data/daily_confirmed_202009_202110.csv")
# data$date <- 0:(nrow(data)-1)
pf_conf <- foreach(i=1:nrep, .packages=c("pfilterCOVID"), .inorder=F) %dopar% {
extract_trace(params = params,
y = y0,
data = data,
data_type = "confirmation",
npart = npart,
tend = nrow(data),
dt = 0.1,
error_pdf = "negbin",
negbin_size = 100,
backward_sampling = TRUE,
stoch = FALSE)
}
parallel::stopCluster(cl)
# saveRDS(pf_conf, "outputs/pf_conf_covid_korea_20230729.rds")Generate Figure 11 (supplmentary figure) Red and green vertical lines represent strengthening and weakening of the social distancing measures. ### Figure 11
dtype <- "confirmation"
tstamp <- format(Sys.Date(), "%Y%m%d")
data <- read_csv("data/daily_confirmed_202009_202110.csv")
pf_korea <- readRDS("outputs/pf_conf_covid_korea_20230729.rds")
parset_chr <- paste0(dtype, "_I0=", params$I0, "_npart=", npart,
"_nrep=", nrep, "_dt=", dt, "_", tstamp)
pr <- c(0.025, 0.25, 0.5, 0.75, 0.975)
Rt_est <- as.data.frame(sapply(pf_korea, function(x) x[, "Rt"]))
Rt_quantile <- as.data.frame(t(apply(Rt_est, 1, function(x) quantile(x, pr))))
cumul_est <- as.data.frame(sapply(pf_korea, function(x) x[, "CR"]))
cumul_quantile <- as.data.frame(t(apply(cumul_est, 1, function(x) quantile(x, pr))))
Rtdf <- cbind(Rt_quantile, data)
Rtdf$date <- data$date
cumuldf <- cbind(cumul_quantile, data)
cumuldf$date <- data$date
npi_data <- data.frame(date=as.Date(c("2020-10-11", "2021-03-12", "2021-06-11",
"2020-11-17", "2020-11-22", "2020-12-05", "2020-12-08", "2020-12-21")),
degree = c("Low","Low","Low","High","High","High","High","High"))
npi_high <- npi_data %>% filter(degree=="High")
npi_low <- npi_data %>% filter(degree=="Low")
# Rt
fig11a <- ggplot(Rtdf, aes(x = date)) +
geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = "steelblue", alpha = 0.5) +
geom_ribbon(aes(ymin = `25%`, ymax = `75%`), fill = "steelblue", alpha = 0.8) +
geom_line(aes(y = `50%`), color = "steelblue", size = 1.2) +
geom_hline(yintercept = 1, color = "black", size = 1, linetype="dotted")+
geom_vline(xintercept=npi_high$date, color="firebrick", size=1)+
geom_vline(xintercept=npi_low$date, color="darkgreen", size=1)+
labs(y = expression(italic(R)[italic(t)]), x="") +
scale_x_date(limits=as.Date(c("2020-10-01","2021-06-30")), date_breaks = "2 week")+
scale_y_continuous(limits=c(0,3),breaks=seq(0,3,by=0.5))+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13),
axis.text.x=element_text(angle=90,hjust=1))
# fig11a
fig11b <- ggplot(cumuldf) +
geom_col(aes(x=date, y=daily_confirmed), fill="grey50", alpha=0.3) +
geom_ribbon(aes(x=date, ymin=`2.5%`, ymax=`97.5%`),
fill="steelblue", alpha=0.3) +
geom_ribbon(aes(x = date, ymin = `25%`, ymax=`75%`),
fill="steelblue", alpha=0.7) +
geom_line(aes(x=date, y=`50%`), color="steelblue", size=1)+
# scale_y_continuous(limits=c(0,2000),breaks=seq(0,1000,by=200))+
labs(title = "", x = "", y = "Daily confirmed case") +
scale_x_date(limits=as.Date(c("2020-10-01","2021-06-30")), date_breaks = "2 week")+
theme(text = element_text(size=16),
axis.text = element_text(size=13),
legend.text=element_text(size=13),
axis.text.x=element_text(angle=90,hjust=1))
# fig11b
figure11 <- plot_grid(fig11a, fig11b, nrow=2, labels="AUTO")
figure11
# ggsave("plots/Rt_confimed_case_korea_20230730.png", plot=figure11,
# width=3.4*2, height=2.7*3, units="in")