In this function we allow it to characterize time-varying immunization among a subset of the population that have been tested positive in an antibody assessment. We expanded the SIR model by adding a time-varying antibody-positive proportion \(\alpha_t\).
eSAIR(
Y,
R,
alpha0 = NULL,
change_time = NULL,
begin_str = "01/13/2020",
T_fin = 200,
nchain = 4,
nadapt = 10000,
M = 500,
thn = 10,
nburnin = 200,
dic = FALSE,
death_in_R = 0.02,
casename = "eSAIR",
beta0 = 0.2586,
gamma0 = 0.0821,
R0 = beta0/gamma0,
gamma0_sd = 0.1,
R0_sd = 1,
file_add = character(0),
add_death = FALSE,
save_files = FALSE,
save_mcmc = FALSE,
save_plot_data = FALSE,
eps = 1e-10
)the time series of daily observed infected compartment proportions.
the time series of daily observed removed compartment proportions, including death and recovered.
a vector of values of the dirac delta function \(\alpha_t\). Each entry denotes the proportion that will be immunized at each change time point. Note that all the entries lie between 0 and 1, its default is NULL.
the change points over time corresponding to alpha0, to formulate the dirac delta function \(\alpha_t\); its defalt value is NULL.
the character of starting time, the default is "01/13/2020".
the end of follow-up time after the beginning date begin_str, the default is 200.
the number of MCMC chains generated by rjags, the default is 4.
the iteration number of adaptation in the MCMC. We recommend using at least the default value 1e4 to obtained fully adapted chains.
the number of draws in each chain, with no thinning. The default is M=5e2 but suggest using 5e5.
the thinning interval between mixing. The total number of draws thus would become round(M/thn)*nchain. The default is 10.
the burn-in period. The default is 2e2 but suggest 2e5.
logical, whether compute the DIC (deviance information criterion) for model selection.
the numeric value of average of cumulative deaths in the removed compartments. The default is 0.4 within Hubei and 0.02 outside Hubei.
the string of the job's name. The default is "eSAIR".
the hyperparameter of average transmission rate, the default is the one estimated from the SARS first-month outbreak (0.2586).
the hyperparameter of average removed rate, the default is the one estimated from the SARS first-month outbreak (0.0821).
the hyperparameter of the mean reproduction number R0. The default is thus the ratio of beta0/gamma0, which can be specified directly.
the standard deviation for the prior distrbution of the removed rate \(\gamma\), the default is 0.1.
the standard deviation for the prior disbution of R0, the default is 1.
the string to denote the location of saving output files and tables.
logical, whether add the approximate death curve to the plot, default is false.
logical, whether to save plots to file.
logical, whether save (TRUE) all the MCMC outputs or not (FALSE).The output file will be an .RData file named by the \(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 theta_p[,,1] and afterwards as theta_pp[,,1] for \(\theta_t^S\), theta_p[,,2] and theta_pp[,,2] for \(\theta_t^I\), and theta_p[,,3] and theta_pp[,,3] for \(\theta_t^R\). The posterior draws of the prevalence process of the antibody-immunized compartment can be obtained via thetaA_p and thetaA_pp. Moreover, the input and predicted proportions Y, Y_pp, R and R_pp can also be retrieved. The prevalence and prediceted proportion matrices have rows for MCMC replicates, and columns for days. The MCMC posterior draws of other parameters including beta, gamma, R0, and variance controllers k_p, lambdaY_p, lambdaR_p are also available.
logical, whether save the plotting data or not.
a non-zero controller so that all the input Y and R values would be bounded above 0 (at least eps). Its default value is 1e-10
the predefined casename.
mean cumulative incidence, the mean prevalence of cumulative confirmed cases at the end of the study.
2.5%, 50%, and 97.5% quantiles of the incidences.
summary tables including the posterior mean of the prevalance processes of the 3 states compartments (\(\theta_t^S,\theta_t^I,\theta_t^R,\theta_t^H\)) at last date of data collected ((\(t^\prime\)) decided by the lengths of your input data Y and R), and their respective credible inctervals (ci); the respective means and ci's of the reporduction number (R0), removed rate (\(\gamma\)), transmission rate (\(\beta\)).
plot of summarizing and forecasting for the infection compartment, in which the vertial blue line denotes the last date of data collected (\(t^\prime\)), the vertial darkgray line denotes the deacceleration point (first turning point) that the posterior mean first-derivative of infection prevalence \(\dot{\theta}_t^I\) achieves the maximum, the vertical purple line denotes the second turning point that the posterior mean first-derivative infection proportion \(\dot{\theta}_t^I\) equals zero, the darkgray line denotes the posterior mean of the infection prevalence \(\theta_t^I\) and the red line denotes its posterior median.
plot of summarizing and forecasting for the removed compartment with lines similar to those in the plot_infection. The vertical lines are identical, but the horizontal mean and median correspond to the posterior mean and median of the removed process \(\theta_t^R\). An additional line indicates the estimated death prevalence from the input death_in_R.
20 randomly selected MCMC draws of the first-order derivative of the posterior prevalence of infection, namely \(\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 plot_infection and plot_removed. Moreover, the 95% credible intervals of these turning points are also highlighted by semi-transparent rectangles.
the date t at which \(\ddot{\theta}_t^I=0\), calculated as the average of the time points with maximum posterior first-order derivatives \(\dot{\theta}_t^I\); this value may be slightly different from the one labeled by the "darkgreen" lines in the two plots plot_infection and plot_removed, which indicate the stationary point such that the first-order derivative of the averaged posterior of \(\theta_t^I\) reaches its maximum.
the date t at which \(\ddot{\theta}_t^I=0\), calculated as the average of the time points with maximum posterior first-order derivatives \(\dot{\theta}_t^I\); this value may be slightly different from the one labeled by the "darkgreen" lines in the two plots plot_infection and plot_removed, which indicate the stationary point such that the first-order derivative of the averaged posterior of \(\theta_t^I\) reaches its maximum.
fwith first_tp_mean, it reports the corresponding credible interval and median.
the date t at which \(\theta_t^I=0\), calculated as the average of the stationary points of all of posterior first-order derivatives \(\dot{\theta}_t^I\); this value may be slightly different from the one labeled by the "pruple" lines in the plots of plot_infection and plot_removed. The latter indicate stationary t at which the first-order derivative of the averaged posterior of \(\theta_t^I\) equals zero.
with second_tp_mean, it reports the corresponding credible interval and median.
the output of dic.samples() in dic.samples, computing deviance information criterion for model comparison.
Since version 0.3.3, we incorporated Gelman And Rubin's Convergence Diagnostic using 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 (T_prime), though it shows to be T_prime+1, which is due to the day 0 for initialization.
# NOT RUN {
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
change_time <- c("02/08/2020")
alpha0 <- c(0.2) # 20% of the susceptible population were found immunized
res.antibody <- eSAIR(Y, R,
begin_str = "01/13/2020", death_in_R = 0.4,
alpha0 = alpha0, change_time = change_time,
casename = "Hubei_antibody", save_files = TRUE, save_mcmc = FALSE,
M = 5e2, nburnin = 2e2
)
res.antibody$plot_infection
# }
# NOT RUN {
change_time <- c("01/16/2020")
alpha0 <- c(0.2)
NI_complete2 <- c(41, 45)
RI_complete2 <- c(1, 1)
N2 <- 1E3
res3 <- eSAIR(
RI_complete2 / N2,
NI_complete2 / N2,
begin_str = "01/13/2020",
T_fin = 4,
alpha0 = alpha0,
change_time = change_time,
dic = FALSE,
casename = "Hubei_q",
save_files = FALSE,
save_mcmc = FALSE,
save_plot_data = FALSE,
M = 50,
nburnin = 1
)
closeAllConnections()
# }
Run the code above in your browser using DataLab