In RussellXu/rtestimate: Estimation Of Reproductive Number

knitr::opts_chunk$set(
  eval = TRUE,
  warning = FALSE,
  collapse = TRUE,
  comment = "#>"
)

This vignette implements three methods in estimating reproduction number(indicating the trend of the epidemic) curves from different input data.This package can adapted for any scenario by changing the parameters or replacing source data.And each method can called easily through instances.So you can observe the predicted epidemic Rt from different perspectives.\ Furthermore,this package, rtestimate ,after be run will generates a series of excel sheets in output, which contain Rt,95 credible intervals on each model predictions and time bar.You also can use Prism graphpad to visualize the inferred Rt curves. Compare this to the true Rt curve to see which methods can do the estimation precisely.And I also offer the interface for calculating R-square, RMSE, etc, to indicate how goodness of this-time-prediction. The input data,output data and parameters requirements of these methods are introduced as follow,respectively.\

# first load the package
library(rtestimate)
library(dplyr)

Branching Process Model

This method developed by Joel Hellewell,Sam Abbott,et al, using a list of parameters to simulation the number of potential secondary cases produced by each individual is drawn from a negative binomial distribution with a mean equal to the reproduction number.\

input

The key parameters of this model are shown below:

• simulation times
• the number of initial cases
• basic reproduction number
• the proportion of transmission that occurred before symptom onset
• maximum cases
• maximum times
• basic reproduction number during isolation
• the serial interval distribution parameters of the delay from symptom onset to isolation
• the proportion of asymptomatic

# call the branching method easily
res <- branching_simulation(n.sim = 3, num.initial.cases = 5, r0community = 2.5,
                            prop.ascertain = 0, cap_cases = 4500, cap_max_days = 90,
                            r0isolated = 0, disp.com = 0.16, prop.asym=0,
                            disp.iso = 1, delay_shape = 1.651524, delay_scale = 4.287786,
                            k = 0, quarantine = TRUE)

output

For this demo, we assume 10 times simulations and 5 initial samples,with R0 equal to 2, each potential new infection was assigned a time of infection drawn from the serial interval distribution.Secondary cases were only created if the person with the infection had not been isolated by the time of infection, if the exposures occur before this time.Each contact is traced with probability $\rho$, with probability $1 - \rho$ they are missed by contact tracing.As so on, until there is no infection.

res_Tableau <- subset(res, select = -c(cases_per_gen)) 
head(res_Tableau)

EpiEstim Module

This method developed by Anne Cori, Neil M. Ferguson, Christophe Fraser, Simon Cauchemez,it is described in the following paper: A New Framework and Software to Estimate Time-Varying Reproduction Numbers During Epidemics.\ The calculation of EpiEstim in estimating Rt is given by the ratio of the number of new infections generated at time step $t$, $I_t$, to the total infectiousness of infected individuals at time $t$, given by $\sum_{s=1}^{t} I_{t-s} w_{s}$⁠, the sum of infection incidence up to time step $t-1$, weighted by the infectivity function $w_s$. Rt is the average number of secondary cases that each infected individual would infect if the conditions remained as they were at time $t$./

input

It need two datasets to to the estimation:

• the daily incidence list of onset
• the mean and standard deviation of the serial interval

As we want running in a stochastic module, so a parametric serial interval is considered.

# load the origin incidence data
load('/Users/rus.xu/Desktop/rtestimate/data/hk_r2_seir_dynamics.rda')
inci <- seir_dynamics$incidence
t_start = seq(2, length(inci) - 6)

res_parametric_si <- estimate_R(inci,
                                method="parametric_si",
                                config = make_config(list(
                                  t_start = t_start,
                                  t_end = t_start + 6,
                                  mean_si = 2.6,
                                  std_si = 1.5))
)


output

para_output <- res_parametric_si$R %>%
  dplyr::select('Quantile.0.025(R)', 'Quantile.0.975(R)', 'Median(R)') %>%
  dplyr::select('lower_band' = 1, 'upper_band' = 2, 'Rt' = 3)
head(para_output)
#plot(res_parametric_si, legend = FALSE, what='R')


---

Viral Load Model

This method developed by James Hay, to infer epidemic Rt by using viral load data.The package name's virosolver,the following works are built based on this package.\ This is a brand new method that uses information inherent in cycle threshold (Ct) values from reverse transcription quantitative polymerase chain reaction(RT-qPCR) tests to estimate the epidemic. Because of variation across individuals, samples, and testing platforms,it is hard to obtain the Ct values easily, so Ct values are generated by using SEIR model, Ct values provide a probabilistic measure of time since infection.The distribution of Ct values across positive specimens at a single time point reflects the epidemic trajectory: A growing epidemic will necessarily have a high proportion of recently infected individuals with high viral loads, whereas a declining epidemic will have more individuals with older infections and thus lower viral loads. Because of these changing proportions, the epidemic trajectory or Rt should be inferable from the distribution of Ct values collected in the Ct generation function.\


# the packages we needed
library(tidyverse)
library(patchwork)
library(lazymcmc)
library(virosolver)
library(foreach)
library(doParallel)
library(writexl)
library(ggplot2)
library(dplyr)
library(extraDistr)
library(ggthemes)
library(odin) ## install from CRAN
library(fitdistrplus)

input

It need two dataset to to the estimation:

• ct_data
• initial parameter table

ct_data needs to be generated using instance viral_load_simulate_ct.R.\ According to the virus variants, set the correct parameter value in the R0_partab.\

ct_data <- read.csv('/Users/rus.xu/Desktop/rtestimate/data/ct_r2_data.csv')
load('/Users/rus.xu/Desktop/rtestimate/data/R0_partab.rda')
#  scenario set 
R0_partab$values[R0_partab$names=='incubation'] <- 5 
R0_partab$values[R0_partab$names=='R0'] <- 2  # IMPORTANT 
# pre parameters 
print(R0_partab)
# Ct values 
print(head(ct_data))


# estimation of Rt values function
# Rt_quants <- viral_load_estimate_rt(ct_data, R0_partab)

# No use of Ct values, just SEIR model to estimate Rt.As the True curve for comparison.
pars <- R0_partab$values
names(pars) <- R0_partab$names
epidemic_process <- virosolver::simulate_seir_process(pars,0:200,8000)
res <- epidemic_process$seir_outputs %>%
  pivot_wider(names_from="variable",values_from="value") %>%
  dplyr::select(time, Rt)
head(res)

output

Use MCMC chains doing the prediction can generate a number of posterior draws, which allows us to construct credible intervals on the model predictions and to visualize the inferred Rt curve. Compare this to the true (blue) Rt curve to see how accurate the estimates are.

readxl::read_excel('/Users/rus.xu/Desktop/rtestimate/output/viral_load_r2output.xlsx', sheet = 1)


The comparison of Rt by EpiEstim and Viral load:

RussellXu/rtestimate documentation built on Jan. 1, 2022, 7:18 p.m.