eSAIR {eSIR}

Extended state-space SIR with a subset of the population showing antibody positivity

Description

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.

Usage

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
)

Arguments

Y

the time series of daily observed infected compartment proportions.

R

the time series of daily observed removed compartment proportions, including death and recovered.

alpha0

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.

change_time

the change points over time corresponding to alpha0, to formulate the dirac delta function \alpha_t; its defalt value is NULL.

begin_str

the character of starting time, the default is "01/13/2020".

T_fin

the end of follow-up time after the beginning date begin_str, the default is 200.

nchain

the number of MCMC chains generated by rjags, the default is 4.

nadapt

the iteration number of adaptation in the MCMC. We recommend using at least the default value 1e4 to obtained fully adapted chains.

M

the number of draws in each chain, with no thinning. The default is M=5e2 but suggest using 5e5.

thn

the thinning interval between mixing. The total number of draws thus would become round(M/thn)*nchain. The default is 10.

nburnin

the burn-in period. The default is 2e2 but suggest 2e5.

dic

logical, whether compute the DIC (deviance information criterion) for model selection.

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.

casename

the string of the job's name. The default is "eSAIR".

beta0

the hyperparameter of average transmission rate, the default is the one estimated from the SARS first-month outbreak (0.2586).

gamma0

the hyperparameter of average removed rate, the default is the one estimated from the SARS first-month outbreak (0.0821).

R0

the hyperparameter of the mean reproduction number R0. The default is thus the ratio of beta0/gamma0, which can be specified directly.

gamma0_sd

the standard deviation for the prior distrbution of the removed rate \gamma, the default is 0.1.

R0_sd

the standard deviation for the prior disbution of R0, the default is 1.

file_add

the string to denote the location of saving output files and tables.

add_death

logical, whether add the approximate death curve to the plot, default is false.

save_files

logical, whether to save plots to file.

save_mcmc

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.

save_plot_data

logical, whether save the plotting data or not.

eps

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

Value

casename

the predefined casename.

incidence_mean

mean cumulative incidence, the mean prevalence of cumulative confirmed cases at the end of the study.

incidence_ci

2.5%, 50%, and 97.5% quantiles of the incidences.

out_table

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_infection

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_removed

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.

spaghetti_plot

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.

first_tp_mean

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.

first_tp_mean

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.

first_tp_ci

fwith first_tp_mean, it reports the corresponding credible interval and median.

second_tp_mean

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.

second_tp_ci

with second_tp_mean, it reports the corresponding credible interval and median.

dic_val

the output of dic.samples() in dic.samples, computing deviance information criterion for model comparison.

gelman_diag_list

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.

Examples

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


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()