Title: | Extended State-Space SIR Models |
---|---|
Description: | An implementation of extended state-space SIR models developed by Song Lab at UM school of Public Health. There are several functions available by 1) including a time-varying transmission modifier, 2) adding a time-dependent quarantine compartment, 3) adding a time-dependent antibody-immunization compartment. Wang L. (2020) <doi:10.6339/JDS.202007_18(3).0003>. |
Authors: | Lili Wang [aut], Song Lab [aut] (<http://websites.umich.edu/~songlab/>), Paul Egeler [ctb] (Spectrum Health Office of Research and Education), Michael Kleinsasser [cre] |
Maintainer: | Michael Kleinsasser <[email protected]> |
License: | CC BY 4.0 |
Version: | 0.4.2 |
Built: | 2024-12-07 06:49:01 UTC |
Source: | CRAN |
Confirmed COVID-19 cases in US states
a list with
Province_State name of the US state
date ... a column for each date
Confirmed COVID-19 deaths in US states
a list with
Province_State name of the US state
date ... a column for each date
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 .
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 )
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 )
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 |
change_time |
the change points over time corresponding to |
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 |
nchain |
the number of MCMC chains generated by |
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 |
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 |
gamma0_sd |
the standard deviation for the prior distrbution of the removed rate |
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 ( |
save_plot_data |
logical, whether save the plotting data or not. |
eps |
a non-zero controller so that all the input |
casename |
the predefined |
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 ( |
plot_infection |
plot of summarizing and forecasting for the infection compartment, in which the vertial blue line denotes the last date of data collected ( |
plot_removed |
plot of summarizing and forecasting for the removed compartment with lines similar to those in the |
spaghetti_plot |
20 randomly selected MCMC draws of the first-order derivative of the posterior prevalence of infection, namely |
first_tp_mean |
the date t at which |
first_tp_mean |
the date t at which |
first_tp_ci |
fwith |
second_tp_mean |
the date t at which |
second_tp_ci |
with |
dic_val |
the output of |
gelman_diag_list |
Since version 0.3.3, we incorporated Gelman And Rubin's Convergence Diagnostic using |
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()
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()
Fit an extended state-space SIR model being reduced by in-home hospitalization.
qh.eSIR( Y, R, phi0 = 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 = "qh.eSIR", 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 )
qh.eSIR( Y, R, phi0 = 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 = "qh.eSIR", 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 )
Y |
the time series of daily observed infected compartment proportions. |
R |
the time series of daily observed removed compartment proportions, including death and recovered. |
phi0 |
a vector of values of the dirac delta function |
change_time |
the change points over time corresponding to |
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 |
nchain |
the number of MCMC chains generated by |
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 |
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 "qh.eSIR". |
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 |
gamma0_sd |
the standard deviation for the prior distrbution of the removed rate |
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 ( |
save_plot_data |
logical, whether save the plotting data or not. |
eps |
a non-zero controller so that all the input |
In this function we allow it to characterize time-varying proportions of susceptible due to government-enforced stringent in-home isolation. We expanded the SIR model by adding a quarantine compartment with a time-varying rate of quarantine , the chance of a susceptible person being willing to take in-home isolation at time t.
casename |
the predefined |
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 prevalence processes of the 3 states compartments ( |
plot_infection |
plot of summarizing and forecasting for the infection compartment, in which the vertical blue line denotes the last date of data collected ( |
plot_removed |
plot of summarizing and forecasting for the removed compartment with lines similar to those in the |
spaghetti_plot |
20 randomly selected MCMC draws of the first-order derivative of the posterior prevalence of infection, namely |
first_tp_mean |
the date t at which |
first_tp_mean |
the date t at which |
first_tp_ci |
fwith |
second_tp_mean |
the date t at which |
second_tp_ci |
with |
dic_val |
the output of |
gelman_diag_list |
Since version 0.3.3, we incorporated Gelman And Rubin's Convergence Diagnostic using |
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("01/23/2020", "02/04/2020", "02/08/2020") phi0 <- c(0.1, 0.4, 0.4) res.q <- qh.eSIR(Y, R, begin_str = "01/13/2020", death_in_R = 0.4, phi0 = phi0, change_time = change_time, casename = "Hubei_q", save_files = TRUE, save_mcmc = FALSE, M = 5e2, nburnin = 2e2 ) res.q$plot_infection # res.q$plot_removed res.noq <- qh.eSIR(Y, R, begin_str = "01/13/2020", death_in_R = 0.4, T_fin = 200, casename = "Hubei_noq", M = 5e2, nburnin = 2e2 ) res.noq$plot_infection change_time <- c("01/16/2020") phi0 <- c(0.1) NI_complete2 <- c(41, 45) RI_complete2 <- c(1, 1) N2 <- 1E3 res2 <- qh.eSIR( RI_complete2 / N2, NI_complete2 / N2, begin_str = "01/13/2020", T_fin = 4, phi0 = phi0, change_time = change_time, dic = FALSE, casename = "Hubei_q", save_files = FALSE, save_mcmc = FALSE, save_plot_data = FALSE, M = 50, nburnin = 1 ) closeAllConnections()
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("01/23/2020", "02/04/2020", "02/08/2020") phi0 <- c(0.1, 0.4, 0.4) res.q <- qh.eSIR(Y, R, begin_str = "01/13/2020", death_in_R = 0.4, phi0 = phi0, change_time = change_time, casename = "Hubei_q", save_files = TRUE, save_mcmc = FALSE, M = 5e2, nburnin = 2e2 ) res.q$plot_infection # res.q$plot_removed res.noq <- qh.eSIR(Y, R, begin_str = "01/13/2020", death_in_R = 0.4, T_fin = 200, casename = "Hubei_noq", M = 5e2, nburnin = 2e2 ) res.noq$plot_infection change_time <- c("01/16/2020") phi0 <- c(0.1) NI_complete2 <- c(41, 45) RI_complete2 <- c(1, 1) N2 <- 1E3 res2 <- qh.eSIR( RI_complete2 / N2, NI_complete2 / N2, begin_str = "01/13/2020", T_fin = 4, phi0 = phi0, change_time = change_time, dic = FALSE, casename = "Hubei_q", save_files = FALSE, save_mcmc = FALSE, save_plot_data = FALSE, M = 50, nburnin = 1 ) closeAllConnections()
Confirmed COVID-19 recovered in US states
a list with
Province_State name of the US state
date ... a column for each date
Fit extended state-space SIR model with pre-specified changes in the transmission rate, either stepwise or continuous, accommodating time-varying quarantine protocols.
tvt.eSIR( Y, R, pi0 = NULL, change_time = NULL, exponential = FALSE, lambda0 = 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, 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 )
tvt.eSIR( Y, R, pi0 = NULL, change_time = NULL, exponential = FALSE, lambda0 = 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, 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 )
Y |
the time series of daily observed infected compartment proportions. |
R |
the time series of daily observed removed compartment proportions, including death and recovered. |
pi0 |
the time-dependent transmission rate modifier |
change_time |
the change points over time for step function pi, defalt value is |
exponential |
logical, whether |
lambda0 |
the rate of decline in the exponential survival function in |
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 |
nchain |
the number of MCMC chains generated by |
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 |
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. |
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 |
gamma0_sd |
the standard deviation for the prior distrbution of the removed rate |
R0_sd |
the standard deviation for the prior disbution of R0, the default is 1. |
casename |
the string of the job's name. The default is "tvt.eSIR". |
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 ( |
save_plot_data |
logical, whether save the plotting data or not. |
eps |
a non-zero controller so that all the input |
We fit a state-space model with extended SIR, in which a time-varying transmission rate modifier (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
. The parameters of the function and change points, if any, should be predefined.
casename |
the predefined |
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 ( |
plot_infection |
plot of summarizing and forecasting for the infection compartment, in which the vertial blue line denotes the last date of data collected ( |
plot_removed |
plot of summarizing and forecasting for the removed compartment with lines similar to those in the |
spaghetti_plot |
20 randomly selected MCMC draws of the first-order derivative of the posterior prevalence of infection, namely |
first_tp_mean |
the date t at which |
first_tp_ci |
fwith |
second_tp_mean |
the date t at which |
second_tp_ci |
with |
dic_val |
the output of |
gelman_diag_list |
Since version 0.3.3, we incorporated Gelman And Rubin's Convergence Diagnostic using |
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()
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()
Fit extended state-space SIR model with pre-specified changes in the transmission rate, either stepwise or continuous, accommodating time-varying quarantine protocols.
tvt.eSIR2( Y, R, pi0 = NULL, change_time = NULL, exponential = FALSE, lambda0 = 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, beta0 = 0.2586, gamma0 = 0.0821, R0 = beta0/gamma0, gamma0_sd = 0.1, R0_sd = 1, casename = "tvt.eSIR2", file_add = character(0), add_death = FALSE, save_files = FALSE, save_mcmc = FALSE, save_plot_data = FALSE, eps = 1e-10, time_unit = 1 )
tvt.eSIR2( Y, R, pi0 = NULL, change_time = NULL, exponential = FALSE, lambda0 = 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, beta0 = 0.2586, gamma0 = 0.0821, R0 = beta0/gamma0, gamma0_sd = 0.1, R0_sd = 1, casename = "tvt.eSIR2", file_add = character(0), add_death = FALSE, save_files = FALSE, save_mcmc = FALSE, save_plot_data = FALSE, eps = 1e-10, time_unit = 1 )
Y |
the time series of daily observed infected compartment proportions. |
R |
the time series of daily observed removed compartment proportions, including death and recovered. |
pi0 |
the time-dependent transmission rate modifier |
change_time |
the change points over time for step function pi, defalt
value is |
exponential |
logical, whether |
lambda0 |
the rate of decline in the exponential survival function in
|
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
|
nchain |
the number of MCMC chains generated by
|
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 |
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. |
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 |
gamma0_sd |
the standard deviation for the prior distrbution of the
removed rate |
R0_sd |
the standard deviation for the prior disbution of R0, the default is 1. |
casename |
the string of the job's name. The default is "tvt.eSIR2". |
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 ( |
save_plot_data |
logical, whether save the plotting data or not. |
eps |
a non-zero controller so that all the input |
time_unit |
numeric, newly added argument, which can be changed to an integer (ceiling) of the input to indicate the time unit for each data point. The default is one-day, i.e., we let each input time series data correspond to each day. If this value is set to be 7, each data point represents one-week's aggregated data. |
We fit a state-space model with extended SIR, in which a time-varying
transmission rate modifier (between 0 and 1) is introduced to
the 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
. The parameters of the function and
change points, if any, should be predefined.
casename |
the predefined |
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
prevalence processes of the 3 states compartments ( |
plot_infection |
plot of summarizing and forecasting for the infection
compartment, in which the vertical blue line denotes the last date of data
collected ( |
plot_removed |
plot of summarizing and forecasting for the removed
compartment with lines similar to those in the |
spaghetti_plot |
20 randomly selected MCMC draws of the first-order
derivative of the posterior prevalence of infection, namely
|
first_tp_mean |
the date t at which |
first_tp_ci |
fwith |
second_tp_mean |
the date t at which |
second_tp_ci |
with |
dic_val |
the output of |
gelman_diag_list |
Since version 0.3.3, we incorporated Gelman And
Rubin's Convergence Diagnostic using |
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( Y = RI_complete2 / N2, R = (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()
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( Y = RI_complete2 / N2, R = (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()
Data frame with the populations of each state.
a list with
Province_State State as a character string
N Population of the state