inst/ground_truth.R

In RussellXu/rtestimate: Estimation Of Reproductive Number

library(readxl)


population_n=8000
solve_times=0:200
sampling_frequency=14
sampling_number=1000

load('./data/R0_partab.rda')
#R0_partab$values[R0_partab$names=='incubation'] <- 5
R0_partab$values[R0_partab$names=='R0'] <- 2
pars <- R0_partab$values
names(pars) <- R0_partab$names

inci_rt <- data.frame()
for(n in 1:100){
  seir_dynamic <- simulate_seir_wrapper_rte(population_n=population_n,solve_times=solve_times,
                                                         pars=pars, ver="odin")

  # add incidence error bar data , 100 times simulation
  temp_inci_rt <- data.frame(incidence = seir_dynamic$incidence, true_rt = seir_dynamic$Rt$value, time = n)
  inci_rt <- rbind(inci_rt, temp_inci_rt)
}
inci_rt_data <- inci_rt %>%
  group_by(time) %>%
  mutate(row = row_number()) %>%
  pivot_wider(names_from = time, values_from = c(incidence, true_rt)) %>%
  dplyr::select(-row)

# calculate true rt and error bar
true_rt <- dplyr::select(inci_rt_data, num_range('true_rt_', 1:100)) %>%
  apply(1, quantile, probs=c(0.5, 0.975, 0.025), na.rm=T) %>%
  t() %>%
  data.frame() %>%
  dplyr::select(true_rt=1,true_rt_upper=2, true_rt_lower=3)
# calculate incidence and error bar
incidence <- dplyr::select(inci_rt_data, num_range('incidence_', 1:100)) %>%
  apply(1, quantile, probs=c(0.5, 0.975, 0.025), na.rm=T) %>%
  t() %>%
  data.frame() %>%
  dplyr::select(incidence=1,incidence_upper=2, incidence_lower=3)
#write.xlsx(true_rt, '~/Desktop/prism_use/true_rt_R10.xlsx')
#write.xlsx(incidence, '~/Desktop/prism_use/incidence_R10.xlsx')