README.md

In RussellXu/rtestimate: Estimation Of Reproductive Number

rtestimate

Estimation of time varying reproduction number during epidemics through multiple strategies; And also it's a tool to consider the best and the most precise way in inferring Rt epidemic curves.

This README provides a package guide for estimating the time-varying reproduction number by three methods which have representative,and include comparison among these methods.It aims to help users to make the best decision.As for the vignette, you can click here to see the details about the methods mentioned above.

This package contains three methods in estimation Rt:

• Branching process model; Developed by Joel Hellewell et al.
• EpiEstim; Developed and maintained by Anne Cori et al.
• Viral load; Developed and maintained by James A. Hay et al.

Note

This packages use three methods above doing this analysis,the most vital thing to do a comparison is to make a equivalence among the different methods despite variant inputs format.I use the SEIR dynamic model to generate input-data for both EpiEstim and Viral load at same time, respectively.\ In order to make comparison convenient and efficient, I adapted several functions from there origin packages rather than using the whole packages.But may be you will notice that I load the virolsolver eventually,because it do not needed to rewrite MCMC again,just using virosolver's API is fine.Also for convenience.

External packages required

The package relies on other R packages:

• lazymacmc; Used for the MCMC procedure in viral load method.
• virosolver; An r package to infer epidemic dynamics from virological data.This package need several functions to estimate Rt curves.

Click on the package name to clone the repository, or just copy the path and then install them.Alternatively,I have download the packages in the out_packages directory, you can install it directly from there.

Quick start

Only the comparison of two method in estimation Rt is demonstrated here.The process about every methods and every details of function are shown here.\

Script

The demonstrations are in the inst directory.Each method will named after METHOD_demo.

Pocess

1. The first step to using the package is to set parameters and generate data the methods need.Open script ground_truth.R to setup the population number, simulation time, and the vital parameter,R0.It will generate ground truth incidence and Rt data with its error band data, respectively. ```r

when R0 == 2, population = 8000,data format

rt_data <- read_xlsx('~/Desktop/prism_use/rt2/true_rt_R2.xlsx') New names: * `` -> ...1 head(rt_data)

A tibble: 6 × 4

...1 true_rt upper_band lower_band 1 1 2 2 2 2 2 2 2 2 3 3 2 2 2 4 4 2 2 2 5 5 2 2 2 6 6 2 2 2

incidence_data <- read_xlsx('~/Desktop/prism_use/rt2/incidence_R2.xlsx') New names: * `` -> ...1 head(incidence_data)

A tibble: 6 × 4

...1 incidence upper_band lower_band 1 1 0 0 0 2 2 0 0 0 3 3 0 0 0 4 4 0 0 0 5 5 0 0 0 6 6 0 0 0 ```

Default of upper band and lower band is 95CI.

incidence Rt

1. Use epiestim_demo.R doing estimation.Load the incidence data from step1, make sure that the first day of input data is import case.Then choose the best window size,I set 14 days as default.Last thing is to set a proper serial interval distribution. ```r

   para_output <- res_parametric_si$R %>%

2. dplyr::select('Median(R)', 'Quantile.0.975(R)', 'Quantile.0.025(R)') %>%
3. dplyr::select('epiestim_Rt' = 1, 'upper_band' = 2, 'lower_band' = 3) %>%

4. transform(sequence(nrow(res_parametric_si$R)))

   write_xlsx(para_output, '~/Desktop/prism_use/epiestim_hk_r2.xlsx')

   head(para_output) epiestim_Rt upper_band lower_band 1 1.607021 2.002294 1.267562 2 1.482412 1.839319 1.174979 3 1.434197 1.765855 1.146933 4 1.465901 1.784860 1.187425 5 1.454328 1.757135 1.188494 6 1.509468 1.802787 1.249841 ```

5. Use viral_load_simulate_ct.R to generate CT values as input for viral load method.Notice that set the corresponding R0 with ground truth process. ```r

   head(ct_data)

A tibble: 6 × 2

  t    ct

1 55 40 2 55 40 3 55 40 4 55 40 5 55 40 6 55 40 ```

1. Use viral_load_demo.R doing estimation.Load the CT values data and R0 parameter table. ```r

   read_excel('~/Desktop/prism_use/rt2/viral_load_r2.xlsx')

A tibble: 201 × 4

time lower95 Rt_d upper95

1 0 1.94 2.33 3.41 2 1 1.94 2.33 3.41 3 2 1.94 2.33 3.41 4 3 1.94 2.33 3.41 5 4 1.94 2.33 3.41 6 5 1.94 2.33 3.41 7 6 1.94 2.33 3.41 8 7 1.94 2.33 3.41 9 8 1.93 2.33 3.41 10 9 1.93 2.33 3.41

… with 191 more rows

```

Which one is better?

• Calculation of RMSE & R squared & pearson, evaluating each estimation with ground truth.

> # the whole line
> # epiestim method: RMSE & R squared
> caret::postResample(pred = rt2_table$epiestim, obs = rt2_table$true_rt)
     RMSE  Rsquared     
0.3958170 0.9748548 
> # pearson
> cor(rt2_table$epiestim, rt2_table$true_rt, method = 'pearson')
[1] 0.9873473
> # viral load method: RMSE & R suqared 
> caret::postResample(pred = rt2_table$viral_load, obs = rt2_table$true_rt)
     RMSE  Rsquared     
0.2253348 0.9423817 
> # pearson
> cor(rt2_table$viral_load, rt2_table$true_rt, method = 'pearson')
[1] 0.9707635

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