inst/viral_load_demo.R

In RussellXu/rtestimate: Estimation Of Reproductive Number

# devtools::install('./out_packages/lazymcmc-master/')
# devtools::install('./out_packages/virosolver-master/')

library(rtestimate)
library(tidyverse)
library(patchwork)
library(lazymcmc)
library(virosolver)
library(foreach)
library(doParallel)
library(writexl)
library(ggplot2)
library(extraDistr)
library(ggthemes)
library(odin) ## install from CRAN
library(fitdistrplus)

# generate ct data and return
#simulate_ct(population = 1000, times = 0:200, samp_freq = 7, samp_num = 100)
ct_data <- read.csv('~/Desktop/prism_use/rt2/ct_r2.csv')
load('data/R0_partab.rda')
#R0_partab$values[R0_partab$names=='t0'] <- 10
#R0_partab$values[R0_partab$names=='incubation'] <- 5
R0_partab$values[R0_partab$names=='R0'] <- 2
#R0_partab$values[R0_partab$names=='I0'] <- 0.002
# R0_widetype_partab <- example_seir_partab
#save(R0_partab, file = 'R0_partab.rda')

# using ct data to estimate rt
Rt_quants <- viral_load_estimate_rt(ct_data, R0_partab)
#write_xlsx(Rt_quants[1], '~/Desktop/prism_use/viral_load_r10.xlsx')


#pars <- R0_partab$values
#names(pars) <- R0_partab$names
## Solve the Ct model over a range of times since infection (referred to as "ages")
#test_ages <- seq(1,50,by=1)
#
# ## This gives the modal Ct value
# cts <- viral_load_func(pars, test_ages)
#
# p_ct_model <- ggplot(data.frame(ct=c(40,cts),t=c(0,test_ages))) +
#   geom_line(aes(x=t,y=ct)) +
#   scale_y_continuous(trans="reverse",
#                      limits=c(40,15)) +
#   theme_bw()+
#   theme(panel.grid=element_blank(),panel.border=element_blank(),axis.line=element_line(size=0.5,colour='black')) +
#   ylab("Ct value") +
#   xlab("Days since infection")
# p_ct_model
# ## Note that this model does not solve for t=0,
# p_ct_data <- ggplot(ct_data %>% filter(ct < 40)) +
#   geom_violin(aes(x=t,group=t,y=ct),scale="width",fill="grey70",draw_quantiles=c(0.025,0.5,0.975)) +
#   geom_jitter(aes(x=t,y=ct),size=0.1,width=2,height=0) +
#   scale_y_continuous(trans="reverse") +
#   theme_bw() +
#   theme(panel.grid=element_blank(),panel.border=element_blank(),axis.line=element_line(size=0.5,colour='black')) +
#   xlim(40, 150) +
#   ylab("Ct value") +
#   xlab("Time")
# p_ct_data