inst/epiestim_demo.R

In RussellXu/rtestimate: Estimation Of Reproductive Number

library(rtestimate)
library(writexl)
library(ggplot2)
library(magrittr)
library(incidence)
library(gridExtra)
library(readxl)
library(dplyr)
#hk_inc <- read.csv('./data/hongkong_rt.csv')
#load('./data/hk_r2_seir_dynamics.rda')
seir_dynamics <- read_excel('~/Desktop/prism_use/rt2/incidence_R2.xlsx')
#inci <- read_excel('./onset_time_data.xlsx', sheet=4)
#hk_inc <- MockRotavirus
#plot(as.incidence(SARS2003$incidence, dates = Flu2009$incidence$dates))
t_start <- seq(2, length(seir_dynamics$incidence[12:125]) - 13)
t_end <- t_start + 13

res_parametric_si <- estimate_R(seir_dynamics$incidence[12:125],
                                method="parametric_si",
                                config = make_config(list(
                                  t_start = t_start,
                                  t_end = t_end,
                                  mean_si = 4.4,
                                  std_si = 3))
)
epiestim_rt_plot(res_parametric_si, legend = FALSE, what='R')
# save the data
para_output <- res_parametric_si$R %>%
  dplyr::select('Median(R)', 'Quantile.0.975(R)', 'Quantile.0.025(R)') %>%
  dplyr::select('epiestim_Rt' = 1, 'upper_band' = 2, 'lower_band' = 3) %>%
  transform(sequence(nrow(res_parametric_si$R)))
write_xlsx(para_output, '~/Desktop/prism_use/modified_rt2/origin_10.xlsx')
head(para_output)








estimate_R_plots <- function(..., legend = FALSE) {
  plot.estimate_R(..., legend = legend)
}








#  find the best std & mean values for estimation
ground_true <- read_excel('~/Desktop/prism_use/rt10/true_rt_R10.xlsx')
seir_dynamics <- read_excel('~/Desktop/prism_use/rt10/incidence_R10.xlsx')
t_start <- seq(2, length(seir_dynamics$incidence[12:60]) - 13)
t_end <- t_start + 13
find_precise_mean_si <- function(mean_si, seir_dynamics, t_start, t_end) {
  res_parametric_si <- estimate_R(seir_dynamics$incidence[12:60],
                                  method="parametric_si",
                                  config = make_config(list(
                                    t_start = t_start,
                                    t_end = t_end,
                                    mean_si = mean_si,
                                    std_si = 3.0))
  )
  predict <- res_parametric_si$R$`Mean(R)`
  return(predict)
}
mean_frame <- data.frame(index = seq(1,35))
for (i in seq(1.1,10,0.1)) {
  mean_temp <- find_precise_mean_si(i, seir_dynamics, t_start, t_end)
  mean_frame[as.character(i)] <- mean_temp
}

res <- mean_frame %>%
  select(-index) %>%
  plyr::adply(2, function(x) caret::postResample(pred = x, obs = ground_true$true_rt[15:60])[1])
res[which.min(res$RMSE),]
#write_xlsx(mean_frame, '~/Desktop/prism_use/modified_rt2/mean_frame_2.xlsx')