perfRestim: Routine to measure the performance of estimR and estimRmcmc

Description

This routine can be used to check the 'statistical performance' of the estimR() and estimRmcmc() routines to estimate the reproduction number \(R_t\). It simulates epidemics using the episim() function and computes the Bias, MSE, coverage probability (CP) and width of \(90\%\) and \(95\%\) credible intervals for \(R_t\) averaged over days \(t=8,...,T\), where \(T\) is the total number of days of the simulated epidemics. As such, it can be used to reproduce part of the results in Gressani et al. (2022) Table 1 and Table 2, respectively. Small differences in results are due to a restructuring of the code since version 1.0.6. If strict reproducible results are required, please refer to version 1.0.6 of the EpiLPS package or visit the GitHub repository https://github.com/oswaldogressani/EpiLPS-ArticleCode.

Usage

perfRestim(nsim = 100, scenario = 1, days = 40, K = 40,
 method = c("LPSMAP", "LPSMALA"), mcmciter = 3000, burnin = 1000,
 si = c("flu", "sars", "mers"), seed = 1325, overdisp = 1000)

Value

A list with the following components:

LPS: Results for the LPS approach.
EpiEstim: Results for the EpiEstim approach with weekly sliding windows.
inciplot: The simulated incidence time series.
Rlpsplot: Estimated \(R_t\) trajectories with LPS.
Repiestimplot: Estimated \(R_t\) trajectories with EpiEstim.

Arguments

nsim: Total number of simulated epidemics.
scenario: The scenario to be used in episim().
days: Number of days for the simulated epidemics.
K: Number of B-splines basis function in the P-spline model.
method: The method for LPS, either LPSMAP or LPSMALA.
mcmciter: Number of MCMC samples for method LPSMALA.
burnin: Burn-in for method LPSMALA.
si: The discrete serial interval distribution. Possible specifications are "flu", "sars" or "mers".
seed: A seed for reproducibility.
overdisp: The value of the overdispersion parameter for the negative binomial model in the episim() routine.

Author

Oswaldo Gressani oswaldo_gressani@hotmail.fr

References

Gressani, O., Wallinga, J., Althaus, C. L., Hens, N. and Faes, C. (2022). EpiLPS: A fast and flexible Bayesian tool for estimation of the time-varying reproduction number. Plos Computational Biology, 18(10): e1010618.

Examples

Run this code

