inst/viral_load_simulate_ct.R

In RussellXu/rtestimate: Estimation Of Reproductive Number

#' using seir partab simulating ct data
library(rtestimate)
library(tidyverse)
library(patchwork)
library(foreach)
library(doParallel)
library(writexl)
library(ggplot2)
library(extraDistr)
library(odin) ## install from CRAN
library(virosolver)
library(dplyr)
library(xlsx)

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

load('./data/R0_partab.rda')


R0_partab$values[R0_partab$names=='R0'] <- 2
#R0_partab$values[R0_partab$names=='I0'] <- 0.002
#R0_partab$values[R0_partab$names=='t0'] <- 10
#R0_partab$values[R0_partab$names=='incubation'] <- 5

pars <- R0_partab$values
names(pars) <- R0_partab$names


## Calculate true Rt from true model parameters
# epidemic_process <- virosolver::simulate_seir_process(pars,solve_times,population_n)
# res <- epidemic_process$seir_outputs %>%
#   pivot_wider(names_from="variable",values_from="value") %>%
#   dplyr::select(time, Rt) %>%
#   write_xlsx('./output/True_rt_r2.xlsx')

# inci_data <- inci %>%
#   group_by(time) %>%
#   mutate(row = row_number())  %>%
#   pivot_wider(names_from = time, values_from = incidence) %>%
#   select(-row)
#inci_data <- data.frame(t(apply(inci_data, 1, quantile, probs=c(0.5, 0.75, 0.25), na.rm=T)))
#colnames(inci_data) <- c('incidence', 'inci_upper_band', 'inci_lower_band')
#write_xlsx(inci_data, '~/Desktop/prism_use/incidence_r2.xlsx')
###########################


#colnames(a) <- c('time', 'lower95', 'median', 'upper95')

# write.csv(as.integer(floor(seir_dynamics$incidence)), file = './data/hk_r10_incidence.csv')

## Simulate onset times, confirmation delays etc
## This returns a tibble with line list entries for **every** individual in the population



## Simulate SEIR dynamics, incidence, growth rates, Rt etc
## Use "ode" for deterministic
## Use "odin" for stochastic
seir_dynamics <- rtestimate::simulate_seir_wrapper_rte(population_n=population_n,solve_times=solve_times,
                                                       pars=pars, ver="odin")
seir_dynamics$plot
#save(seir_dynamics, file = 'hk_r2_seir_dynamics.rda')

#incidence <- readxl::read_excel('~/Desktop/prism_use/incidence_R2.xlsx')
complete_linelist <- simulate_observations_wrapper(seir_dynamics$incidence,times=solve_times,
                                                               population_n=population_n)

########################################
## Simulate observation process
########################################
## This bit of code generates linelist data for some specified observation process.
## Here, we simulate sampling 1000 people at random from the population every 14 days
## arguments changed above
sample_probs <- c(rep(0, sampling_frequency-1),sampling_number/population_n)
sample_probs <- rep(sample_probs, length(solve_times)/sampling_frequency +1)
sample_probs <- sample_probs[1:length(solve_times)]
frac_report <- tibble(t=solve_times,prob=sample_probs)
frac_report <- frac_report %>% filter(t >= 14 & t <= 100)


# have plot
observed_linelist <- simulate_reporting(complete_linelist,
                                        frac_report=NULL,
                                        timevarying_prob=frac_report,
                                        solve_times=solve_times,
                                        symptomatic=FALSE)

## Simulate viral load/Ct value for each person in the line list
simulated_viral_loads <- simulate_viral_loads_wrapper(observed_linelist$sampled_individuals,
                                                      kinetics_pars=pars)

## Clean data to form expected by virosolver
obs_dat <- simulated_viral_loads %>% dplyr::select(sampled_time, ct_obs) %>%
  rename(t = sampled_time, ct=ct_obs) %>% arrange(t)

## Save simulated line list, Cts and SEIR dynamics
#write_csv(seir_dynamics$seir_outputs, path=paste0(save_file,"_seir_outputs.csv"))
#write_csv(complete_linelist, path=paste0(save_file,"_full_linelist.csv"))
obs_dat <- drop_na(obs_dat)
#write_csv(obs_dat, "~/Desktop/prism_use/ct_r10.csv")

#save(obs_dat, file = ('ct_data.Rda'))
#ggsave(paste0(save_file,"_ct_dist.pdf"),p_dat,width=7,height=4,dpi=150)
#ggsave(paste0(save_file,"_seir_model.pdf"),seir_dynamics$plot,width=7,height=7,dpi=150)
#ggsave(paste0(save_file,"_growth_rates.pdf"),seir_dynamics$growth_rate_p,width=7,height=6,dpi=150)

#example_ct_data <- read.csv('~/Desktop/prism_use/rt2/ct_r2.csv')
ct_data <- read_csv('~/Desktop/prism_use/rt2/ct_r2.csv')
head(ct_data)