Nothing
R/tvt.eSIR.R
In eSIR: Extended State-Space SIR Models
Defines functions
tvt.eSIR
Documented in
tvt.eSIR
# This is the source code for R package eSIR: extended-SIR
# Built on Feb 13, 2020, and last edited on Oct 27, 2021
# Correspondence : Peter X.K. Song, Ph.D. (pxsong@umich.edu)
# Creator: Lili Wang, Ph.D. (lilywang@umich.edu)
# Model 1: An extended SIR model with a time-varying transmission rate modifier
# pi(t)
#' Fit extended state-space SIR model with time-varying transmission rates
#'
#' Fit extended state-space SIR model with pre-specified changes in the transmission rate, either stepwise or continuous, accommodating time-varying quarantine protocols.
#'
#' We fit a state-space model with extended SIR, in which a time-varying transmission rate modifier \eqn{\pi(t)} (between 0 and 1) is introduced to model. This allows us to accommodate quarantine protocol changes and ignored resources of hospitalization. The form of reducing rate may be a step-function with jumps at times of big policy changes or a smooth exponential survival function \eqn{\exp(-\lambda_0t)}. The parameters of the function and change points, if any, should be predefined.
#'
#'
#' @param Y the time series of daily observed infected compartment proportions.
#' @param R the time series of daily observed removed compartment proportions, including death and recovered.
#' @param pi0 the time-dependent transmission rate modifier \eqn{\pi(t)} between 0 and 1.
#' @param change_time the change points over time for step function pi, defalt value is \code{NULL}.
#' @param exponential logical, whether \eqn{\pi(t)} is exponential \eqn{\exp(-\lambda_0t)} or not; the default is \code{FALSE}.
#' @param lambda0 the rate of decline in the exponential survival function in \eqn{\exp(-\lambda_0t)}.
#' @param begin_str the character of starting time, the default is "01/13/2020".
#' @param T_fin the end of follow-up time after the beginning date \code{begin_str}, the default is 200.
#' @param nchain the number of MCMC chains generated by \code{\link[rjags]{rjags}}, the default is 4.
#' @param nadapt the iteration number of adaptation in the MCMC. We recommend using at least the default value 1e4 to obtained fully adapted chains.
#' @param M the number of draws in each chain, with no thinning. The default is M=5e2 but suggest using 5e5.
#' @param thn the thinning interval between mixing. The total number of draws thus would become \code{round(M/thn)*nchain}. The default is 10.
#' @param nburnin the burn-in period. The default is 2e2 but suggest 2e5.
#' @param dic logical, whether compute the DIC (deviance information criterion) for model selection.
#' @param death_in_R the numeric value of average of cumulative deaths in the removed compartments. The default is 0.4 within Hubei and 0.02 outside Hubei.
#' @param casename the string of the job's name. The default is "tvt.eSIR".
#' @param beta0 the hyperparameter of average transmission rate, the default is the one estimated from the SARS first-month outbreak (0.2586).
#' @param gamma0 the hyperparameter of average removed rate, the default is the one estimated from the SARS first-month outbreak (0.0821).
#' @param R0 the hyperparameter of the mean reproduction number R0. The default is thus the ratio of \code{beta0/gamma0}, which can be specified directly.
#' @param gamma0_sd the standard deviation for the prior distrbution of the removed rate \eqn{\gamma}, the default is 0.1.
#' @param R0_sd the standard deviation for the prior disbution of R0, the default is 1.
#' @param file_add the string to denote the location of saving output files and tables.
#' @param save_files logical, whether save (\code{TRUE}) results or not (\code{FALSE}). This enables to save summary tables, trace plots, and plots of the posterior means of the first-order derivatives of the infection prevalence process \eqn{\theta_t^I}.
#' @param save_mcmc logical, whether save (\code{TRUE}) all the MCMC outputs or not (\code{FALSE}).The output file will be an \code{.RData} file named by the \eqn{casename}. We include arrays of prevalence values of the three compartments with their matrices of posterior draws up to the last date of the collected data as \code{theta_p[,,1]} and afterwards as \code{theta_pp[,,1]} for \eqn{\theta_t^S}, \code{theta_p[,,2]} and \code{theta_pp[,,2]} for \eqn{\theta_t^I}, and \code{theta_p[,,3]} and \code{theta_pp[,,3]} for \eqn{\theta_t^R}. Moreover, the input and predicted proportions \code{Y}, \code{Y_pp}, \code{R} and \code{R_pp} can also be retrieved. The prevalence and predicted proportion matrices have rows for MCMC replicates, and columns for days. The MCMC posterior draws of other parameters including \code{beta_p}, \code{gamma_p}, \code{R0_p}, and variance controllers \code{k_p}, \code{lambdaY_p}, \code{lambdaR_p} are also available.
#' @param save_plot_data logical, whether save the plotting data or not.
#' @param save_files logical, whether to save plots to file.
#' @param add_death logical, whether add the approximate death curve to the plot, default is false.
#' @param eps a non-zero controller so that all the input \code{Y} and \code{R} values would be bounded above 0 (at least \code{eps}). Its default value is 1e-10
#'
#' @return
#' \item{casename}{the predefined \code{casename}.}
#' \item{incidence_mean}{mean cumulative incidence, the mean prevalence of cumulative confirmed cases at the end of the study.}
#' \item{incidence_ci}{2.5\%, 50\%, and 97.5\% quantiles of the incidences.}
#' \item{out_table}{summary tables including the posterior mean of the prevalance processes of the 3 states compartments (\eqn{\theta_t^S,\theta_t^I,\theta_t^R}) at last date of data collected ((\eqn{t^\prime}) decided by the lengths of your input data \code{Y} and \code{R}), and their respective credible inctervals (ci); the respective means and ci's of the reporduction number (R0), removed rate (\eqn{\gamma}), transmission rate (\eqn{\beta}).}
#' \item{plot_infection}{plot of summarizing and forecasting for the infection compartment, in which the vertial blue line denotes the last date of data collected (\eqn{t^\prime}), the vertial darkgray line denotes the deacceleration point (first turning point) that the posterior mean first-derivative of infection prevalence \eqn{\dot{\theta}_t^I} achieves the maximum, the vertical purple line denotes the second turning point that the posterior mean first-derivative infection proportion \eqn{\dot{\theta}_t^I} equals zero, the darkgray line denotes the posterior mean of the infection prevalence \eqn{\theta_t^I} and the red line denotes its posterior median. }
#' \item{plot_removed}{plot of summarizing and forecasting for the removed compartment with lines similar to those in the \code{plot_infection}. The vertical lines are identical, but the horizontal mean and median correspond to the posterior mean and median of the removed process \eqn{\theta_t^R}. An additional line indicates the estimated death prevalence from the input \code{death_in_R}.}
#' \item{spaghetti_plot}{20 randomly selected MCMC draws of the first-order derivative of the posterior prevalence of infection, namely \eqn{\dot{\theta}_t^I}. The black curve is the posterior mean of the derivative, and the vertical lines mark times of turning points corresponding respectively to those shown in \code{plot_infection} and \code{plot_removed}. Moreover, the 95\% credible intervals of these turning points are also highlighted by semi-transparent rectangles. }
#' \item{first_tp_mean}{the date t at which \eqn{\ddot{\theta}_t^I=0}, calculated as the average of the time points with maximum posterior first-order derivatives \eqn{\dot{\theta}_t^I}; this value may be slightly different from the one labeled by the "darkgreen" lines in the two plots \code{plot_infection} and \code{plot_removed}, which indicate the stationary point such that the first-order derivative of the averaged posterior of \eqn{\theta_t^I} reaches its maximum.}
#'
#' \item{first_tp_ci}{fwith \code{first_tp_mean}, it reports the corresponding credible interval and median.}
#' \item{second_tp_mean}{the date t at which \eqn{\theta_t^I=0}, calculated as the average of the stationary points of all of posterior first-order derivatives \eqn{\dot{\theta}_t^I}; this value may be slightly different from the one labeled by the "pruple" lines in the plots of \code{plot_infection} and \code{plot_removed}. The latter indicate stationary t at which the first-order derivative of the averaged posterior of \eqn{\theta_t^I} equals zero.}
#' \item{second_tp_ci}{with \code{second_tp_mean}, it reports the corresponding credible interval and median.}
#' \item{dic_val}{the output of \code{dic.samples()} in \code{\link[rjags]{dic.samples}}, computing deviance information criterion for model comparison.}
#' \item{gelman_diag_list}{ Since version 0.3.3, we incorporated Gelman And Rubin's Convergence Diagnostic using \code{\link[coda]{gelman.diag}}. We included both the statistics and their upper C.I. limits. Values substantially above 1 indicate lack of convergence. Error messages would be printed as they are. This would be only valid for multiple chains (e.g. nchain > 1). Note that for time dependent processes, we only compute the convergence of the last observation data (\code{T_prime}), though it shows to be \code{T_prime+1}, which is due to the day 0 for initialization.}
#' @examples
#' \donttest{
#' NI_complete <- c(
#' 41, 41, 41, 45, 62, 131, 200, 270, 375, 444, 549, 729,
#' 1052, 1423, 2714, 3554, 4903, 5806, 7153, 9074, 11177,
#' 13522, 16678, 19665, 22112, 24953, 27100, 29631, 31728, 33366
#' )
#' RI_complete <- c(
#' 1, 1, 7, 10, 14, 20, 25, 31, 34, 45, 55, 71, 94, 121, 152, 213,
#' 252, 345, 417, 561, 650, 811, 1017, 1261, 1485, 1917, 2260,
#' 2725, 3284, 3754
#' )
#' N <- 58.5e6
#' R <- RI_complete / N
#' Y <- NI_complete / N - R # Jan13->Feb 11
#' ### Step function of pi(t)
#' change_time <- c("01/23/2020", "02/04/2020", "02/08/2020")
#' pi0 <- c(1.0, 0.9, 0.5, 0.1)
#' res.step <- tvt.eSIR(Y, R,
#' begin_str = "01/13/2020", death_in_R = 0.4,
#' T_fin = 200, pi0 = pi0, change_time = change_time, dic = TRUE,
#' casename = "Hubei_step", save_files = TRUE,
#' save_mcmc = FALSE, M = 5e2, nburnin = 2e2
#' )
#' res.step$plot_infection
#' res.step$plot_removed
#' res.step$dic_val
#'
#' ### continuous exponential function of pi(t)
#' res.exp <- tvt.eSIR(Y, R,
#' begin_str = "01/13/2020", death_in_R = 0.4,
#' T_fin = 200, exponential = TRUE, dic = FALSE, lambda0 = 0.05,
#' casename = "Hubei_exp", save_files = FALSE, save_mcmc = FALSE,
#' M = 5e2, nburnin = 2e2
#' )
#' res.exp$plot_infection
#' # res.exp$plot_removed
#'
#' ### without pi(t), the standard state-space SIR model without intervention
#' res.nopi <- tvt.eSIR(Y, R,
#' begin_str = "01/13/2020", death_in_R = 0.4,
#' T_fin = 200, casename = "Hubei_nopi", save_files = FALSE,
#' M = 5e2, nburnin = 2e2
#' )
#' res.nopi$plot_infection
#' # res.nopi$plot_removed
#' }
#'
#' change_time <- c("01/18/2020")
#' pi0<- c(1.0, 0.9)
#' NI_complete2 <- c(41, 45, 62, 131)
#' RI_complete2 <- c(1, 1, 7, 10)
#' N2 <- 1E3
#' res1 <- tvt.eSIR(
#' RI_complete2 / N2,
#' (NI_complete2 - RI_complete2) / N2,
#' begin_str = "01/10/2020",
#' T_fin =10,
#' pi0 = pi0,
#' change_time = change_time,
#' dic = FALSE,
#' casename = "Hubei_step",
#' save_files = FALSE,
#' save_mcmc = FALSE,
#' save_plot_data = FALSE,
#' M = 50,
#' nburnin = 1
#' )
#' closeAllConnections()
#'
#' @export
tvt.eSIR <- function(Y,
R,
pi0 = NULL,
change_time = NULL,
exponential = FALSE,
lambda0 = NULL,
begin_str = "01/13/2020",
T_fin = 200,
nchain = 4,
nadapt = 1e4,
M = 5e2,
thn = 10,
nburnin = 2e2,
dic = FALSE,
death_in_R = 0.02,
beta0 = 0.2586,
gamma0 = 0.0821,
R0 = beta0 / gamma0,
gamma0_sd = 0.1,
R0_sd = 1,
casename = "tvt.eSIR",
file_add = character(0),
add_death = FALSE,
save_files = FALSE,
save_mcmc = FALSE,
save_plot_data = FALSE,
eps = 1e-10) {
oldpar <- par(no.readonly = TRUE)
on.exit(par(oldpar))
beta0 <- R0 * gamma0
len <- round(M / thn) * nchain # number of MCMC draws in total
T_prime <- length(Y)
if (T_prime != length(R)) stop("Y and R should be matched.")
if(T_prime >= T_fin) {
warning("T_fin must be larger than input time series Y and R,",
"we will automatically add 100 time units to it: T_fin = length(Y) + 10.")
T_fin = T_prime + 10
}
Y <- pmax(Y, eps)
R <- pmax(R, eps)
if (add_death == TRUE && death_in_R == 0.02) {
message("use the default death_in_R which is equal to 0.02",
" to plot the death curve in the removed process forecast plot")}
begin <- chron(dates. = begin_str)
chron_ls <- chron(begin:(begin + T_fin))
end <- chron(begin:(begin + T_fin))[T_fin]
message(
"The follow-up is from ", begin, " to ", end,
" and the last observed date is ", chron_ls[T_prime], "."
)
# current data up to this date
gamma_var <- gamma0_sd^2
lognorm_gamma_parm <- lognorm.parm(gamma0, gamma_var)
R0_var <- R0_sd^2
lognorm_R0_parm <- lognorm.parm(R0, R0_var)
if (exponential == FALSE & !is.null(change_time) & !is.null(pi0)) {
message("Running for step-function pi(t)")
if (length(change_time) != length(pi0) - 1) {
stop("We need the length of vector change_time",
" to be the length of pi0 minus 1.")
}
change_time_chorn <- c(begin - 1, chron(dates. = change_time), end)
pi <- rep(pi0, diff(change_time_chorn))
} else if (exponential == TRUE & !is.null(lambda0)) {
message("Running for exponential-function pi(t)")
pi <- exp(-lambda0 * (0:(T_fin - 1)))
} else {
message("Running without pi(t)")
pi <- rep(1, T_fin)
}
################ MCMC ##########
model1.string <- paste0("
model{
for(t in 2:(T_prime+1)){
Km[t-1,1] <- -beta*pi[t-1]*theta[t-1,1]*theta[t-1,2]
Km[t-1,9] <- gamma*theta[t-1,2]
Km[t-1,5] <- -Km[t-1,1]-Km[t-1,9]
Km[t-1,2] <- -beta*pi[t-1]*(theta[t-1,1]+0.5*Km[t-1,1])*(theta[t-1,2]+0.5*Km[t-1,5])
Km[t-1,10] <- gamma*(theta[t-1,2]+0.5*Km[t-1,5])
Km[t-1,6] <- -Km[t-1,2]-Km[t-1,10]
Km[t-1,3] <- -beta*pi[t-1]*(theta[t-1,1]+0.5*Km[t-1,2])*(theta[t-1,2]+0.5*Km[t-1,6])
Km[t-1,11] <- gamma*(theta[t-1,2]+0.5*Km[t-1,6])
Km[t-1,7] <- -Km[t-1,3]-Km[t-1,11]
Km[t-1,4] <- -beta*pi[t-1]*(theta[t-1,1]+Km[t-1,3])*(theta[t-1,2]+Km[t-1,7])
Km[t-1,12] <- gamma*(theta[t-1,2]+Km[t-1,7])
Km[t-1,8] <- -Km[t-1,4]-Km[t-1,12]
alpha[t-1,1] <- theta[t-1,1]+(Km[t-1,1]+2*Km[t-1,2]+2*Km[t-1,3]+Km[t-1,4])/6
alpha[t-1,2] <- theta[t-1,2]+(Km[t-1,5]+2*Km[t-1,6]+2*Km[t-1,7]+Km[t-1,8])/6
alpha[t-1,3] <- theta[t-1,3]+(Km[t-1,9]+2*Km[t-1,10]+2*Km[t-1,11]+Km[t-1,12])/6
theta[t,1:3] ~ ddirch(k*alpha[t-1,1:3])
Y[t-1] ~ dbeta(lambdaY*theta[t,2],lambdaY*(1-theta[t,2]))
R[t-1] ~ dbeta(lambdaR*theta[t,3],lambdaR*(1-theta[t,3]))
}
theta0[1:3]<-c(", 1 - Y[1] - R[1], ",", Y[1], ",", R[1], ")
theta[1,1:3] ~ ddirch(theta0[1:3])
gamma ~ dlnorm(", lognorm_gamma_parm$mu, ",", 1 / lognorm_gamma_parm$var, ")
R0 ~ dlnorm(", lognorm_R0_parm$mu, ",", 1 / lognorm_R0_parm$var, ")
beta <- R0*gamma
k ~ dgamma(2,0.0001)
lambdaY ~ dgamma(2,0.0001)
lambdaR ~ dgamma(2,0.0001)
}
")
model.spec <- textConnection(model1.string)
posterior <- suppressWarnings(
jags.model(
model.spec,
data = list(
"Y" = Y,
"R" = R,
"T_prime" = T_prime,
"pi" = pi
),
n.chains = nchain,
n.adapt = nadapt
)
)
suppressWarnings(
update(posterior, nburnin)
) # burn-in
jags_sample <- suppressWarnings(
jags.samples(
posterior,
c("theta",
"gamma",
"R0",
"beta",
"Y",
"lambdaY",
"lambdaR",
"k"),
n.iter = M,
thin = thn
)
)
if (dic) {
dic_val <- suppressWarnings(
dic.samples(posterior, n.iter = M, thin = thn)
)
} else {
dic_val <- NULL
}
if (save_files) {
png(paste0(file_add, casename, "theta_p.png"), width = 700, height = 900)
plot(as.mcmc.list(jags_sample$theta)[, (1:3) * (T_prime + 1)])
# posterior true probabilities
dev.off()
png(paste0(file_add, casename, "R0_p.png"), width = 700, height = 350)
R0_p <- unlist(as.mcmc.list(jags_sample$R0))
plot(as.mcmc.list(jags_sample$R0))
dev.off()
png(paste0(file_add, casename, "gamma_p.png"), width = 700, height = 350)
gamma_p <- unlist(as.mcmc.list(jags_sample$gamma))
plot(as.mcmc.list(jags_sample$gamma))
dev.off()
png(paste0(file_add, casename, "beta_p.png"), width = 700, height = 350)
beta_p <- unlist(as.mcmc.list(jags_sample$beta))
plot(as.mcmc.list(jags_sample$beta))
dev.off()
png(paste0(file_add, casename, "lambdaY_p.png"), width = 700, height = 350)
lambdaY_p <- unlist(as.mcmc.list(jags_sample$lambdaY))
plot(as.mcmc.list(jags_sample$lambdaY))
dev.off()
png(paste0(file_add, casename, "lambdaR_p.png"), width = 700, height = 350)
lambdaR_p <- unlist(as.mcmc.list(jags_sample$lambdaR))
plot(as.mcmc.list(jags_sample$lambdaR))
dev.off()
png(paste0(file_add, casename, "k_p.png"), width = 700, height = 350)
k_p <- unlist(as.mcmc.list(jags_sample$k))
plot(as.mcmc.list(jags_sample$k))
dev.off()
} else {
R0_p <- unlist(as.mcmc.list(jags_sample$R0))
gamma_p <- unlist(as.mcmc.list(jags_sample$gamma))
beta_p <- unlist(as.mcmc.list(jags_sample$beta))
lambdaY_p <- unlist(as.mcmc.list(jags_sample$lambdaY))
lambdaR_p <- unlist(as.mcmc.list(jags_sample$lambdaR))
k_p <- unlist(as.mcmc.list(jags_sample$k))
}
theta_p <- array(
Reduce(
rbind, as.mcmc.list(jags_sample$theta)
),
dim = c(len, T_prime + 1, 3)
)
theta_p_mean <- apply(theta_p[, T_prime + 1, ], 2, mean)
theta_p_ci <- as.vector(
apply(
theta_p[, T_prime + 1, ],
2,
quantile,
c(0.025, 0.5, 0.975)
)
)
R0_p_mean <- mean(R0_p)
R0_p_ci <- quantile(
R0_p,
c(0.025, 0.5, 0.975)
)
gamma_p_mean <- mean(gamma_p)
gamma_p_ci <- quantile(
gamma_p,
c(0.025, 0.5, 0.975)
)
beta_p_mean <- mean(beta_p)
beta_p_ci <- quantile(
beta_p,
c(0.025, 0.5, 0.975)
)
lambdaY_p_mean <- mean(lambdaY_p)
lambdaY_p_ci <- quantile(
lambdaY_p,
c(0.025, 0.5, 0.975)
)
lambdaR_p_mean <- mean(lambdaR_p)
lambdaR_p_ci <- quantile(
lambdaR_p,
c(0.025, 0.5, 0.975)
)
k_p_mean <- mean(k_p)
k_p_ci <- quantile(
k_p,
c(0.025, 0.5, 0.975)
)
#### Forecast ####
theta_pp <- array(0, dim = c(len, T_fin - T_prime, 3))
Y_pp <- matrix(NA, nrow = len, ncol = T_fin - T_prime)
R_pp <- matrix(NA, nrow = len, ncol = T_fin - T_prime)
for (l in 1:len) {
thetalt1 <- theta_p[l, T_prime + 1, 1]
thetalt2 <- theta_p[l, T_prime + 1, 2]
thetalt3 <- theta_p[l, T_prime + 1, 3]
betal <- c(beta_p)[l]
gammal <- c(gamma_p)[l]
kt <- c(k_p)[l]
lambdaYl <- c(lambdaY_p)[l]
lambdaRl <- c(lambdaR_p)[l]
if (betal < 0 | gammal < 0 | thetalt1 < 0 |
thetalt2 < 0 | thetalt3 < 0) next
for (t in 1:(T_fin - T_prime)) {
Km <- NULL
alpha_pp <- NULL
Km[1] <- -betal * pi[t + T_prime] * thetalt1 * thetalt2
Km[9] <- gammal * thetalt2
Km[5] <- -Km[1] - Km[9]
Km[2] <- -betal *
pi[t + T_prime] * (thetalt1 + 0.5 * Km[1]) * (thetalt2 + 0.5 * Km[5])
Km[10] <- gammal * (thetalt2 + 0.5 * Km[5])
Km[6] <- -Km[2] - Km[10]
Km[3] <- -betal *
pi[t + T_prime] * (thetalt1 + 0.5 * Km[2]) * (thetalt2 + 0.5 * Km[6])
Km[11] <- gammal * (thetalt2 + 0.5 * Km[6])
Km[7] <- -Km[3] - Km[11]
Km[4] <- - betal *
pi[t + T_prime] * (thetalt1 + Km[3]) * (thetalt2 + Km[7])
Km[12] <- gammal * (thetalt2 + Km[7])
Km[8] <- -Km[4] - Km[12]
alpha_pp[1] <- thetalt1 + (Km[1] + 2 * Km[2] + 2 * Km[3] + Km[4]) / 6
alpha_pp[2] <- thetalt2 + (Km[5] + 2 * Km[6] + 2 * Km[7] + Km[8]) / 6
alpha_pp[3] <- thetalt3 + (Km[9] + 2 * Km[10] + 2 * Km[11] + Km[12]) / 6
thetalt_tmp <- rdirichlet(1, kt * c(alpha_pp))
thetalt1 <- theta_pp[l, t, 1] <- thetalt_tmp[1]
thetalt2 <- theta_pp[l, t, 2] <- thetalt_tmp[2]
thetalt3 <- theta_pp[l, t, 3] <- thetalt_tmp[3]
Y_pp[l, t] <- rbeta(1, lambdaYl * thetalt2, lambdaYl * (1 - thetalt2))
R_pp[l, t] <- rbeta(1, lambdaRl * thetalt3, lambdaRl * (1 - thetalt3))
}
}
par(mfrow = c(1, 1))
col2 <- gg_color_hue(2)
Y_band <- data.frame(
t(apply(Y_pp,
2,
quantile,
probs = c(0.025, 0.975),
na.rm = TRUE
))
)
thetaI_band <- data.frame(
t(apply(theta_p[, -1, 2],
2,
quantile,
probs = c(0.025, 0.975),
na.rm = TRUE
))
)
Y_mean <- c(colMeans(Y_pp, na.rm = TRUE))
thetaI_mean <- c(
colMeans(theta_p[, -1, 2], na.rm = TRUE),
colMeans(theta_pp[, , 2], na.rm = TRUE)
)
thetaI_median <- c(
apply(theta_p[, -1, 2], 2, median, na.rm = TRUE),
apply(theta_pp[, , 2], 2, median, na.rm = TRUE)
)
colnames(Y_band) <- c("lower", "upper")
colnames(thetaI_band) <- c("lower", "upper")
data_pre <- data.frame(time = 1:T_prime, Y)
data_post <- data.frame(time = 1:T_prime, thetaI_band)
data_fore <- data.frame(time = (T_prime + 1):T_fin, Y_band, Y_mean)
data_comp <- data.frame(
time = 1:T_fin,
rbind(thetaI_band, Y_band),
phase = c(
rep("pre", nrow(thetaI_band)),
rep("post", nrow(Y_band))
),
mean = thetaI_mean,
median = thetaI_median
)
data_poly <- data.frame(
y = c(thetaI_band$upper,
rev(thetaI_band$lower),
Y_band$upper,
rev(Y_band$lower)),
x = c(1:T_prime,
T_prime:1,
(T_prime + 1):T_fin,
T_fin:(T_prime + 1)),
phase = c(rep("pre", T_prime * 2),
rep("post", (T_fin - T_prime) * 2)),
value = c(rep(col2[1], T_prime * 2),
rep(col2[2], (T_fin - T_prime) * 2))
)
## First-order derivative check
thetaS_mat <- cbind(theta_p[, -1, 1], theta_pp[,, 1])
thetaI_mat <- cbind(theta_p[, -1, 2], theta_pp[,, 2])
thetaR_mat <- cbind(theta_p[, -1, 3], theta_pp[,, 3])
dthetaI_mat_post <- (theta_pp[,, 1] * theta_pp[,, 2]) *
((c(beta_p)) %o% (pi[(T_prime + 1):T_fin])) -
theta_pp[,, 2] * replicate(T_fin - T_prime, c(gamma_p))
dthetaI_mat_pre <- t(
apply(
theta_p[,, 2],
1,
function(v) {
diff(smooth(v))
}
)
)
dthetaI_mat <- cbind(dthetaI_mat_pre, dthetaI_mat_post)
dthetaI <- colMeans(dthetaI_mat, na.rm = TRUE)
dthetaI_tp1 <- (1:T_fin)[which.max(dthetaI)] # first second order derivative=0
dthetaI_tp2 <- (dthetaI_tp1:T_fin)[which.min(dthetaI[dthetaI_tp1:T_fin] > 0)] # first order derivative=0
dthetaI_tp1_rd <- max(round(dthetaI_tp1), 1)
if (dthetaI_tp1_rd > T_prime) {
thetatI_tp1_vec <- thetaI_mat[, dthetaI_tp1_rd]
thetaI_tp1_mean <- mean(thetatI_tp1_vec, na.rm = TRUE)
thetaI_tp1_ci <- quantile(
thetatI_tp1_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
Y_tp1_vec <- Y_pp[, dthetaI_tp1_rd - T_prime]
Y_tp1_mean <- mean(Y_tp1_vec, na.rm = TRUE)
Y_tp1_ci <- quantile(
Y_tp1_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
thetatR_tp1_vec <- thetaR_mat[, dthetaI_tp1_rd]
thetaR_tp1_mean <- mean(thetatR_tp1_vec, na.rm = TRUE)
thetaR_tp1_ci <- quantile(
thetatR_tp1_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
R_tp1_vec <- R_pp[, dthetaI_tp1_rd - T_prime]
R_tp1_mean <- mean(R_tp1_vec, na.rm = TRUE)
R_tp1_ci <- quantile(
R_tp1_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
} else {
thetatI_tp1_vec <- thetaI_mat[, dthetaI_tp1_rd]
thetaI_tp1_mean <- mean(thetatI_tp1_vec, na.rm = TRUE)
thetaI_tp1_ci <- quantile(
thetatI_tp1_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
Y_tp1_vec <- NA
Y_tp1_mean <- mean(Y_tp1_vec, na.rm = TRUE)
Y_tp1_ci <- quantile(
Y_tp1_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
thetatR_tp1_vec <- thetaR_mat[, dthetaI_tp1_rd]
thetaR_tp1_mean <- mean(thetatR_tp1_vec, na.rm = TRUE)
thetaR_tp1_ci <- quantile(
thetatR_tp1_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
R_tp1_vec <- R_pp[, dthetaI_tp1_rd - T_prime]
R_tp1_mean <- mean(R_tp1_vec, na.rm = TRUE)
R_tp1_ci <- quantile(
R_tp1_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
}
dthetaI_tp2_rd <- max(round(dthetaI_tp2), 1)
if (dthetaI_tp1_rd == dthetaI_tp2_rd) {
thetatI_tp2_vec <- NA
thetaI_tp2_mean <- mean(thetatI_tp2_vec, na.rm = TRUE)
thetaI_tp2_ci <- quantile(
thetatI_tp2_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
Y_tp2_vec <- NA
Y_tp2_mean <- mean(Y_tp2_vec, na.rm = TRUE)
Y_tp2_ci <- quantile(
Y_tp2_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
thetatR_tp2_vec <- NA
thetaR_tp2_mean <- mean(thetatR_tp2_vec, na.rm = TRUE)
thetaR_tp2_ci <- quantile(
thetatR_tp2_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
R_tp2_vec <- NA
R_tp2_mean <- mean(R_tp2_vec, na.rm = TRUE)
R_tp2_ci <- quantile(
R_tp2_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
} else if (dthetaI_tp2_rd > T_prime) {
thetatI_tp2_vec <- thetaI_mat[, dthetaI_tp2_rd]
thetaI_tp2_mean <- mean(thetatI_tp2_vec, na.rm = TRUE)
thetaI_tp2_ci <- quantile(
thetatI_tp2_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
Y_tp2_vec <- Y_pp[, dthetaI_tp2_rd - T_prime]
Y_tp2_mean <- mean(Y_tp2_vec, na.rm = TRUE)
Y_tp2_ci <- quantile(
Y_tp2_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
thetatR_tp2_vec <- thetaR_mat[, dthetaI_tp2_rd]
thetaR_tp2_mean <- mean(thetatR_tp2_vec, na.rm = TRUE)
thetaR_tp2_ci <- quantile(
thetatR_tp2_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
R_tp2_vec <- R_pp[, dthetaI_tp2_rd - T_prime]
R_tp2_mean <- mean(R_tp2_vec, na.rm = TRUE)
R_tp2_ci <- quantile(
R_tp2_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
} else {
thetatI_tp2_vec <- thetaI_mat[, dthetaI_tp2_rd]
thetaI_tp2_mean <- mean(thetatI_tp2_vec, na.rm = TRUE)
thetaI_tp2_ci <- quantile(
thetatI_tp2_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
Y_tp2_vec <- NA
Y_tp2_mean <- mean(Y_tp2_vec, na.rm = TRUE)
Y_tp2_ci <- quantile(
Y_tp2_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
thetatR_tp2_vec <- thetaR_mat[, dthetaI_tp2_rd]
thetaR_tp2_mean <- mean(thetatR_tp2_vec, na.rm = TRUE)
thetaR_tp2_ci <- quantile(
thetatR_tp2_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
R_tp2_vec <- NA
R_tp2_mean <- mean(R_tp2_vec, na.rm = TRUE)
R_tp2_ci <- quantile(
R_tp2_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
}
thetaR_max_vec <- thetaR_mat[, T_fin]
thetaR_max_mean <- mean(thetaR_max_vec)
thetaR_max_ci <- quantile(
thetaR_max_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
cumInf_vec <- thetaR_mat[, T_fin] + thetaI_mat[, T_fin]
cumInf_mean <- mean(cumInf_vec)
cumInf_ci <- quantile(
cumInf_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
dthetaI_tp1_date <- chron_ls[dthetaI_tp1]
dthetaI_tp2_date <- chron_ls[dthetaI_tp2]
incidence_vec <- 1 - thetaS_mat[, T_fin]
incidence_mean <- mean(incidence_vec, na.rm = TRUE)
incidence_ci <- quantile(
incidence_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
first_tp_vec <- (1:T_fin)[apply(dthetaI_mat, 1, which.max)] # first second order derivative=0
second_tp_vec <- sapply(
1:len,
function(l) {
(first_tp_vec[l]:T_fin)[
which.min(dthetaI_mat[l, first_tp_vec[l]:T_fin] > 0)
]
}
)
end_p_vec <- sapply(
1:len,
function(l) {
if (any(thetaI_mat[l, first_tp_vec[l]:T_fin] <= eps, na.rm = TRUE)) {
(first_tp_vec[l]:T_fin)[
which.max(thetaI_mat[l, first_tp_vec[l]:T_fin] <= eps)
]
} else {
T_fin
}
}
)
# first order derivative=0
first_tp_mean <- mean(first_tp_vec, na.rm = TRUE)
second_tp_mean <- mean(second_tp_vec, na.rm = TRUE)
end_p_mean <- mean(end_p_vec, na.rm = TRUE)
first_tp_ci <- quantile(
first_tp_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
second_tp_ci <- quantile(
second_tp_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
end_p_ci <- quantile(
end_p_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
first_tp_date_mean <- chron_ls[first_tp_mean]
second_tp_date_mean <- chron_ls[second_tp_mean]
end_p_date_mean <- chron_ls[end_p_mean]
first_tp_date_ci <- chron_ls[first_tp_ci]
second_tp_date_ci <- chron_ls[second_tp_ci]
end_p_date_ci <- chron_ls[end_p_ci]
names(first_tp_date_ci) <- c("2.5%", "50%", "97.5%")
names(second_tp_date_ci) <- c("2.5%", "50%", "97.5%")
names(end_p_date_ci) <- c("2.5%", "50%", "97.5%")
if (save_files) {
png(paste0(file_add, casename, "deriv.png"), width = 700, height = 350)
plot(
y = dthetaI,
x = chron_ls,
type = "l",
ylab = "1st order derivative",
main = "Infection Prevalence"
)
abline(h = 0, col = 2)
abline(v = chron_ls[T_prime], col = "blue")
if (exponential == FALSE & !is.null(change_time) & !is.null(pi0)) {
abline(
v = change_time_chorn[-c(1, length(change_time_chorn))],
col = "gray"
)
}
legend(
"topright",
legend = c("Last observation", "change point"),
col = c("blue", "gray"),
lty = 1, title = "",
bty = "n"
)
dev.off()
}
## Prepare the Spaghetti plot
sample_dthetaI_mat <- cbind(
dthetaI_mat[sample.int(len, 20, replace = FALSE), ]
)
colnames(sample_dthetaI_mat) <- c(as.character(chron_ls)[-1])
sample_dthetaI_mat_long <- reshape2::melt(sample_dthetaI_mat)
colnames(sample_dthetaI_mat_long) <- c("id", "date", "dthetaI")
sample_dthetaI_mat_long$date <- (chron(as.character(
sample_dthetaI_mat_long$date)))
dthetaI_mean_data <- data.frame(dthetaI, date = chron_ls[-1])
spaghetti_ht <- mean(range(sample_dthetaI_mat)) / 2
spaghetti_plot <- ggplot() +
theme_bw() +
theme(
plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
plot.subtitle = element_text(hjust = 0.5, size = 12),
axis.text.x = element_text(angle = 45, hjust = 1),
axis.text = element_text(size = 15),
axis.title = element_text(size = 14)
) +
geom_line(
data = sample_dthetaI_mat_long,
aes(x = date, y = dthetaI, group = id, color = id)
) +
scale_color_gradientn(colours = rainbow(5, alpha = 0.5)) +
labs(
title = "spaghetti plot of infection prevalence",
x = "time", y = "1st order derivative"
) +
geom_line(
data = dthetaI_mean_data,
aes(x = date, y = dthetaI), color = 1
) +
scale_x_continuous(
labels = as.character(chron_ls)[seq(1, T_fin, 30)],
breaks = as.numeric(chron_ls[-1][seq(1, T_fin, 30)])
) +
annotate(
geom = "text",
label = as.character(chron(chron_ls[T_prime]), format = "mon day"),
x = as.numeric(chron_ls[T_prime]) + 12,
y = spaghetti_ht,
color = "blue"
) +
annotate(
geom = "text",
label = as.character(chron(dthetaI_tp1_date, format = "mon day")),
x = as.numeric(dthetaI_tp1_date) + 12,
y = spaghetti_ht * 1.25,
color = "darkgreen"
) +
geom_vline(
xintercept = as.numeric(chron_ls[T_prime]),
color = "blue",
show.legend = TRUE
) +
geom_vline(
xintercept = as.numeric(dthetaI_tp1_date),
color = "darkgreen",
show.legend = TRUE
)
# if(dthetaI_tp1>T_prime)
spaghetti_plot <- spaghetti_plot +
geom_rect(
data = data.frame(
xmin = as.numeric(first_tp_date_ci[1]),
xmax = as.numeric(first_tp_date_ci[3]),
ymin = -Inf,
ymax = Inf, ci = "first tp"
),
aes(
xmin = xmin,
xmax = xmax,
ymin = ymin,
ymax = ymax
),
fill = "darkgreen",
alpha = 0.15
)
# if(dthetaI_tp2>T_prime)
spaghetti_plot <- spaghetti_plot +
geom_rect(
data = data.frame(
xmin = as.numeric(second_tp_date_ci[1]),
xmax = as.numeric(second_tp_date_ci[3], ci = "second tp"),
ymin = -Inf, ymax = Inf
),
aes(
xmin = xmin,
xmax = xmax,
ymin = ymin,
ymax = ymax
),
fill = "purple",
alpha = 0.15
)
if (dthetaI_tp2_date > dthetaI_tp1_date) {
spaghetti_plot <- spaghetti_plot +
geom_vline(
xintercept = as.numeric(dthetaI_tp2_date),
color = "purple",
show.legend = TRUE
) +
annotate(
geom = "text",
label = as.character(chron(dthetaI_tp2_date, format = "mon day")),
x = as.numeric(dthetaI_tp2_date) + 12,
y = spaghetti_ht * 1.5, color = "purple"
)
}
if (save_files) {
ggsave(paste0(
file_add, casename, "_spaghetti.png"
),
width = 12,
height = 10
)
}
###################
y_text_ht <- max(rbind(thetaI_band, Y_band),
na.rm = TRUE) / 2
plot1 <- ggplot(
data = data_poly,
aes(x = x, y = y)
) +
geom_polygon(
alpha = 0.5,
aes(fill = value, group = phase)
) +
labs(
title = substitute(
paste(
casename, ": infection forecast with prior ",
beta[0], "=", v1, ",",
gamma[0], "=", v2, " and ",
R[0], "=", v3
),
list(
casename = casename,
v1 = format(beta0, digits = 3),
v2 = format(gamma0, digits = 3),
v3 = format(R0, digits = 3)
)
),
subtitle = substitute(
paste(
"Posterior ",
beta[p], "=", v1, ",",
gamma[p], "=", v2, " and ",
R[0], "=", v3
),
list(
v1 = format(beta_p_mean, digits = 3),
v2 = format(gamma_p_mean, digits = 3),
v3 = format(R0_p_mean, digits = 3)
)
),
x = "time", y = "P(Infected)"
) +
geom_line(data = data_comp, aes(x = time, y = median), color = "red") +
geom_vline(xintercept = T_prime, color = "blue", show.legend = TRUE) +
geom_vline(
xintercept = dthetaI_tp1,
color = "darkgreen",
show.legend = TRUE
) +
geom_line(data = data_comp, aes(x = time, y = mean), color = "darkgray") +
geom_point(data = data_pre, aes(x = time, y = Y)) +
theme_bw() +
theme(
plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
plot.subtitle = element_text(hjust = 0.5, size = 12),
axis.text.x = element_text(angle = 45, hjust = 1),
axis.text = element_text(size = 15),
axis.title = element_text(size = 14),
legend.title = element_text(size = 15),
legend.text = element_text(size = 12)
) +
scale_x_continuous(
labels = as.character(chron_ls)[seq(1, T_fin, 30)],
breaks = seq(1, T_fin, 30)
) +
scale_fill_discrete(
name = "Posterior",
labels = c(
expression(paste(y[t[0] + 1:T]^I, " | ", y[1:t[0]]^I, ", ", y[1:t[0]]^R)),
expression(paste(theta[1:t[0]]^I, " | ", y[1:t[0]]^I, ", ", y[1:t[0]]^R))
)
) +
annotate(
geom = "text",
label = as.character(
chron(chron_ls[T_prime]),
format = "mon day"
),
x = T_prime + 12,
y = y_text_ht,
color = "blue"
) +
annotate(
geom = "text",
label = as.character(chron(
dthetaI_tp1_date,
format = "mon day"
)),
x = dthetaI_tp1 + 12,
y = y_text_ht * 1.25,
color = "darkgreen"
)
if (dthetaI_tp2 > dthetaI_tp1) {
plot1 <- plot1 + geom_vline(
xintercept = dthetaI_tp2,
color = "purple",
show.legend = TRUE
) +
annotate(
geom = "text",
label = as.character(chron(dthetaI_tp2_date, format = "mon day")),
x = dthetaI_tp2 + 12, y = y_text_ht * 1.5, color = "purple"
)
}
if (save_files) ggsave(paste0(file_add, casename, "_forecast.png"), width = 12, height = 10)
### Removed
R_band <- data.frame(t(apply(R_pp, 2, quantile, probs = c(0.025, 0.975), na.rm = TRUE)))
thetaR_band <- data.frame(t(apply(theta_p[, -1, 3], 2, quantile, probs = c(0.025, 0.975), na.rm = TRUE)))
R_mean <- c(colMeans(R_pp, na.rm = TRUE))
thetaR_mean <- c(colMeans(theta_p[, -1, 3], na.rm = TRUE), colMeans(theta_pp[, , 3], na.rm = TRUE))
thetaR_med <- c(apply(theta_p[, -1, 3], 2, median, na.rm = TRUE), apply(theta_pp[, , 3], 2, median, na.rm = TRUE))
colnames(R_band) <- c("lower", "upper")
colnames(thetaR_band) <- c("lower", "upper")
data_pre_R <- data.frame(time = 1:T_prime, R) # previous data
data_post_R <- data.frame(time = 1:T_prime, thetaR_band) # posterior of theta^R
data_fore_R <- data.frame(time = (T_prime + 1):T_fin, R_band, R_mean) # The forecast of R after T_prime
data_comp_R <- data.frame(time = 1:T_fin, rbind(thetaR_band, R_band), phase = c(rep("pre", nrow(thetaR_band)), rep("post", nrow(R_band))), mean = thetaR_mean, median = thetaR_med, dead = thetaR_mean * death_in_R, dead_med = thetaR_med * death_in_R) # the filled area--polygon
data_poly_R <- data.frame(y = c(thetaR_band$upper, rev(thetaR_band$lower), R_band$upper, rev(R_band$lower)), x = c(1:T_prime, T_prime:1, (T_prime + 1):T_fin, T_fin:(T_prime + 1)), phase = c(rep("pre", T_prime * 2), rep("post", (T_fin - T_prime) * 2)), value = c(rep(col2[1], T_prime * 2), rep(col2[2], (T_fin - T_prime) * 2)))
r_text_ht <- max(rbind(thetaR_band, R_band), na.rm = TRUE) / 2
plot2 <- ggplot(data = data_poly_R, aes(x = x, y = y)) +
geom_polygon(alpha = 0.5, aes(fill = value, group = phase)) +
labs(
title = substitute(
paste(
casename,
": removed forecast with prior ", beta[0], "=", v1, ",",
gamma[0], "=", v2, " and ", R[0], "=", v3
),
list(
casename = casename, v1 = format(beta0, digits = 3),
v2 = format(gamma0, digits = 3), v3 = format(R0, digits = 3)
)
),
subtitle = substitute(
paste(
"posterior: ", beta[p], "=", v1, ",",
gamma[p], "=", v2, " and ", R[0], "=", v3
),
list(
v1 = format(beta_p_mean, digits = 3), v2 = format(gamma_p_mean, digits = 3),
v3 = format(R0_p_mean, digits = 3)
)
), x = "time", y = "P(Removed)"
) +
geom_line(data = data_comp_R, aes(x = time, y = median), color = "red", linetype = 1) +
geom_vline(xintercept = T_prime, color = "blue") +
geom_vline(xintercept = dthetaI_tp1, color = "darkgreen") +
geom_line(data = data_comp_R, aes(x = time, y = mean), color = "darkgray") +
geom_point(data = data_pre_R, aes(x = time, y = R)) +
theme_bw() +
theme(
plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
plot.subtitle = element_text(hjust = 0.5, size = 12),
axis.text.x = element_text(angle = 45, hjust = 1),
axis.text = element_text(size = 15),
axis.title = element_text(size = 14),
legend.title = element_text(size = 15),
legend.text = element_text(size = 12)
) +
scale_x_continuous(labels = as.character(chron_ls)[seq(1, T_fin, 30)], breaks = seq(1, T_fin, 30)) +
scale_fill_discrete(
name = "Posterior",
labels = c(
expression(paste(y[t[0] + 1:T]^R, " | ", y[1:t[0]]^I, ", ", y[1:t[0]]^R)),
expression(paste(theta[1:t[0]]^R, " | ", y[1:t[0]]^I, ", ", y[1:t[0]]^R))
)
) +
annotate(geom = "text", label = as.character(chron(chron_ls[T_prime]), format = "mon day"), x = T_prime + 12, y = r_text_ht, color = "blue") +
annotate(geom = "text", label = as.character(chron(dthetaI_tp1_date, format = "mon day")), x = dthetaI_tp1 + 12, y = r_text_ht * 1.25, color = "darkgreen")
if (dthetaI_tp2 > dthetaI_tp1) {
plot2 <- plot2 + geom_vline(xintercept = dthetaI_tp2, color = "purple", show.legend = TRUE) + annotate(geom = "text", label = as.character(chron(dthetaI_tp2_date, format = "mon day")), x = dthetaI_tp2 + 12, y = r_text_ht * 1.5, color = "purple")
}
if (add_death) plot2 <- plot2 + geom_line(data = data_comp_R, aes(x = time, y = dead), color = "black", linetype = 1) + geom_line(data = data_comp_R, aes(x = time, y = dead_med), color = "black", linetype = 2)
if (save_files) ggsave(paste0(file_add, casename, "_forecast2.png"), width = 12, height = 10)
out_table1 <- data.frame(matrix(c(theta_p_mean, theta_p_ci, R0_p_mean, R0_p_ci, gamma_p_mean, gamma_p_ci, beta_p_mean, beta_p_ci, incidence_mean = incidence_mean, incidence_ci = incidence_ci, thetaI_tp1_mean = thetaI_tp1_mean, thetaI_tp1_ci = thetaI_tp1_ci, thetaR_tp1_mean = thetaR_tp1_mean, thetaR_tp1_ci = thetaR_tp1_ci, Y_tp1_mean = Y_tp1_mean, Y_tp1_ci = Y_tp1_ci, R_tp1_mean = R_tp1_mean, R_tp1_ci = R_tp1_ci, thetaI_tp2_mean = thetaI_tp2_mean, thetaI_tp2_ci = thetaI_tp2_ci, thetaR_tp2_mean = thetaR_tp2_mean, thetaR_tp2_ci = thetaR_tp2_ci, Y_tp2_mean = Y_tp2_mean, Y_tp2_ci = Y_tp2_ci, R_tp2_mean = R_tp2_mean, R_tp2_ci = R_tp2_ci, thetaR_max_mean, thetaR_max_ci, cumInf_mean = cumInf_mean, cumInf_ci = cumInf_ci), nrow = 1))
out_table2 <- data.frame(matrix(c(dthetaI_tp1_date = as.character(dthetaI_tp1_date), first_tp_mean = as.character(first_tp_date_mean), first_tp_ci = as.character(first_tp_date_ci), dthetaI_tp2_date = as.character(dthetaI_tp2_date), second_tp_mean = as.character(second_tp_date_mean), second_tp_ci = as.character(second_tp_date_ci), end_p_date_mean = as.character(end_p_date_mean), end_p_date_ci = as.character(end_p_date_ci), begin_str = begin_str), nrow = 1))
out_table <- cbind(out_table1, out_table2)
colnames(out_table) <- c(
"thetaS_last_obs_p_mean",
"thetaI_last_obs_p_mean",
"thetaR_last_obs_p_mean",
"thetaS_last_obs_p_ci_low",
"thetaS_last_obs_p_ci_med",
"thetaS_last_obs_p_ci_up",
"thetaI_last_obs_p_ci_low",
"thetaI_last_obs_p_ci_med",
"thetaI_last_obs_p_ci_up",
"thetaR_last_obs_p_ci_low",
"thetaR_last_obs_p_ci_med",
"thetaR_last_obs_p_ci_up",
"R0_p_mean",
"R0_p_ci_low",
"R0_p_ci_med",
"R0_p_ci_up",
"gamma_p_mean",
"gamma_p_ci_low",
"gamma_p_ci_med",
"gamma_p_ci_up",
"beta_p_mean",
"beta_p_ci_low",
"beta_p_ci_med",
"beta_p_ci_up",
"incidence_mean",
"incidence_ci_low",
"incidence_ci_median",
"incidence_ci_up",
"thetaI_tp1_mean",
"thetaI_tp1_ci_low",
"thetaI_tp1_ci_med",
"thetaI_tp1_ci_up",
"thetaR_tp1_mean",
"thetaR_tp1_ci_low",
"thetaR_tp1_ci_med",
"thetaR_tp1_ci_up",
"Y_tp1_mean",
"Y_tp1_ci_low",
"Y_tp1_ci_med",
"Y_tp1_ci_up",
"R_tp1_mean",
"R_tp1_ci_low",
"R_tp1_ci_med",
"R_tp1_ci_up",
"thetaI_tp2_mean",
"thetaI_tp2_ci_low",
"thetaI_tp2_ci_med",
"thetaI_tp2_ci_up",
"thetaR_tp2_mean",
"thetaR_tp2_ci_low",
"thetaR_tp2_ci_med",
"thetaR_tp2_ci_up",
"Y_tp2_mean",
"Y_tp2_ci_low",
"Y_tp2_ci_med",
"Y_tp2_ci_up",
"R_tp2_mean",
"R_tp2_ci_low",
"R_tp2_ci_med",
"R_tp2_ci_up",
"thetaR_max_mean",
"thetaR_max_ci_low",
"thetaR_max_ci_med",
"thetaR_max_ci_up",
"cumInf_mean",
"cumInf_ci_low",
"cumInf_ci_med",
"cumInf_ci_up",
"dthetaI_tp1_date",
"first_tp_mean",
"first_tp_ci_low",
"first_tp_ci_med",
"first_tp_ci_up",
"dthetaI_tp2_date",
"second_tp_mean",
"second_tp_ci_low",
"second_tp_ci_med",
"second_tp_ci_up",
"end_p_mean",
"end_p_ci_low",
"end_p_ci_med",
"end_p_ci_up",
"begin_str"
)
if (save_files) write.csv(out_table, file = paste0(file_add, casename, "_summary.csv"))
if (save_mcmc) {
save(
theta_pp,
Y_pp,
R_pp,
file = paste0(file_add, casename, "_forecast_MCMC.RData")
)
save(
jags_sample,
file = paste0(file_add, casename, "_rjags_MCMC.RData")
)
}
if (save_plot_data) {
other_plot <- list(T_prime = T_prime, T_fin = T_fin, chron_ls = chron_ls, dthetaI_tp1 = dthetaI_tp1, dthetaI_tp2 = dthetaI_tp2, dthetaI_tp1_date = dthetaI_tp1_date, dthetaI_tp2_date = dthetaI_tp2_date, beta_p_mean = beta_p_mean, gamma_p_mean = gamma_p_mean, R0_p_mean = R0_p_mean)
spaghetti_plot_ls <- list(spaghetti_ht = spaghetti_ht, dthetaI_mean_data = dthetaI_mean_data, sample_dthetaI_mat_long = sample_dthetaI_mat_long, first_tp_date_ci = first_tp_date_ci, second_tp_date_ci = second_tp_date_ci)
infection_plot_ls <- list(y_text_ht = y_text_ht, data_poly = data_poly, data_comp = data_comp, data_pre = data_pre)
removed_plot_ls <- list(r_text_ht = r_text_ht, data_poly_R = data_poly_R, data_comp_R = data_comp_R, data_pre_R = data_pre_R)
plot_data_ls <- list(casename = casename, other_plot = other_plot, spaghetti_plot_ls = spaghetti_plot_ls, infection_plot_ls = infection_plot_ls, removed_plot_ls = removed_plot_ls)
save(plot_data_ls, file = paste0(file_add, casename, "_plot_data.RData"))
}
# Gelman-Rubin convergence diagnostics
gelman_out_name <- c("R0", "gamma", "beta", "theta", "lambdaY", "lambdaR", "k")
gelman_diag_list <- lapply(gelman_out_name, function(name) {
if (name == "theta") {
tryCatch(
gelman.diag(as.mcmc.list(jags_sample[[name]])[, (1:3) * (T_prime + 1)],
autoburnin = FALSE
),
error = function(e) e
)
} else {
tryCatch(
gelman.diag(as.mcmc.list(jags_sample[[name]]),
autoburnin = FALSE
),
error = function(e) e
)
}
})
if (save_files) {
sink(paste0(file_add, casename, "_Gelman_diag.txt"))
print(gelman_diag_list)
sink()
}
res <- list(
casename = casename,
incidence_mean = incidence_mean,
incidence_ci = incidence_ci,
out_table = out_table,
plot_infection = plot1,
plot_removed = plot2,
spaghetti_plot = spaghetti_plot,
first_tp_mean = as.character(first_tp_date_mean),
first_tp_ci = as.character(first_tp_date_ci),
second_tp_mean = as.character(second_tp_date_mean),
second_tp_ci = as.character(second_tp_date_ci),
dic_val = dic_val,
gelman_diag_list = gelman_diag_list
)
res
}
if (FALSE) {
NI_complete <- c(
41, 41, 41, 45, 62, 131, 200, 270, 375, 444, 549, 729,
1052, 1423, 2714, 3554, 4903, 5806, 7153, 9074, 11177,
13522, 16678, 19665, 22112, 24953, 27100, 29631, 31728, 33366
)
RI_complete <- c(
1, 1, 7, 10, 14, 20, 25, 31, 34, 45, 55, 71, 94, 121, 152, 213,
252, 345, 417, 561, 650, 811, 1017, 1261, 1485, 1917, 2260,
2725, 3284, 3754
)
N <- 58.5e6
R <- RI_complete / N
Y <- NI_complete / N - R # Jan13->Feb 11
### Step function of pi(t)
change_time <- c("01/23/2020", "02/04/2020", "02/08/2020")
pi0 <- c(1.0, 0.9, 0.5, 0.1)
res.step <- tvt.eSIR(Y, R,
begin_str = "01/13/2020", death_in_R = 0.4, T_fin = 200,
pi0 = pi0, change_time = change_time, dic = FALSE, casename = "Hubei_step", save_files = FALSE, save_mcmc = FALSE,
M = 5e2, nburnin = 2e2
)
res.step$plot_infection
res.step$plot_removed
# res.step$dic_val
res.step$gelman_diag_list
### continuous exponential function of pi(t)
res.exp <- tvt.eSIR(Y, R,
begin_str = "01/13/2020", death_in_R = 0.4, T_fin = 200, exponential = TRUE, dic = FALSE, lambda0 = 0.05,
casename = "Hubei_exp", save_files = TRUE, save_mcmc = FALSE,
M = 5e2, nburnin = 2e2
)
res.exp$plot_infection
# res.exp$plot_removed
### without pi(t), the standard state-space SIR model without intervention
res.nopi <- tvt.eSIR(Y, R,
begin_str = "01/13/2020", death_in_R = 0.4, T_fin = 200,
casename = "Hubei_nopi", save_files = FALSE,
M = 5e2, nburnin = 2e2
)
res.nopi$plot_infection
# res.nopi$plot_removed
}
Try the eSIR package in your browser
Any scripts or data that you put into this service are public.
eSIR documentation built on Dec. 11, 2021, 5:07 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions tvt.eSIR
Documented in tvt.eSIR
# This is the source code for R package eSIR: extended-SIR
# Built on Feb 13, 2020, and last edited on Oct 27, 2021
# Correspondence : Peter X.K. Song, Ph.D. (pxsong@umich.edu)
# Creator: Lili Wang, Ph.D. (lilywang@umich.edu)
# Model 1: An extended SIR model with a time-varying transmission rate modifier
# pi(t)
#' Fit extended state-space SIR model with time-varying transmission rates
#'
#' Fit extended state-space SIR model with pre-specified changes in the transmission rate, either stepwise or continuous, accommodating time-varying quarantine protocols.
#'
#' We fit a state-space model with extended SIR, in which a time-varying transmission rate modifier \eqn{\pi(t)} (between 0 and 1) is introduced to model. This allows us to accommodate quarantine protocol changes and ignored resources of hospitalization. The form of reducing rate may be a step-function with jumps at times of big policy changes or a smooth exponential survival function \eqn{\exp(-\lambda_0t)}. The parameters of the function and change points, if any, should be predefined.
#'
#'
#' @param Y the time series of daily observed infected compartment proportions.
#' @param R the time series of daily observed removed compartment proportions, including death and recovered.
#' @param pi0 the time-dependent transmission rate modifier \eqn{\pi(t)} between 0 and 1.
#' @param change_time the change points over time for step function pi, defalt value is \code{NULL}.
#' @param exponential logical, whether \eqn{\pi(t)} is exponential \eqn{\exp(-\lambda_0t)} or not; the default is \code{FALSE}.
#' @param lambda0 the rate of decline in the exponential survival function in \eqn{\exp(-\lambda_0t)}.
#' @param begin_str the character of starting time, the default is "01/13/2020".
#' @param T_fin the end of follow-up time after the beginning date \code{begin_str}, the default is 200.
#' @param nchain the number of MCMC chains generated by \code{\link[rjags]{rjags}}, the default is 4.
#' @param nadapt the iteration number of adaptation in the MCMC. We recommend using at least the default value 1e4 to obtained fully adapted chains.
#' @param M the number of draws in each chain, with no thinning. The default is M=5e2 but suggest using 5e5.
#' @param thn the thinning interval between mixing. The total number of draws thus would become \code{round(M/thn)*nchain}. The default is 10.
#' @param nburnin the burn-in period. The default is 2e2 but suggest 2e5.
#' @param dic logical, whether compute the DIC (deviance information criterion) for model selection.
#' @param death_in_R the numeric value of average of cumulative deaths in the removed compartments. The default is 0.4 within Hubei and 0.02 outside Hubei.
#' @param casename the string of the job's name. The default is "tvt.eSIR".
#' @param beta0 the hyperparameter of average transmission rate, the default is the one estimated from the SARS first-month outbreak (0.2586).
#' @param gamma0 the hyperparameter of average removed rate, the default is the one estimated from the SARS first-month outbreak (0.0821).
#' @param R0 the hyperparameter of the mean reproduction number R0. The default is thus the ratio of \code{beta0/gamma0}, which can be specified directly.
#' @param gamma0_sd the standard deviation for the prior distrbution of the removed rate \eqn{\gamma}, the default is 0.1.
#' @param R0_sd the standard deviation for the prior disbution of R0, the default is 1.
#' @param file_add the string to denote the location of saving output files and tables.
#' @param save_files logical, whether save (\code{TRUE}) results or not (\code{FALSE}). This enables to save summary tables, trace plots, and plots of the posterior means of the first-order derivatives of the infection prevalence process \eqn{\theta_t^I}.
#' @param save_mcmc logical, whether save (\code{TRUE}) all the MCMC outputs or not (\code{FALSE}).The output file will be an \code{.RData} file named by the \eqn{casename}. We include arrays of prevalence values of the three compartments with their matrices of posterior draws up to the last date of the collected data as \code{theta_p[,,1]} and afterwards as \code{theta_pp[,,1]} for \eqn{\theta_t^S}, \code{theta_p[,,2]} and \code{theta_pp[,,2]} for \eqn{\theta_t^I}, and \code{theta_p[,,3]} and \code{theta_pp[,,3]} for \eqn{\theta_t^R}. Moreover, the input and predicted proportions \code{Y}, \code{Y_pp}, \code{R} and \code{R_pp} can also be retrieved. The prevalence and predicted proportion matrices have rows for MCMC replicates, and columns for days. The MCMC posterior draws of other parameters including \code{beta_p}, \code{gamma_p}, \code{R0_p}, and variance controllers \code{k_p}, \code{lambdaY_p}, \code{lambdaR_p} are also available.
#' @param save_plot_data logical, whether save the plotting data or not.
#' @param save_files logical, whether to save plots to file.
#' @param add_death logical, whether add the approximate death curve to the plot, default is false.
#' @param eps a non-zero controller so that all the input \code{Y} and \code{R} values would be bounded above 0 (at least \code{eps}). Its default value is 1e-10
#'
#' @return
#' \item{casename}{the predefined \code{casename}.}
#' \item{incidence_mean}{mean cumulative incidence, the mean prevalence of cumulative confirmed cases at the end of the study.}
#' \item{incidence_ci}{2.5\%, 50\%, and 97.5\% quantiles of the incidences.}
#' \item{out_table}{summary tables including the posterior mean of the prevalance processes of the 3 states compartments (\eqn{\theta_t^S,\theta_t^I,\theta_t^R}) at last date of data collected ((\eqn{t^\prime}) decided by the lengths of your input data \code{Y} and \code{R}), and their respective credible inctervals (ci); the respective means and ci's of the reporduction number (R0), removed rate (\eqn{\gamma}), transmission rate (\eqn{\beta}).}
#' \item{plot_infection}{plot of summarizing and forecasting for the infection compartment, in which the vertial blue line denotes the last date of data collected (\eqn{t^\prime}), the vertial darkgray line denotes the deacceleration point (first turning point) that the posterior mean first-derivative of infection prevalence \eqn{\dot{\theta}_t^I} achieves the maximum, the vertical purple line denotes the second turning point that the posterior mean first-derivative infection proportion \eqn{\dot{\theta}_t^I} equals zero, the darkgray line denotes the posterior mean of the infection prevalence \eqn{\theta_t^I} and the red line denotes its posterior median. }
#' \item{plot_removed}{plot of summarizing and forecasting for the removed compartment with lines similar to those in the \code{plot_infection}. The vertical lines are identical, but the horizontal mean and median correspond to the posterior mean and median of the removed process \eqn{\theta_t^R}. An additional line indicates the estimated death prevalence from the input \code{death_in_R}.}
#' \item{spaghetti_plot}{20 randomly selected MCMC draws of the first-order derivative of the posterior prevalence of infection, namely \eqn{\dot{\theta}_t^I}. The black curve is the posterior mean of the derivative, and the vertical lines mark times of turning points corresponding respectively to those shown in \code{plot_infection} and \code{plot_removed}. Moreover, the 95\% credible intervals of these turning points are also highlighted by semi-transparent rectangles. }
#' \item{first_tp_mean}{the date t at which \eqn{\ddot{\theta}_t^I=0}, calculated as the average of the time points with maximum posterior first-order derivatives \eqn{\dot{\theta}_t^I}; this value may be slightly different from the one labeled by the "darkgreen" lines in the two plots \code{plot_infection} and \code{plot_removed}, which indicate the stationary point such that the first-order derivative of the averaged posterior of \eqn{\theta_t^I} reaches its maximum.}
#'
#' \item{first_tp_ci}{fwith \code{first_tp_mean}, it reports the corresponding credible interval and median.}
#' \item{second_tp_mean}{the date t at which \eqn{\theta_t^I=0}, calculated as the average of the stationary points of all of posterior first-order derivatives \eqn{\dot{\theta}_t^I}; this value may be slightly different from the one labeled by the "pruple" lines in the plots of \code{plot_infection} and \code{plot_removed}. The latter indicate stationary t at which the first-order derivative of the averaged posterior of \eqn{\theta_t^I} equals zero.}
#' \item{second_tp_ci}{with \code{second_tp_mean}, it reports the corresponding credible interval and median.}
#' \item{dic_val}{the output of \code{dic.samples()} in \code{\link[rjags]{dic.samples}}, computing deviance information criterion for model comparison.}
#' \item{gelman_diag_list}{ Since version 0.3.3, we incorporated Gelman And Rubin's Convergence Diagnostic using \code{\link[coda]{gelman.diag}}. We included both the statistics and their upper C.I. limits. Values substantially above 1 indicate lack of convergence. Error messages would be printed as they are. This would be only valid for multiple chains (e.g. nchain > 1). Note that for time dependent processes, we only compute the convergence of the last observation data (\code{T_prime}), though it shows to be \code{T_prime+1}, which is due to the day 0 for initialization.}
#' @examples
#' \donttest{
#' NI_complete <- c(
#' 41, 41, 41, 45, 62, 131, 200, 270, 375, 444, 549, 729,
#' 1052, 1423, 2714, 3554, 4903, 5806, 7153, 9074, 11177,
#' 13522, 16678, 19665, 22112, 24953, 27100, 29631, 31728, 33366
#' )
#' RI_complete <- c(
#' 1, 1, 7, 10, 14, 20, 25, 31, 34, 45, 55, 71, 94, 121, 152, 213,
#' 252, 345, 417, 561, 650, 811, 1017, 1261, 1485, 1917, 2260,
#' 2725, 3284, 3754
#' )
#' N <- 58.5e6
#' R <- RI_complete / N
#' Y <- NI_complete / N - R # Jan13->Feb 11
#' ### Step function of pi(t)
#' change_time <- c("01/23/2020", "02/04/2020", "02/08/2020")
#' pi0 <- c(1.0, 0.9, 0.5, 0.1)
#' res.step <- tvt.eSIR(Y, R,
#' begin_str = "01/13/2020", death_in_R = 0.4,
#' T_fin = 200, pi0 = pi0, change_time = change_time, dic = TRUE,
#' casename = "Hubei_step", save_files = TRUE,
#' save_mcmc = FALSE, M = 5e2, nburnin = 2e2
#' )
#' res.step$plot_infection
#' res.step$plot_removed
#' res.step$dic_val
#'
#' ### continuous exponential function of pi(t)
#' res.exp <- tvt.eSIR(Y, R,
#' begin_str = "01/13/2020", death_in_R = 0.4,
#' T_fin = 200, exponential = TRUE, dic = FALSE, lambda0 = 0.05,
#' casename = "Hubei_exp", save_files = FALSE, save_mcmc = FALSE,
#' M = 5e2, nburnin = 2e2
#' )
#' res.exp$plot_infection
#' # res.exp$plot_removed
#'
#' ### without pi(t), the standard state-space SIR model without intervention
#' res.nopi <- tvt.eSIR(Y, R,
#' begin_str = "01/13/2020", death_in_R = 0.4,
#' T_fin = 200, casename = "Hubei_nopi", save_files = FALSE,
#' M = 5e2, nburnin = 2e2
#' )
#' res.nopi$plot_infection
#' # res.nopi$plot_removed
#' }
#'
#' change_time <- c("01/18/2020")
#' pi0<- c(1.0, 0.9)
#' NI_complete2 <- c(41, 45, 62, 131)
#' RI_complete2 <- c(1, 1, 7, 10)
#' N2 <- 1E3
#' res1 <- tvt.eSIR(
#' RI_complete2 / N2,
#' (NI_complete2 - RI_complete2) / N2,
#' begin_str = "01/10/2020",
#' T_fin =10,
#' pi0 = pi0,
#' change_time = change_time,
#' dic = FALSE,
#' casename = "Hubei_step",
#' save_files = FALSE,
#' save_mcmc = FALSE,
#' save_plot_data = FALSE,
#' M = 50,
#' nburnin = 1
#' )
#' closeAllConnections()
#'
#' @export
tvt.eSIR <- function(Y,
R,
pi0 = NULL,
change_time = NULL,
exponential = FALSE,
lambda0 = NULL,
begin_str = "01/13/2020",
T_fin = 200,
nchain = 4,
nadapt = 1e4,
M = 5e2,
thn = 10,
nburnin = 2e2,
dic = FALSE,
death_in_R = 0.02,
beta0 = 0.2586,
gamma0 = 0.0821,
R0 = beta0 / gamma0,
gamma0_sd = 0.1,
R0_sd = 1,
casename = "tvt.eSIR",
file_add = character(0),
add_death = FALSE,
save_files = FALSE,
save_mcmc = FALSE,
save_plot_data = FALSE,
eps = 1e-10) {
oldpar <- par(no.readonly = TRUE)
on.exit(par(oldpar))
beta0 <- R0 * gamma0
len <- round(M / thn) * nchain # number of MCMC draws in total
T_prime <- length(Y)
if (T_prime != length(R)) stop("Y and R should be matched.")
if(T_prime >= T_fin) {
warning("T_fin must be larger than input time series Y and R,",
"we will automatically add 100 time units to it: T_fin = length(Y) + 10.")
T_fin = T_prime + 10
}
Y <- pmax(Y, eps)
R <- pmax(R, eps)
if (add_death == TRUE && death_in_R == 0.02) {
message("use the default death_in_R which is equal to 0.02",
" to plot the death curve in the removed process forecast plot")}
begin <- chron(dates. = begin_str)
chron_ls <- chron(begin:(begin + T_fin))
end <- chron(begin:(begin + T_fin))[T_fin]
message(
"The follow-up is from ", begin, " to ", end,
" and the last observed date is ", chron_ls[T_prime], "."
)
# current data up to this date
gamma_var <- gamma0_sd^2
lognorm_gamma_parm <- lognorm.parm(gamma0, gamma_var)
R0_var <- R0_sd^2
lognorm_R0_parm <- lognorm.parm(R0, R0_var)
if (exponential == FALSE & !is.null(change_time) & !is.null(pi0)) {
message("Running for step-function pi(t)")
if (length(change_time) != length(pi0) - 1) {
stop("We need the length of vector change_time",
" to be the length of pi0 minus 1.")
}
change_time_chorn <- c(begin - 1, chron(dates. = change_time), end)
pi <- rep(pi0, diff(change_time_chorn))
} else if (exponential == TRUE & !is.null(lambda0)) {
message("Running for exponential-function pi(t)")
pi <- exp(-lambda0 * (0:(T_fin - 1)))
} else {
message("Running without pi(t)")
pi <- rep(1, T_fin)
}
################ MCMC ##########
model1.string <- paste0("
model{
for(t in 2:(T_prime+1)){
Km[t-1,1] <- -beta*pi[t-1]*theta[t-1,1]*theta[t-1,2]
Km[t-1,9] <- gamma*theta[t-1,2]
Km[t-1,5] <- -Km[t-1,1]-Km[t-1,9]
Km[t-1,2] <- -beta*pi[t-1]*(theta[t-1,1]+0.5*Km[t-1,1])*(theta[t-1,2]+0.5*Km[t-1,5])
Km[t-1,10] <- gamma*(theta[t-1,2]+0.5*Km[t-1,5])
Km[t-1,6] <- -Km[t-1,2]-Km[t-1,10]
Km[t-1,3] <- -beta*pi[t-1]*(theta[t-1,1]+0.5*Km[t-1,2])*(theta[t-1,2]+0.5*Km[t-1,6])
Km[t-1,11] <- gamma*(theta[t-1,2]+0.5*Km[t-1,6])
Km[t-1,7] <- -Km[t-1,3]-Km[t-1,11]
Km[t-1,4] <- -beta*pi[t-1]*(theta[t-1,1]+Km[t-1,3])*(theta[t-1,2]+Km[t-1,7])
Km[t-1,12] <- gamma*(theta[t-1,2]+Km[t-1,7])
Km[t-1,8] <- -Km[t-1,4]-Km[t-1,12]
alpha[t-1,1] <- theta[t-1,1]+(Km[t-1,1]+2*Km[t-1,2]+2*Km[t-1,3]+Km[t-1,4])/6
alpha[t-1,2] <- theta[t-1,2]+(Km[t-1,5]+2*Km[t-1,6]+2*Km[t-1,7]+Km[t-1,8])/6
alpha[t-1,3] <- theta[t-1,3]+(Km[t-1,9]+2*Km[t-1,10]+2*Km[t-1,11]+Km[t-1,12])/6
theta[t,1:3] ~ ddirch(k*alpha[t-1,1:3])
Y[t-1] ~ dbeta(lambdaY*theta[t,2],lambdaY*(1-theta[t,2]))
R[t-1] ~ dbeta(lambdaR*theta[t,3],lambdaR*(1-theta[t,3]))
}
theta0[1:3]<-c(", 1 - Y[1] - R[1], ",", Y[1], ",", R[1], ")
theta[1,1:3] ~ ddirch(theta0[1:3])
gamma ~ dlnorm(", lognorm_gamma_parm$mu, ",", 1 / lognorm_gamma_parm$var, ")
R0 ~ dlnorm(", lognorm_R0_parm$mu, ",", 1 / lognorm_R0_parm$var, ")
beta <- R0*gamma
k ~ dgamma(2,0.0001)
lambdaY ~ dgamma(2,0.0001)
lambdaR ~ dgamma(2,0.0001)
}
")
model.spec <- textConnection(model1.string)
posterior <- suppressWarnings(
jags.model(
model.spec,
data = list(
"Y" = Y,
"R" = R,
"T_prime" = T_prime,
"pi" = pi
),
n.chains = nchain,
n.adapt = nadapt
)
)
suppressWarnings(
update(posterior, nburnin)
) # burn-in
jags_sample <- suppressWarnings(
jags.samples(
posterior,
c("theta",
"gamma",
"R0",
"beta",
"Y",
"lambdaY",
"lambdaR",
"k"),
n.iter = M,
thin = thn
)
)
if (dic) {
dic_val <- suppressWarnings(
dic.samples(posterior, n.iter = M, thin = thn)
)
} else {
dic_val <- NULL
}
if (save_files) {
png(paste0(file_add, casename, "theta_p.png"), width = 700, height = 900)
plot(as.mcmc.list(jags_sample$theta)[, (1:3) * (T_prime + 1)])
# posterior true probabilities
dev.off()
png(paste0(file_add, casename, "R0_p.png"), width = 700, height = 350)
R0_p <- unlist(as.mcmc.list(jags_sample$R0))
plot(as.mcmc.list(jags_sample$R0))
dev.off()
png(paste0(file_add, casename, "gamma_p.png"), width = 700, height = 350)
gamma_p <- unlist(as.mcmc.list(jags_sample$gamma))
plot(as.mcmc.list(jags_sample$gamma))
dev.off()
png(paste0(file_add, casename, "beta_p.png"), width = 700, height = 350)
beta_p <- unlist(as.mcmc.list(jags_sample$beta))
plot(as.mcmc.list(jags_sample$beta))
dev.off()
png(paste0(file_add, casename, "lambdaY_p.png"), width = 700, height = 350)
lambdaY_p <- unlist(as.mcmc.list(jags_sample$lambdaY))
plot(as.mcmc.list(jags_sample$lambdaY))
dev.off()
png(paste0(file_add, casename, "lambdaR_p.png"), width = 700, height = 350)
lambdaR_p <- unlist(as.mcmc.list(jags_sample$lambdaR))
plot(as.mcmc.list(jags_sample$lambdaR))
dev.off()
png(paste0(file_add, casename, "k_p.png"), width = 700, height = 350)
k_p <- unlist(as.mcmc.list(jags_sample$k))
plot(as.mcmc.list(jags_sample$k))
dev.off()
} else {
R0_p <- unlist(as.mcmc.list(jags_sample$R0))
gamma_p <- unlist(as.mcmc.list(jags_sample$gamma))
beta_p <- unlist(as.mcmc.list(jags_sample$beta))
lambdaY_p <- unlist(as.mcmc.list(jags_sample$lambdaY))
lambdaR_p <- unlist(as.mcmc.list(jags_sample$lambdaR))
k_p <- unlist(as.mcmc.list(jags_sample$k))
}
theta_p <- array(
Reduce(
rbind, as.mcmc.list(jags_sample$theta)
),
dim = c(len, T_prime + 1, 3)
)
theta_p_mean <- apply(theta_p[, T_prime + 1, ], 2, mean)
theta_p_ci <- as.vector(
apply(
theta_p[, T_prime + 1, ],
2,
quantile,
c(0.025, 0.5, 0.975)
)
)
R0_p_mean <- mean(R0_p)
R0_p_ci <- quantile(
R0_p,
c(0.025, 0.5, 0.975)
)
gamma_p_mean <- mean(gamma_p)
gamma_p_ci <- quantile(
gamma_p,
c(0.025, 0.5, 0.975)
)
beta_p_mean <- mean(beta_p)
beta_p_ci <- quantile(
beta_p,
c(0.025, 0.5, 0.975)
)
lambdaY_p_mean <- mean(lambdaY_p)
lambdaY_p_ci <- quantile(
lambdaY_p,
c(0.025, 0.5, 0.975)
)
lambdaR_p_mean <- mean(lambdaR_p)
lambdaR_p_ci <- quantile(
lambdaR_p,
c(0.025, 0.5, 0.975)
)
k_p_mean <- mean(k_p)
k_p_ci <- quantile(
k_p,
c(0.025, 0.5, 0.975)
)
#### Forecast ####
theta_pp <- array(0, dim = c(len, T_fin - T_prime, 3))
Y_pp <- matrix(NA, nrow = len, ncol = T_fin - T_prime)
R_pp <- matrix(NA, nrow = len, ncol = T_fin - T_prime)
for (l in 1:len) {
thetalt1 <- theta_p[l, T_prime + 1, 1]
thetalt2 <- theta_p[l, T_prime + 1, 2]
thetalt3 <- theta_p[l, T_prime + 1, 3]
betal <- c(beta_p)[l]
gammal <- c(gamma_p)[l]
kt <- c(k_p)[l]
lambdaYl <- c(lambdaY_p)[l]
lambdaRl <- c(lambdaR_p)[l]
if (betal < 0 | gammal < 0 | thetalt1 < 0 |
thetalt2 < 0 | thetalt3 < 0) next
for (t in 1:(T_fin - T_prime)) {
Km <- NULL
alpha_pp <- NULL
Km[1] <- -betal * pi[t + T_prime] * thetalt1 * thetalt2
Km[9] <- gammal * thetalt2
Km[5] <- -Km[1] - Km[9]
Km[2] <- -betal *
pi[t + T_prime] * (thetalt1 + 0.5 * Km[1]) * (thetalt2 + 0.5 * Km[5])
Km[10] <- gammal * (thetalt2 + 0.5 * Km[5])
Km[6] <- -Km[2] - Km[10]
Km[3] <- -betal *
pi[t + T_prime] * (thetalt1 + 0.5 * Km[2]) * (thetalt2 + 0.5 * Km[6])
Km[11] <- gammal * (thetalt2 + 0.5 * Km[6])
Km[7] <- -Km[3] - Km[11]
Km[4] <- - betal *
pi[t + T_prime] * (thetalt1 + Km[3]) * (thetalt2 + Km[7])
Km[12] <- gammal * (thetalt2 + Km[7])
Km[8] <- -Km[4] - Km[12]
alpha_pp[1] <- thetalt1 + (Km[1] + 2 * Km[2] + 2 * Km[3] + Km[4]) / 6
alpha_pp[2] <- thetalt2 + (Km[5] + 2 * Km[6] + 2 * Km[7] + Km[8]) / 6
alpha_pp[3] <- thetalt3 + (Km[9] + 2 * Km[10] + 2 * Km[11] + Km[12]) / 6
thetalt_tmp <- rdirichlet(1, kt * c(alpha_pp))
thetalt1 <- theta_pp[l, t, 1] <- thetalt_tmp[1]
thetalt2 <- theta_pp[l, t, 2] <- thetalt_tmp[2]
thetalt3 <- theta_pp[l, t, 3] <- thetalt_tmp[3]
Y_pp[l, t] <- rbeta(1, lambdaYl * thetalt2, lambdaYl * (1 - thetalt2))
R_pp[l, t] <- rbeta(1, lambdaRl * thetalt3, lambdaRl * (1 - thetalt3))
}
}
par(mfrow = c(1, 1))
col2 <- gg_color_hue(2)
Y_band <- data.frame(
t(apply(Y_pp,
2,
quantile,
probs = c(0.025, 0.975),
na.rm = TRUE
))
)
thetaI_band <- data.frame(
t(apply(theta_p[, -1, 2],
2,
quantile,
probs = c(0.025, 0.975),
na.rm = TRUE
))
)
Y_mean <- c(colMeans(Y_pp, na.rm = TRUE))
thetaI_mean <- c(
colMeans(theta_p[, -1, 2], na.rm = TRUE),
colMeans(theta_pp[, , 2], na.rm = TRUE)
)
thetaI_median <- c(
apply(theta_p[, -1, 2], 2, median, na.rm = TRUE),
apply(theta_pp[, , 2], 2, median, na.rm = TRUE)
)
colnames(Y_band) <- c("lower", "upper")
colnames(thetaI_band) <- c("lower", "upper")
data_pre <- data.frame(time = 1:T_prime, Y)
data_post <- data.frame(time = 1:T_prime, thetaI_band)
data_fore <- data.frame(time = (T_prime + 1):T_fin, Y_band, Y_mean)
data_comp <- data.frame(
time = 1:T_fin,
rbind(thetaI_band, Y_band),
phase = c(
rep("pre", nrow(thetaI_band)),
rep("post", nrow(Y_band))
),
mean = thetaI_mean,
median = thetaI_median
)
data_poly <- data.frame(
y = c(thetaI_band$upper,
rev(thetaI_band$lower),
Y_band$upper,
rev(Y_band$lower)),
x = c(1:T_prime,
T_prime:1,
(T_prime + 1):T_fin,
T_fin:(T_prime + 1)),
phase = c(rep("pre", T_prime * 2),
rep("post", (T_fin - T_prime) * 2)),
value = c(rep(col2[1], T_prime * 2),
rep(col2[2], (T_fin - T_prime) * 2))
)
## First-order derivative check
thetaS_mat <- cbind(theta_p[, -1, 1], theta_pp[,, 1])
thetaI_mat <- cbind(theta_p[, -1, 2], theta_pp[,, 2])
thetaR_mat <- cbind(theta_p[, -1, 3], theta_pp[,, 3])
dthetaI_mat_post <- (theta_pp[,, 1] * theta_pp[,, 2]) *
((c(beta_p)) %o% (pi[(T_prime + 1):T_fin])) -
theta_pp[,, 2] * replicate(T_fin - T_prime, c(gamma_p))
dthetaI_mat_pre <- t(
apply(
theta_p[,, 2],
1,
function(v) {
diff(smooth(v))
}
)
)
dthetaI_mat <- cbind(dthetaI_mat_pre, dthetaI_mat_post)
dthetaI <- colMeans(dthetaI_mat, na.rm = TRUE)
dthetaI_tp1 <- (1:T_fin)[which.max(dthetaI)] # first second order derivative=0
dthetaI_tp2 <- (dthetaI_tp1:T_fin)[which.min(dthetaI[dthetaI_tp1:T_fin] > 0)] # first order derivative=0
dthetaI_tp1_rd <- max(round(dthetaI_tp1), 1)
if (dthetaI_tp1_rd > T_prime) {
thetatI_tp1_vec <- thetaI_mat[, dthetaI_tp1_rd]
thetaI_tp1_mean <- mean(thetatI_tp1_vec, na.rm = TRUE)
thetaI_tp1_ci <- quantile(
thetatI_tp1_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
Y_tp1_vec <- Y_pp[, dthetaI_tp1_rd - T_prime]
Y_tp1_mean <- mean(Y_tp1_vec, na.rm = TRUE)
Y_tp1_ci <- quantile(
Y_tp1_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
thetatR_tp1_vec <- thetaR_mat[, dthetaI_tp1_rd]
thetaR_tp1_mean <- mean(thetatR_tp1_vec, na.rm = TRUE)
thetaR_tp1_ci <- quantile(
thetatR_tp1_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
R_tp1_vec <- R_pp[, dthetaI_tp1_rd - T_prime]
R_tp1_mean <- mean(R_tp1_vec, na.rm = TRUE)
R_tp1_ci <- quantile(
R_tp1_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
} else {
thetatI_tp1_vec <- thetaI_mat[, dthetaI_tp1_rd]
thetaI_tp1_mean <- mean(thetatI_tp1_vec, na.rm = TRUE)
thetaI_tp1_ci <- quantile(
thetatI_tp1_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
Y_tp1_vec <- NA
Y_tp1_mean <- mean(Y_tp1_vec, na.rm = TRUE)
Y_tp1_ci <- quantile(
Y_tp1_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
thetatR_tp1_vec <- thetaR_mat[, dthetaI_tp1_rd]
thetaR_tp1_mean <- mean(thetatR_tp1_vec, na.rm = TRUE)
thetaR_tp1_ci <- quantile(
thetatR_tp1_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
R_tp1_vec <- R_pp[, dthetaI_tp1_rd - T_prime]
R_tp1_mean <- mean(R_tp1_vec, na.rm = TRUE)
R_tp1_ci <- quantile(
R_tp1_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
}
dthetaI_tp2_rd <- max(round(dthetaI_tp2), 1)
if (dthetaI_tp1_rd == dthetaI_tp2_rd) {
thetatI_tp2_vec <- NA
thetaI_tp2_mean <- mean(thetatI_tp2_vec, na.rm = TRUE)
thetaI_tp2_ci <- quantile(
thetatI_tp2_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
Y_tp2_vec <- NA
Y_tp2_mean <- mean(Y_tp2_vec, na.rm = TRUE)
Y_tp2_ci <- quantile(
Y_tp2_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
thetatR_tp2_vec <- NA
thetaR_tp2_mean <- mean(thetatR_tp2_vec, na.rm = TRUE)
thetaR_tp2_ci <- quantile(
thetatR_tp2_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
R_tp2_vec <- NA
R_tp2_mean <- mean(R_tp2_vec, na.rm = TRUE)
R_tp2_ci <- quantile(
R_tp2_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
} else if (dthetaI_tp2_rd > T_prime) {
thetatI_tp2_vec <- thetaI_mat[, dthetaI_tp2_rd]
thetaI_tp2_mean <- mean(thetatI_tp2_vec, na.rm = TRUE)
thetaI_tp2_ci <- quantile(
thetatI_tp2_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
Y_tp2_vec <- Y_pp[, dthetaI_tp2_rd - T_prime]
Y_tp2_mean <- mean(Y_tp2_vec, na.rm = TRUE)
Y_tp2_ci <- quantile(
Y_tp2_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
thetatR_tp2_vec <- thetaR_mat[, dthetaI_tp2_rd]
thetaR_tp2_mean <- mean(thetatR_tp2_vec, na.rm = TRUE)
thetaR_tp2_ci <- quantile(
thetatR_tp2_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
R_tp2_vec <- R_pp[, dthetaI_tp2_rd - T_prime]
R_tp2_mean <- mean(R_tp2_vec, na.rm = TRUE)
R_tp2_ci <- quantile(
R_tp2_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
} else {
thetatI_tp2_vec <- thetaI_mat[, dthetaI_tp2_rd]
thetaI_tp2_mean <- mean(thetatI_tp2_vec, na.rm = TRUE)
thetaI_tp2_ci <- quantile(
thetatI_tp2_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
Y_tp2_vec <- NA
Y_tp2_mean <- mean(Y_tp2_vec, na.rm = TRUE)
Y_tp2_ci <- quantile(
Y_tp2_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
thetatR_tp2_vec <- thetaR_mat[, dthetaI_tp2_rd]
thetaR_tp2_mean <- mean(thetatR_tp2_vec, na.rm = TRUE)
thetaR_tp2_ci <- quantile(
thetatR_tp2_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
R_tp2_vec <- NA
R_tp2_mean <- mean(R_tp2_vec, na.rm = TRUE)
R_tp2_ci <- quantile(
R_tp2_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
}
thetaR_max_vec <- thetaR_mat[, T_fin]
thetaR_max_mean <- mean(thetaR_max_vec)
thetaR_max_ci <- quantile(
thetaR_max_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
cumInf_vec <- thetaR_mat[, T_fin] + thetaI_mat[, T_fin]
cumInf_mean <- mean(cumInf_vec)
cumInf_ci <- quantile(
cumInf_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
dthetaI_tp1_date <- chron_ls[dthetaI_tp1]
dthetaI_tp2_date <- chron_ls[dthetaI_tp2]
incidence_vec <- 1 - thetaS_mat[, T_fin]
incidence_mean <- mean(incidence_vec, na.rm = TRUE)
incidence_ci <- quantile(
incidence_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
first_tp_vec <- (1:T_fin)[apply(dthetaI_mat, 1, which.max)] # first second order derivative=0
second_tp_vec <- sapply(
1:len,
function(l) {
(first_tp_vec[l]:T_fin)[
which.min(dthetaI_mat[l, first_tp_vec[l]:T_fin] > 0)
]
}
)
end_p_vec <- sapply(
1:len,
function(l) {
if (any(thetaI_mat[l, first_tp_vec[l]:T_fin] <= eps, na.rm = TRUE)) {
(first_tp_vec[l]:T_fin)[
which.max(thetaI_mat[l, first_tp_vec[l]:T_fin] <= eps)
]
} else {
T_fin
}
}
)
# first order derivative=0
first_tp_mean <- mean(first_tp_vec, na.rm = TRUE)
second_tp_mean <- mean(second_tp_vec, na.rm = TRUE)
end_p_mean <- mean(end_p_vec, na.rm = TRUE)
first_tp_ci <- quantile(
first_tp_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
second_tp_ci <- quantile(
second_tp_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
end_p_ci <- quantile(
end_p_vec,
c(0.025, 0.5, 0.975),
na.rm = TRUE
)
first_tp_date_mean <- chron_ls[first_tp_mean]
second_tp_date_mean <- chron_ls[second_tp_mean]
end_p_date_mean <- chron_ls[end_p_mean]
first_tp_date_ci <- chron_ls[first_tp_ci]
second_tp_date_ci <- chron_ls[second_tp_ci]
end_p_date_ci <- chron_ls[end_p_ci]
names(first_tp_date_ci) <- c("2.5%", "50%", "97.5%")
names(second_tp_date_ci) <- c("2.5%", "50%", "97.5%")
names(end_p_date_ci) <- c("2.5%", "50%", "97.5%")
if (save_files) {
png(paste0(file_add, casename, "deriv.png"), width = 700, height = 350)
plot(
y = dthetaI,
x = chron_ls,
type = "l",
ylab = "1st order derivative",
main = "Infection Prevalence"
)
abline(h = 0, col = 2)
abline(v = chron_ls[T_prime], col = "blue")
if (exponential == FALSE & !is.null(change_time) & !is.null(pi0)) {
abline(
v = change_time_chorn[-c(1, length(change_time_chorn))],
col = "gray"
)
}
legend(
"topright",
legend = c("Last observation", "change point"),
col = c("blue", "gray"),
lty = 1, title = "",
bty = "n"
)
dev.off()
}
## Prepare the Spaghetti plot
sample_dthetaI_mat <- cbind(
dthetaI_mat[sample.int(len, 20, replace = FALSE), ]
)
colnames(sample_dthetaI_mat) <- c(as.character(chron_ls)[-1])
sample_dthetaI_mat_long <- reshape2::melt(sample_dthetaI_mat)
colnames(sample_dthetaI_mat_long) <- c("id", "date", "dthetaI")
sample_dthetaI_mat_long$date <- (chron(as.character(
sample_dthetaI_mat_long$date)))
dthetaI_mean_data <- data.frame(dthetaI, date = chron_ls[-1])
spaghetti_ht <- mean(range(sample_dthetaI_mat)) / 2
spaghetti_plot <- ggplot() +
theme_bw() +
theme(
plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
plot.subtitle = element_text(hjust = 0.5, size = 12),
axis.text.x = element_text(angle = 45, hjust = 1),
axis.text = element_text(size = 15),
axis.title = element_text(size = 14)
) +
geom_line(
data = sample_dthetaI_mat_long,
aes(x = date, y = dthetaI, group = id, color = id)
) +
scale_color_gradientn(colours = rainbow(5, alpha = 0.5)) +
labs(
title = "spaghetti plot of infection prevalence",
x = "time", y = "1st order derivative"
) +
geom_line(
data = dthetaI_mean_data,
aes(x = date, y = dthetaI), color = 1
) +
scale_x_continuous(
labels = as.character(chron_ls)[seq(1, T_fin, 30)],
breaks = as.numeric(chron_ls[-1][seq(1, T_fin, 30)])
) +
annotate(
geom = "text",
label = as.character(chron(chron_ls[T_prime]), format = "mon day"),
x = as.numeric(chron_ls[T_prime]) + 12,
y = spaghetti_ht,
color = "blue"
) +
annotate(
geom = "text",
label = as.character(chron(dthetaI_tp1_date, format = "mon day")),
x = as.numeric(dthetaI_tp1_date) + 12,
y = spaghetti_ht * 1.25,
color = "darkgreen"
) +
geom_vline(
xintercept = as.numeric(chron_ls[T_prime]),
color = "blue",
show.legend = TRUE
) +
geom_vline(
xintercept = as.numeric(dthetaI_tp1_date),
color = "darkgreen",
show.legend = TRUE
)
# if(dthetaI_tp1>T_prime)
spaghetti_plot <- spaghetti_plot +
geom_rect(
data = data.frame(
xmin = as.numeric(first_tp_date_ci[1]),
xmax = as.numeric(first_tp_date_ci[3]),
ymin = -Inf,
ymax = Inf, ci = "first tp"
),
aes(
xmin = xmin,
xmax = xmax,
ymin = ymin,
ymax = ymax
),
fill = "darkgreen",
alpha = 0.15
)
# if(dthetaI_tp2>T_prime)
spaghetti_plot <- spaghetti_plot +
geom_rect(
data = data.frame(
xmin = as.numeric(second_tp_date_ci[1]),
xmax = as.numeric(second_tp_date_ci[3], ci = "second tp"),
ymin = -Inf, ymax = Inf
),
aes(
xmin = xmin,
xmax = xmax,
ymin = ymin,
ymax = ymax
),
fill = "purple",
alpha = 0.15
)
if (dthetaI_tp2_date > dthetaI_tp1_date) {
spaghetti_plot <- spaghetti_plot +
geom_vline(
xintercept = as.numeric(dthetaI_tp2_date),
color = "purple",
show.legend = TRUE
) +
annotate(
geom = "text",
label = as.character(chron(dthetaI_tp2_date, format = "mon day")),
x = as.numeric(dthetaI_tp2_date) + 12,
y = spaghetti_ht * 1.5, color = "purple"
)
}
if (save_files) {
ggsave(paste0(
file_add, casename, "_spaghetti.png"
),
width = 12,
height = 10
)
}
###################
y_text_ht <- max(rbind(thetaI_band, Y_band),
na.rm = TRUE) / 2
plot1 <- ggplot(
data = data_poly,
aes(x = x, y = y)
) +
geom_polygon(
alpha = 0.5,
aes(fill = value, group = phase)
) +
labs(
title = substitute(
paste(
casename, ": infection forecast with prior ",
beta[0], "=", v1, ",",
gamma[0], "=", v2, " and ",
R[0], "=", v3
),
list(
casename = casename,
v1 = format(beta0, digits = 3),
v2 = format(gamma0, digits = 3),
v3 = format(R0, digits = 3)
)
),
subtitle = substitute(
paste(
"Posterior ",
beta[p], "=", v1, ",",
gamma[p], "=", v2, " and ",
R[0], "=", v3
),
list(
v1 = format(beta_p_mean, digits = 3),
v2 = format(gamma_p_mean, digits = 3),
v3 = format(R0_p_mean, digits = 3)
)
),
x = "time", y = "P(Infected)"
) +
geom_line(data = data_comp, aes(x = time, y = median), color = "red") +
geom_vline(xintercept = T_prime, color = "blue", show.legend = TRUE) +
geom_vline(
xintercept = dthetaI_tp1,
color = "darkgreen",
show.legend = TRUE
) +
geom_line(data = data_comp, aes(x = time, y = mean), color = "darkgray") +
geom_point(data = data_pre, aes(x = time, y = Y)) +
theme_bw() +
theme(
plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
plot.subtitle = element_text(hjust = 0.5, size = 12),
axis.text.x = element_text(angle = 45, hjust = 1),
axis.text = element_text(size = 15),
axis.title = element_text(size = 14),
legend.title = element_text(size = 15),
legend.text = element_text(size = 12)
) +
scale_x_continuous(
labels = as.character(chron_ls)[seq(1, T_fin, 30)],
breaks = seq(1, T_fin, 30)
) +
scale_fill_discrete(
name = "Posterior",
labels = c(
expression(paste(y[t[0] + 1:T]^I, " | ", y[1:t[0]]^I, ", ", y[1:t[0]]^R)),
expression(paste(theta[1:t[0]]^I, " | ", y[1:t[0]]^I, ", ", y[1:t[0]]^R))
)
) +
annotate(
geom = "text",
label = as.character(
chron(chron_ls[T_prime]),
format = "mon day"
),
x = T_prime + 12,
y = y_text_ht,
color = "blue"
) +
annotate(
geom = "text",
label = as.character(chron(
dthetaI_tp1_date,
format = "mon day"
)),
x = dthetaI_tp1 + 12,
y = y_text_ht * 1.25,
color = "darkgreen"
)
if (dthetaI_tp2 > dthetaI_tp1) {
plot1 <- plot1 + geom_vline(
xintercept = dthetaI_tp2,
color = "purple",
show.legend = TRUE
) +
annotate(
geom = "text",
label = as.character(chron(dthetaI_tp2_date, format = "mon day")),
x = dthetaI_tp2 + 12, y = y_text_ht * 1.5, color = "purple"
)
}
if (save_files) ggsave(paste0(file_add, casename, "_forecast.png"), width = 12, height = 10)
### Removed
R_band <- data.frame(t(apply(R_pp, 2, quantile, probs = c(0.025, 0.975), na.rm = TRUE)))
thetaR_band <- data.frame(t(apply(theta_p[, -1, 3], 2, quantile, probs = c(0.025, 0.975), na.rm = TRUE)))
R_mean <- c(colMeans(R_pp, na.rm = TRUE))
thetaR_mean <- c(colMeans(theta_p[, -1, 3], na.rm = TRUE), colMeans(theta_pp[, , 3], na.rm = TRUE))
thetaR_med <- c(apply(theta_p[, -1, 3], 2, median, na.rm = TRUE), apply(theta_pp[, , 3], 2, median, na.rm = TRUE))
colnames(R_band) <- c("lower", "upper")
colnames(thetaR_band) <- c("lower", "upper")
data_pre_R <- data.frame(time = 1:T_prime, R) # previous data
data_post_R <- data.frame(time = 1:T_prime, thetaR_band) # posterior of theta^R
data_fore_R <- data.frame(time = (T_prime + 1):T_fin, R_band, R_mean) # The forecast of R after T_prime
data_comp_R <- data.frame(time = 1:T_fin, rbind(thetaR_band, R_band), phase = c(rep("pre", nrow(thetaR_band)), rep("post", nrow(R_band))), mean = thetaR_mean, median = thetaR_med, dead = thetaR_mean * death_in_R, dead_med = thetaR_med * death_in_R) # the filled area--polygon
data_poly_R <- data.frame(y = c(thetaR_band$upper, rev(thetaR_band$lower), R_band$upper, rev(R_band$lower)), x = c(1:T_prime, T_prime:1, (T_prime + 1):T_fin, T_fin:(T_prime + 1)), phase = c(rep("pre", T_prime * 2), rep("post", (T_fin - T_prime) * 2)), value = c(rep(col2[1], T_prime * 2), rep(col2[2], (T_fin - T_prime) * 2)))
r_text_ht <- max(rbind(thetaR_band, R_band), na.rm = TRUE) / 2
plot2 <- ggplot(data = data_poly_R, aes(x = x, y = y)) +
geom_polygon(alpha = 0.5, aes(fill = value, group = phase)) +
labs(
title = substitute(
paste(
casename,
": removed forecast with prior ", beta[0], "=", v1, ",",
gamma[0], "=", v2, " and ", R[0], "=", v3
),
list(
casename = casename, v1 = format(beta0, digits = 3),
v2 = format(gamma0, digits = 3), v3 = format(R0, digits = 3)
)
),
subtitle = substitute(
paste(
"posterior: ", beta[p], "=", v1, ",",
gamma[p], "=", v2, " and ", R[0], "=", v3
),
list(
v1 = format(beta_p_mean, digits = 3), v2 = format(gamma_p_mean, digits = 3),
v3 = format(R0_p_mean, digits = 3)
)
), x = "time", y = "P(Removed)"
) +
geom_line(data = data_comp_R, aes(x = time, y = median), color = "red", linetype = 1) +
geom_vline(xintercept = T_prime, color = "blue") +
geom_vline(xintercept = dthetaI_tp1, color = "darkgreen") +
geom_line(data = data_comp_R, aes(x = time, y = mean), color = "darkgray") +
geom_point(data = data_pre_R, aes(x = time, y = R)) +
theme_bw() +
theme(
plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
plot.subtitle = element_text(hjust = 0.5, size = 12),
axis.text.x = element_text(angle = 45, hjust = 1),
axis.text = element_text(size = 15),
axis.title = element_text(size = 14),
legend.title = element_text(size = 15),
legend.text = element_text(size = 12)
) +
scale_x_continuous(labels = as.character(chron_ls)[seq(1, T_fin, 30)], breaks = seq(1, T_fin, 30)) +
scale_fill_discrete(
name = "Posterior",
labels = c(
expression(paste(y[t[0] + 1:T]^R, " | ", y[1:t[0]]^I, ", ", y[1:t[0]]^R)),
expression(paste(theta[1:t[0]]^R, " | ", y[1:t[0]]^I, ", ", y[1:t[0]]^R))
)
) +
annotate(geom = "text", label = as.character(chron(chron_ls[T_prime]), format = "mon day"), x = T_prime + 12, y = r_text_ht, color = "blue") +
annotate(geom = "text", label = as.character(chron(dthetaI_tp1_date, format = "mon day")), x = dthetaI_tp1 + 12, y = r_text_ht * 1.25, color = "darkgreen")
if (dthetaI_tp2 > dthetaI_tp1) {
plot2 <- plot2 + geom_vline(xintercept = dthetaI_tp2, color = "purple", show.legend = TRUE) + annotate(geom = "text", label = as.character(chron(dthetaI_tp2_date, format = "mon day")), x = dthetaI_tp2 + 12, y = r_text_ht * 1.5, color = "purple")
}
if (add_death) plot2 <- plot2 + geom_line(data = data_comp_R, aes(x = time, y = dead), color = "black", linetype = 1) + geom_line(data = data_comp_R, aes(x = time, y = dead_med), color = "black", linetype = 2)
if (save_files) ggsave(paste0(file_add, casename, "_forecast2.png"), width = 12, height = 10)
out_table1 <- data.frame(matrix(c(theta_p_mean, theta_p_ci, R0_p_mean, R0_p_ci, gamma_p_mean, gamma_p_ci, beta_p_mean, beta_p_ci, incidence_mean = incidence_mean, incidence_ci = incidence_ci, thetaI_tp1_mean = thetaI_tp1_mean, thetaI_tp1_ci = thetaI_tp1_ci, thetaR_tp1_mean = thetaR_tp1_mean, thetaR_tp1_ci = thetaR_tp1_ci, Y_tp1_mean = Y_tp1_mean, Y_tp1_ci = Y_tp1_ci, R_tp1_mean = R_tp1_mean, R_tp1_ci = R_tp1_ci, thetaI_tp2_mean = thetaI_tp2_mean, thetaI_tp2_ci = thetaI_tp2_ci, thetaR_tp2_mean = thetaR_tp2_mean, thetaR_tp2_ci = thetaR_tp2_ci, Y_tp2_mean = Y_tp2_mean, Y_tp2_ci = Y_tp2_ci, R_tp2_mean = R_tp2_mean, R_tp2_ci = R_tp2_ci, thetaR_max_mean, thetaR_max_ci, cumInf_mean = cumInf_mean, cumInf_ci = cumInf_ci), nrow = 1))
out_table2 <- data.frame(matrix(c(dthetaI_tp1_date = as.character(dthetaI_tp1_date), first_tp_mean = as.character(first_tp_date_mean), first_tp_ci = as.character(first_tp_date_ci), dthetaI_tp2_date = as.character(dthetaI_tp2_date), second_tp_mean = as.character(second_tp_date_mean), second_tp_ci = as.character(second_tp_date_ci), end_p_date_mean = as.character(end_p_date_mean), end_p_date_ci = as.character(end_p_date_ci), begin_str = begin_str), nrow = 1))
out_table <- cbind(out_table1, out_table2)
colnames(out_table) <- c(
"thetaS_last_obs_p_mean",
"thetaI_last_obs_p_mean",
"thetaR_last_obs_p_mean",
"thetaS_last_obs_p_ci_low",
"thetaS_last_obs_p_ci_med",
"thetaS_last_obs_p_ci_up",
"thetaI_last_obs_p_ci_low",
"thetaI_last_obs_p_ci_med",
"thetaI_last_obs_p_ci_up",
"thetaR_last_obs_p_ci_low",
"thetaR_last_obs_p_ci_med",
"thetaR_last_obs_p_ci_up",
"R0_p_mean",
"R0_p_ci_low",
"R0_p_ci_med",
"R0_p_ci_up",
"gamma_p_mean",
"gamma_p_ci_low",
"gamma_p_ci_med",
"gamma_p_ci_up",
"beta_p_mean",
"beta_p_ci_low",
"beta_p_ci_med",
"beta_p_ci_up",
"incidence_mean",
"incidence_ci_low",
"incidence_ci_median",
"incidence_ci_up",
"thetaI_tp1_mean",
"thetaI_tp1_ci_low",
"thetaI_tp1_ci_med",
"thetaI_tp1_ci_up",
"thetaR_tp1_mean",
"thetaR_tp1_ci_low",
"thetaR_tp1_ci_med",
"thetaR_tp1_ci_up",
"Y_tp1_mean",
"Y_tp1_ci_low",
"Y_tp1_ci_med",
"Y_tp1_ci_up",
"R_tp1_mean",
"R_tp1_ci_low",
"R_tp1_ci_med",
"R_tp1_ci_up",
"thetaI_tp2_mean",
"thetaI_tp2_ci_low",
"thetaI_tp2_ci_med",
"thetaI_tp2_ci_up",
"thetaR_tp2_mean",
"thetaR_tp2_ci_low",
"thetaR_tp2_ci_med",
"thetaR_tp2_ci_up",
"Y_tp2_mean",
"Y_tp2_ci_low",
"Y_tp2_ci_med",
"Y_tp2_ci_up",
"R_tp2_mean",
"R_tp2_ci_low",
"R_tp2_ci_med",
"R_tp2_ci_up",
"thetaR_max_mean",
"thetaR_max_ci_low",
"thetaR_max_ci_med",
"thetaR_max_ci_up",
"cumInf_mean",
"cumInf_ci_low",
"cumInf_ci_med",
"cumInf_ci_up",
"dthetaI_tp1_date",
"first_tp_mean",
"first_tp_ci_low",
"first_tp_ci_med",
"first_tp_ci_up",
"dthetaI_tp2_date",
"second_tp_mean",
"second_tp_ci_low",
"second_tp_ci_med",
"second_tp_ci_up",
"end_p_mean",
"end_p_ci_low",
"end_p_ci_med",
"end_p_ci_up",
"begin_str"
)
if (save_files) write.csv(out_table, file = paste0(file_add, casename, "_summary.csv"))
if (save_mcmc) {
save(
theta_pp,
Y_pp,
R_pp,
file = paste0(file_add, casename, "_forecast_MCMC.RData")
)
save(
jags_sample,
file = paste0(file_add, casename, "_rjags_MCMC.RData")
)
}
if (save_plot_data) {
other_plot <- list(T_prime = T_prime, T_fin = T_fin, chron_ls = chron_ls, dthetaI_tp1 = dthetaI_tp1, dthetaI_tp2 = dthetaI_tp2, dthetaI_tp1_date = dthetaI_tp1_date, dthetaI_tp2_date = dthetaI_tp2_date, beta_p_mean = beta_p_mean, gamma_p_mean = gamma_p_mean, R0_p_mean = R0_p_mean)
spaghetti_plot_ls <- list(spaghetti_ht = spaghetti_ht, dthetaI_mean_data = dthetaI_mean_data, sample_dthetaI_mat_long = sample_dthetaI_mat_long, first_tp_date_ci = first_tp_date_ci, second_tp_date_ci = second_tp_date_ci)
infection_plot_ls <- list(y_text_ht = y_text_ht, data_poly = data_poly, data_comp = data_comp, data_pre = data_pre)
removed_plot_ls <- list(r_text_ht = r_text_ht, data_poly_R = data_poly_R, data_comp_R = data_comp_R, data_pre_R = data_pre_R)
plot_data_ls <- list(casename = casename, other_plot = other_plot, spaghetti_plot_ls = spaghetti_plot_ls, infection_plot_ls = infection_plot_ls, removed_plot_ls = removed_plot_ls)
save(plot_data_ls, file = paste0(file_add, casename, "_plot_data.RData"))
}
# Gelman-Rubin convergence diagnostics
gelman_out_name <- c("R0", "gamma", "beta", "theta", "lambdaY", "lambdaR", "k")
gelman_diag_list <- lapply(gelman_out_name, function(name) {
if (name == "theta") {
tryCatch(
gelman.diag(as.mcmc.list(jags_sample[[name]])[, (1:3) * (T_prime + 1)],
autoburnin = FALSE
),
error = function(e) e
)
} else {
tryCatch(
gelman.diag(as.mcmc.list(jags_sample[[name]]),
autoburnin = FALSE
),
error = function(e) e
)
}
})
if (save_files) {
sink(paste0(file_add, casename, "_Gelman_diag.txt"))
print(gelman_diag_list)
sink()
}
res <- list(
casename = casename,
incidence_mean = incidence_mean,
incidence_ci = incidence_ci,
out_table = out_table,
plot_infection = plot1,
plot_removed = plot2,
spaghetti_plot = spaghetti_plot,
first_tp_mean = as.character(first_tp_date_mean),
first_tp_ci = as.character(first_tp_date_ci),
second_tp_mean = as.character(second_tp_date_mean),
second_tp_ci = as.character(second_tp_date_ci),
dic_val = dic_val,
gelman_diag_list = gelman_diag_list
)
res
}
if (FALSE) {
NI_complete <- c(
41, 41, 41, 45, 62, 131, 200, 270, 375, 444, 549, 729,
1052, 1423, 2714, 3554, 4903, 5806, 7153, 9074, 11177,
13522, 16678, 19665, 22112, 24953, 27100, 29631, 31728, 33366
)
RI_complete <- c(
1, 1, 7, 10, 14, 20, 25, 31, 34, 45, 55, 71, 94, 121, 152, 213,
252, 345, 417, 561, 650, 811, 1017, 1261, 1485, 1917, 2260,
2725, 3284, 3754
)
N <- 58.5e6
R <- RI_complete / N
Y <- NI_complete / N - R # Jan13->Feb 11
### Step function of pi(t)
change_time <- c("01/23/2020", "02/04/2020", "02/08/2020")
pi0 <- c(1.0, 0.9, 0.5, 0.1)
res.step <- tvt.eSIR(Y, R,
begin_str = "01/13/2020", death_in_R = 0.4, T_fin = 200,
pi0 = pi0, change_time = change_time, dic = FALSE, casename = "Hubei_step", save_files = FALSE, save_mcmc = FALSE,
M = 5e2, nburnin = 2e2
)
res.step$plot_infection
res.step$plot_removed
# res.step$dic_val
res.step$gelman_diag_list
### continuous exponential function of pi(t)
res.exp <- tvt.eSIR(Y, R,
begin_str = "01/13/2020", death_in_R = 0.4, T_fin = 200, exponential = TRUE, dic = FALSE, lambda0 = 0.05,
casename = "Hubei_exp", save_files = TRUE, save_mcmc = FALSE,
M = 5e2, nburnin = 2e2
)
res.exp$plot_infection
# res.exp$plot_removed
### without pi(t), the standard state-space SIR model without intervention
res.nopi <- tvt.eSIR(Y, R,
begin_str = "01/13/2020", death_in_R = 0.4, T_fin = 200,
casename = "Hubei_nopi", save_files = FALSE,
M = 5e2, nburnin = 2e2
)
res.nopi$plot_infection
# res.nopi$plot_removed
}
Try the eSIR package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.