# # FLU serial interval (Scenarios 1-4)
# S1 <- perfRestim(si = "flu", scenario = 1, seed = 1325)
# S1mcmc <- perfRestim(si = "flu", scenario = 1, seed = 1325, method = "LPSMALA")
# suppressWarnings(gridExtra::grid.arrange(S1$inciplot, S1$Rlpsplot, S1$Repiestimplot, nrow = 1))
# S2 <- perfRestim(si = "flu", scenario = 2, seed = 1123)
# S2mcmc <- perfRestim(si = "flu", scenario = 2, seed = 1123, method = "LPSMALA")
# suppressWarnings(gridExtra::grid.arrange(S2$inciplot, S2$Rlpsplot, S2$Repiestimplot, nrow = 1))
# S3 <- perfRestim(si = "flu", scenario = 3, seed = 1314)
# S3mcmc <- perfRestim(si = "flu", scenario = 3, seed = 1314, method = "LPSMALA")
# suppressWarnings(gridExtra::grid.arrange(S3$inciplot, S3$Rlpsplot, S3$Repiestimplot, nrow = 1))
# S4 <- perfRestim(si = "flu", scenario = 4, seed = 1966)
# S4mcmc <- perfRestim(si = "flu", scenario = 4, seed = 1966, method = "LPSMALA")
# suppressWarnings(gridExtra::grid.arrange(S4$inciplot, S4$Rlpsplot, S4$Repiestimplot, nrow = 1))
#
# # SARS serial interval (Scenarios 5-8)
# S5 <- perfRestim(si = "sars", scenario = 1, seed = 1998, overdisp = 5)
# S5mcmc <- perfRestim(si = "sars", scenario = 1, seed = 1998, overdisp = 5, method = "LPSMALA")
# suppressWarnings(gridExtra::grid.arrange(S5$inciplot, S5$Rlpsplot, S5$Repiestimplot, nrow = 1))
# S6 <- perfRestim(si = "sars", scenario = 2, seed = 1870, overdisp = 5)
# S6mcmc <- perfRestim(si = "sars", scenario = 2, seed = 1870, overdisp = 5, method = "LPSMALA")
# suppressWarnings(gridExtra::grid.arrange(S6$inciplot, S6$Rlpsplot, S6$Repiestimplot, nrow = 1))
# S7 <- perfRestim(si = "sars", scenario = 3, seed = 115,  overdisp = 5)
# S7mcmc <- perfRestim(si = "sars", scenario = 3, seed = 115,  overdisp = 5, method = "LPSMALA")
# suppressWarnings(gridExtra::grid.arrange(S7$inciplot, S7$Rlpsplot, S7$Repiestimplot, nrow = 1))
# S8 <- perfRestim(si = "sars", scenario = 4, seed = 1464, overdisp = 5)
# S8mcmc <- perfRestim(si = "sars", scenario = 4, seed = 1464, overdisp = 5, method = "LPSMALA")
# suppressWarnings(gridExtra::grid.arrange(S8$inciplot, S8$Rlpsplot, S8$Repiestimplot, nrow = 1))
#
# # MERS serial interval (Scenario 9)
# S9 <- perfRestim(si = "mers", scenario = 5, days = 60,  seed = 1905, overdisp = 50)
# S9mcmc <- perfRestim(si = "mers", scenario = 5, days = 60,
# seed = 1905, overdisp = 50, method = "LPSMALA")
# suppressWarnings(gridExtra::grid.arrange(S9$inciplot, S9$Rlpsplot,
# S9$Repiestimplot, nrow = 1))
#
# #(Partially recovering Table 2 and Table 3 of Gressani et al. 2022)
# simsummary <- matrix(0, nrow = 36, ncol = 7)
# colnames(simsummary) <- c("Method", "Bias", "MSE", "CP90%", "CP95%",
#                           "CIwidth90%", "CIwidth95%")
# simsummary <- as.data.frame(simsummary)
#
# # Scenario 1
# simsummary[1,] <- c(rownames(S1$LPS),S1$LPS)
# simsummary[2,] <- c(rownames(S1mcmc$LPS),S1mcmc$LPS)
# simsummary[3,] <- c(rownames(S1$EpiEstim),S1$EpiEstim)
# simsummary[4,] <- rep("--",7)
# # Scenario 2
# simsummary[5,] <- c(rownames(S2$LPS),S2$LPS)
# simsummary[6,] <- c(rownames(S2mcmc$LPS),S2mcmc$LPS)
# simsummary[7,] <- c(rownames(S2$EpiEstim),S2$EpiEstim)
# simsummary[8,] <- rep("--",7)
# # Scenario 3
# simsummary[9,] <- c(rownames(S3$LPS),S3$LPS)
# simsummary[10,] <- c(rownames(S3mcmc$LPS),S3mcmc$LPS)
# simsummary[11,] <- c(rownames(S3$EpiEstim),S3$EpiEstim)
# simsummary[12,] <- rep("--",7)
# # Scenario 4
# simsummary[13,] <- c(rownames(S4$LPS),S4$LPS)
# simsummary[14,] <- c(rownames(S4mcmc$LPS),S4mcmc$LPS)
# simsummary[15,] <- c(rownames(S4$EpiEstim),S4$EpiEstim)
# simsummary[16,] <- rep("--",7)
# # Scenario 5
# simsummary[17,] <- c(rownames(S5$LPS),S5$LPS)
# simsummary[18,] <- c(rownames(S5mcmc$LPS),S5mcmc$LPS)
# simsummary[19,] <- c(rownames(S5$EpiEstim),S5$EpiEstim)
# simsummary[20,] <- rep("--",7)
# # Scenario 6
# simsummary[21,] <- c(rownames(S6$LPS),S6$LPS)
# simsummary[22,] <- c(rownames(S6mcmc$LPS),S6mcmc$LPS)
# simsummary[23,] <- c(rownames(S6$EpiEstim),S6$EpiEstim)
# simsummary[24,] <- rep("--",7)
# # Scenario 7
# simsummary[25,] <- c(rownames(S7$LPS),S7$LPS)
# simsummary[26,] <- c(rownames(S7mcmc$LPS),S7mcmc$LPS)
# simsummary[27,] <- c(rownames(S7$EpiEstim),S7$EpiEstim)
# simsummary[28,] <- rep("--",7)
# # Scenario 8
# simsummary[29,] <- c(rownames(S8$LPS),S8$LPS)
# simsummary[30,] <- c(rownames(S8mcmc$LPS),S8mcmc$LPS)
# simsummary[31,] <- c(rownames(S8$EpiEstim),S8$EpiEstim)
# simsummary[32,] <- rep("--",7)
# # Scenario 9
# simsummary[33,] <- c(rownames(S9$LPS),S9$LPS)
# simsummary[34,] <- c(rownames(S9mcmc$LPS),S9mcmc$LPS)
# simsummary[35,] <- c(rownames(S9$EpiEstim),S9$EpiEstim)
# simsummary[36,] <- rep("--",7)
# simsummary

Run the code above in your browser using DataLab