Rt Methods
Dhihram Tenrisau
2023-08-07
Introduction
Di sini akan ditunjukkan beberapa perhitungan Reproduksi Efektif atau
(Rt) dengan beberapa package di R. seperti EpiEstim,
EpiInvert, R0, dan EpiNow2
Library
Membuka File
Di metode ini, kita akan menggunakan file yang bersumber dari data COVID-19 di Kementerian Kesehatan Indonesia pada tanggal 23 Juni 2023.
data_nas <- read_excel("C:/Users/dhihr/Downloads/230608_nasional_jb_njb.xlsx")
data_olah_nasional <- filter(data_nas, Provinsi == "Indonesia")
DT::datatable(data_nas)Epiestim
EDA
Rt_Indonesia <- data.frame(I = data_olah_nasional$`Kasus Harian`, dates = data_olah_nasional$Tanggal, Provinsi = data_olah_nasional$Provinsi)
Rt_Indonesia$I <- Rt_Indonesia$I %>% replace_na(0)
Rt_Indonesia$dates <- as.Date(Rt_Indonesia$dates)
DT::datatable(Rt_Indonesia)str(Rt_Indonesia)## 'data.frame': 1076 obs. of 3 variables:
## $ I : num 1193 1095 1287 1254 1710 ...
## $ dates : Date, format: "2020-06-28" "2020-06-29" ...
## $ Provinsi: chr "Indonesia" "Indonesia" "Indonesia" "Indonesia" ...kasus_tanggal_Indonesia <- Rt_Indonesia[-c(1:7),]Utilizing
res_parametric_si_Indonesia <- estimate_R(Rt_Indonesia$I,
method="parametric_si",
config = make_config(list(
mean_si = 3.5,
std_si = 2.8))
)## Default config will estimate R on weekly sliding windows.
## To change this change the t_start and t_end arguments.plot(res_parametric_si_Indonesia)## Warning: The `guide` argument in `scale_*()` cannot be `FALSE`. This was deprecated in
## ggplot2 3.3.4.
## ℹ Please use "none" instead.
## ℹ The deprecated feature was likely used in the incidence package.
## Please report the issue at <https://github.com/reconhub/incidence/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
R_epiestim_Indonesia <- data.frame(tanggal = kasus_tanggal_Indonesia$dates, median_r = round(res_parametric_si_Indonesia$R$`Median(R)`,2),
Provinsi = kasus_tanggal_Indonesia$Provinsi)EpiNow
EDA
epinow1 <- data.frame(date = data_olah_nasional$Tanggal, confirm = data_olah_nasional$`Kasus Harian`)
epinow_template <- tail(epinow1, n = 100)
reported_cases <- epinow_template
reported_cases$date <- as.Date(reported_cases$date, format='%Y-%m-%d')
str(reported_cases)## 'data.frame': 100 obs. of 2 variables:
## $ date : Date, format: "2023-03-01" "2023-03-02" ...
## $ confirm: num 252 278 211 207 165 144 303 312 307 295 ...reporting_delay <- estimate_delay(rlnorm(1000, log(3), 1),
max_value = 15, bootstraps = 1)
generation_time <- get_generation_time(disease = "SARS-CoV-2", source = "ganyani")
incubation_period <- get_incubation_period(disease = "SARS-CoV-2", source = "lauer")
estimates <- epinow(reported_cases = reported_cases,
generation_time = generation_time,
delays = delay_opts(incubation_period, reporting_delay),
rt = rt_opts(prior = list(mean = 2, sd = 0.2)),
stan = stan_opts(cores = 4))## Logging threshold set at INFO for the EpiNow2 logger## Writing EpiNow2 logs to the console and: C:\Users\dhihr\AppData\Local\Temp\RtmpKK6pQE/regional-epinow/2023-06-08.log## Logging threshold set at INFO for the EpiNow2.epinow logger## Writing EpiNow2.epinow logs to the console and: C:\Users\dhihr\AppData\Local\Temp\RtmpKK6pQE/epinow/2023-06-08.log## WARN [2023-10-14 01:21:09] epinow: There were 9 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them. -
## WARN [2023-10-14 01:21:09] epinow: Examine the pairs() plot to diagnose sampling problems
## -knitr::kable(summary(estimates))| measure | estimate |
|---|---|
| New confirmed cases by infection date | 102 (30 – 369) |
| Expected change in daily cases | Likely decreasing |
| Effective reproduction no. | 0.68 (0.33 – 1.3) |
| Rate of growth | -0.093 (-0.21 – 0.084) |
| Doubling/halving time (days) | -7.5 (8.2 – -3.3) |
head(summary(estimates, type = "parameters", params = "R"))## date variable strat type median mean sd lower_90
## 1: 2023-03-06 R <NA> estimate 1.694045 1.699673 0.16081840 1.4392897
## 2: 2023-03-07 R <NA> estimate 1.542666 1.545543 0.11322794 1.3666038
## 3: 2023-03-08 R <NA> estimate 1.406980 1.410477 0.09198956 1.2671839
## 4: 2023-03-09 R <NA> estimate 1.295742 1.295820 0.08747250 1.1505835
## 5: 2023-03-10 R <NA> estimate 1.207527 1.201811 0.08822609 1.0492017
## 6: 2023-03-11 R <NA> estimate 1.133180 1.127810 0.08780674 0.9701602
## lower_50 lower_20 upper_20 upper_50 upper_90
## 1: 1.588868 1.656258 1.735909 1.804089 1.979385
## 2: 1.464988 1.514602 1.569707 1.618907 1.739946
## 3: 1.347918 1.383712 1.430038 1.467035 1.570285
## 4: 1.239400 1.275433 1.317432 1.351645 1.437212
## 5: 1.146851 1.185349 1.226309 1.261037 1.337892
## 6: 1.071065 1.111330 1.156132 1.188944 1.261446head(summary(estimates, output = "estimated_reported_cases"))## date type median mean sd lower_90 lower_50 lower_20
## 1: 2023-03-06 gp_rt 148.0 152.4410 42.12408 93.00 123.00 138
## 2: 2023-03-07 gp_rt 282.0 286.9100 73.95496 178.00 234.75 263
## 3: 2023-03-08 gp_rt 270.5 276.5150 73.04927 168.00 226.00 252
## 4: 2023-03-09 gp_rt 248.0 256.5785 69.57393 155.00 208.00 231
## 5: 2023-03-10 gp_rt 250.0 255.9080 68.75323 155.00 207.00 234
## 6: 2023-03-11 gp_rt 267.0 274.2875 71.06700 170.95 223.00 251
## upper_20 upper_50 upper_90
## 1: 158 177.00 229.05
## 2: 300 332.00 417.05
## 3: 288 320.25 405.05
## 4: 267 301.00 380.00
## 5: 265 301.00 374.05
## 6: 286 317.00 403.00R_epinow1 <- summary(estimates, type = "parameters", params = "R")
R_epinow2 <- summary(estimates, output = "estimated_reported_cases")
R_epinow3 <- summary(estimates, type = "parameters", params = "infections")
plot(estimates)
TD Bettencourt & Ribiero
EDA
data_kasus <- kasus_tanggal_Indonesia$I
names(data_kasus) <- kasus_tanggal_Indonesia$dates
data_kasus <- tail(data_kasus, 1000)
tail(data_kasus)## 2023-06-03 2023-06-04 2023-06-05 2023-06-06 2023-06-07 2023-06-08
## 252 229 231 362 247 254Utilizing
mGT <- generation.time("gamma", c(3.5,2.8))
SB <- est.R0.SB(data_kasus, begin=900, end=1000, mGT)
SB$R## [1] 0.80 1.09 0.73 0.78 0.68 0.66 1.09 1.09 1.07 1.04 1.08 0.95 0.95 1.12 1.09
## [16] 1.09 1.05 1.07 1.00 0.98 1.13 1.11 1.02 1.05 1.08 1.08 1.03 1.10 1.13 1.11
## [31] 1.07 1.08 1.05 1.02 1.11 1.12 1.09 1.09 1.04 1.07 1.05 1.16 1.16 1.15 1.15
## [46] 1.12 1.11 1.08 1.17 1.15 1.13 1.09 1.07 1.02 1.05 1.11 1.14 1.19 1.20 1.19
## [61] 1.13 1.06 1.10 1.21 1.18 1.15 1.12 1.08 1.06 1.11 1.10 1.08 1.08 1.06 1.04
## [76] 1.03 1.02 1.04 1.02 1.01 1.02 1.01 1.01 1.04 1.03 1.01 1.02 1.01 1.00 1.00
## [91] 1.04 1.01 1.01 1.00 1.00 1.00 1.00 1.00 1.00 1.00plot(SB)
SB$end## [1] "2023-06-08"SB$begin## [1] "2023-02-28"data_bettencourt <- tail(data_kasus,101)
data_bettencourt <- as.data.frame(data_bettencourt)
data_bettencourt$date <- row.names(data_bettencourt)
data_bettencourt <- data_bettencourt %>% filter(row_number() <= n()-1)
data_bettencourt$bettencourt <- SB$RTD Teunis & Wallinga
Utilizing
TD <- est.R0.TD(data_kasus, mGT, begin=900, end=1000, nsim=80)## Warning in est.R0.TD(data_kasus, mGT, begin = 900, end = 1000, nsim = 80):
## Accurate confidence interval for R(t) requires a large number of simulations.
## Consider increasing 'nsim'## Warning in est.R0.TD(data_kasus, mGT, begin = 900, end = 1000, nsim = 80):
## Using initial incidence as initial number of cases.plot(TD)
J <- data.frame(R0 = TD$R)
J$tanggal <- row.names(J)
tail(J)## R0 tanggal
## 2023-06-03 0.7984722 2023-06-03
## 2023-06-04 0.8819473 2023-06-04
## 2023-06-05 0.9519575 2023-06-05
## 2023-06-06 0.8241786 2023-06-06
## 2023-06-07 0.8555524 2023-06-07
## 2023-06-08 0.0000000 2023-06-08J$tanggal <- as.Date(J$tanggal)EpiInvert
Utilizing
kasus_tanggal_Indonesia$dates <- as.Date(kasus_tanggal_Indonesia$dates)
EpiInvert <- EpiInvert(kasus_tanggal_Indonesia$I,
"2023-06-25",kasus_tanggal_Indonesia$dates,
select_params(list(mean_si = 3.5,sd_si=2.8, shift_si=-1)))## EpiInvert parameters used:
## Incidence tail : i[1063]=252, i[1064]=229, i[1065]=231, i[1066]=362, i[1067]=247, i[1068]=254,
## Last incidence date 2023-06-25
## Festive days tail : 2023-06-03, 2023-06-04, 2023-06-05, 2023-06-06, 2023-06-07, 2023-06-08,
## max_time_interval=150
## Shifted log-normal serial interval parameters:
## mean_si=3.500000
## sd_si=2.800000
## shift_si=-1.000000
## Rt_regularization_weight=5.000000
## seasonality_regularization_weight=5.000000
## EXECUTION TIME : 1.635000 SECONDSEpiInvert_plot(EpiInvert)
Gabung Seluruh Metode
EDA
R_epinow <- data.frame(Date = as.Date(R_epinow1$date), Epinow = round(R_epinow1$median,2), treshold = 1)
data_ends <- R_epinow %>% filter(Date == "2023-06-24")
R_teunis <- data.frame(Date = as.Date(J$tanggal), Teunis = round(J$R0,2), treshold = 1)
R_teunis <- R_teunis %>% filter(row_number() <= n()-1)
data_ends2 <- R_teunis %>% filter(Date == "2023-06-24")
R_epiestim <-data.frame(Date = as.Date(R_epiestim_Indonesia$tanggal), Epiestim = round(R_epiestim_Indonesia$median_r,2), treshold = 1)
data_ends3 <- R_epiestim %>% filter(Date == "2023-06-24")
R_bettencourt <- data.frame(Date = as.Date(data_bettencourt$date), bettencourt = data_bettencourt$bettencourt, treshold = 1)
data_ends4 <- R_bettencourt %>% filter(Date == "2023-06-24")
Rt_Alvarez <- data.frame(Date = as.Date(EpiInvert$dates), EpiInvert = round(EpiInvert$Rt,2), treshold = 1)
data_ends5 <- Rt_Alvarez %>% filter(Date == "2023-06-24")Plot Gabung
p <- ggplot() + geom_line(data = R_epinow, aes(x = Date, y = Epinow), size = 1.5, color = "turquoise4") +
geom_line(data = R_epinow, aes(x = Date, y = treshold), size = 0.5, color = "black", linetype = "dashed") +
theme(plot.title=element_text(size=16), panel.grid.minor = element_blank()) + theme(strip.text = element_text(face="bold", size=9)) + labs(title ="Rt Epinow") + ylab("Rt") + xlab("Tanggal") +
scale_x_date(date_breaks="1 week", date_labels="%d-%b", limits = as.Date(c("2023-05-20", "2023-06-24")))
p + geom_text_repel(
aes(Date, Epinow, label = Epinow), data = data_ends,
fontface ="bold", color = "black", size = 4) + theme_minimal()
p2 <- ggplot() + geom_line(data = R_teunis, aes(x = Date, y = Teunis), size = 1.5, color = "turquoise4") +
geom_line(data = R_teunis, aes(x = Date, y = treshold), size = 0.5, color = "black", linetype = "dashed") +
theme(plot.title=element_text(size=16), panel.grid.minor = element_blank()) + theme(strip.text = element_text(face="bold", size=9)) + labs(title ="Rt Wallinga & Teunis") + ylab("Rt") + xlab("Tanggal") +
scale_x_date(date_breaks="1 week", date_labels="%d-%b", limits = as.Date(c("2023-05-20", "2023-06-24")))
p2 + geom_text_repel(
aes(Date, Teunis, label = Teunis), data = data_ends2,
fontface ="bold", color = "black", size = 4) + theme_minimal()
p3 <- ggplot() + geom_line(data = R_epiestim, aes(x = Date, y = Epiestim), size = 1.5, color = "turquoise4") +
geom_line(data = R_epiestim, aes(x = Date, y = treshold), size = 0.5, color = "black", linetype = "dashed") +
theme(plot.title=element_text(size=16), panel.grid.minor = element_blank()) + theme(strip.text = element_text(face="bold", size=9)) + labs(title ="Rt Epiestim") + ylab("Rt") + xlab("Tanggal") +
scale_x_date(date_breaks="1 week", date_labels="%d-%b", limits = as.Date(c("2023-05-20", "2023-06-24")))
p3 + geom_text_repel(
aes(Date, Epiestim, label = Epiestim), data = data_ends3,
fontface ="bold", color = "black", size = 4) + theme_minimal()
p4 <- ggplot() + geom_line(data = R_bettencourt, aes(x = Date, y = bettencourt), size = 1.5, color = "turquoise4") +
geom_line(data = R_bettencourt, aes(x = Date, y = treshold), size = 0.5, color = "black", linetype = "dashed") +
theme(plot.title=element_text(size=16), panel.grid.minor = element_blank()) + theme(strip.text = element_text(face="bold", size=9)) + labs(title ="Rt Bettencourt & Ribiero") + ylab("Rt") + xlab("Tanggal") +
scale_x_date(date_breaks="1 week", date_labels="%d-%b", limits = as.Date(c("2023-05-20", "2023-06-24"))) + ylim(0.6,1.3)
p4 + geom_text_repel(
aes(Date, bettencourt, label = bettencourt), data = data_ends4,
fontface ="bold", color = "black", size = 4) + theme_minimal()
p5 <- ggplot() + geom_line(data = Rt_Alvarez, aes(x = Date, y = EpiInvert), size = 1.5, color = "turquoise4") +
geom_line(data = Rt_Alvarez, aes(x = Date, y = treshold), size = 0.5, color = "black", linetype = "dashed") +
theme(plot.title=element_text(size=16), panel.grid.minor = element_blank()) + theme(strip.text = element_text(face="bold", size=9)) + labs(title ="Rt Epilinvert") + ylab("Rt") + xlab("Tanggal") +
scale_x_date(date_breaks="1 week", date_labels="%d-%b", limits = as.Date(c("2023-05-20", "2023-06-24")))
p5 + geom_text_repel(
aes(Date, EpiInvert, label = EpiInvert), data = data_ends5,
fontface ="bold", color = "black", size = 4) + theme_minimal()
Referensi
- Cori A (2021). EpiEstim: Estimate Time Varying Reproduction Numbers from Epidemic Curves. R package version 2.2-4, https://CRAN.R-project.org/package=EpiEstim.
- Sam Abbott, Joel Hellewell, Katharine Sherratt, Katelyn Gostic, Joe Hickson, Hamada S. Badr, Michael DeWitt, Robin Thompson, EpiForecasts, Sebastian Funk (2020). EpiNow2: Estimate Real-Time Case Counts and Time-Varying Epidemiological Parameters. doi:10.5281/zenodo.3957489 https://doi.org/10.5281/zenodo.3957489.
- Boelle P, Obadia T (2022). R0: Estimation of R0 and Real-Time Reproduction Number from Epidemics. R package version 1.2-10, https://CRAN.R-project.org/package=R0.
- Alvarez L (2022). EpiInvert: Variational Techniques in Epidemiology. R package version 0.3.1, https://CRAN.R-project.org/package=EpiInvert.