cobeylab/rt_estimation

Code to test Rt estimates against synthetic data

Code for analyses and figures shown in:

Practical considerations for detecting changes in the effective reproductive number, Rt.

Last updated 27-Aug-2020.

Code

This directory contains functions and wrappers used to perform analyses and generate figures:

• simulation.R - code to generate synthetic data using an SIR or SEIR-type model, deterministic or stochastic.
• funs_simulation-sweep.R - wrapper functions for epidemic simulation.
• infer_times_of_infection_observation.R - functions to infer times of observation from SEIR times of infection, and to infer times of infection from times of observation by (1) drawing samples from a known delay distribution, or (2) shifting back in time by the mean delay to observation.
• rtlive.R and rtlive.stan together provide code to reproduce an adaptation of the Bettencourt & Ribeiro method for Rt estimation popularized by rt.live.
• util.R - various utility functions, including wrappers to estimate Rt using the methods of Cori et al., Wallinga & Teunis, and using methods adapted from Bettencourt & Ribeiro by rt.live. The first two methods are implemented in the package EpiEstim. The final method uses the rstan implementation above.
• caseR.R - Functions to calculate the exact case reproductive number within the synthetic data (dashed black lines shown in Fig. 2 and Fig. B.2).
• Rc_math.Rmd - Notes on the math used to calculate the case reproductive number exactly.
• Richardson_Lucy.R - Performs deconvolution.

Analyses

This directory contains scripts and notebooks used to run analyses and generate figures:

Workflow:

• 01-simulate_data.R - Specify inputs, generate synthetic data and save to a directory called R0-xx/.
• 02-... - Various notebooks estimate Rt from synthetic data and generate plots.
• Run_all_scripts.R - Runs the entire workflow. Comments within indicate which notebooks generate which figures.

rt_estimation's Issues

Use user's own p_delay distribution in "rtlive.stan"

I'm wondering how to eventually add any user's own p_delay_distribution in the stan model.

Because I see from here, line 280
cumulative_p_delay <- rep(1, length(onset_frac))
that the cumulative is set to 1, no matter one. This has the result, as pointed out, that once a person shows symptoms it suddenly is tested and found positive. But in reality, this is not likely to occur, therefore I would like to add my own cumulative delay distribution (already computed), which takes into account the average onset time for symptoms + the delay of being tested.
I think there should be some sort of convolution being done somewhere, but I cannot figure out exactly where.