In mrc-ide/reactidd: Package for the temporal analysis of the REACT study

knitr::opts_chunk$set(echo = TRUE)

library(reactidd)

Loading the Data

First we will load the example REACT data available using reactidd::load_example_data()

phe <- load_example_phe_data()

In order to subset the data we define the dates of each round of the study

min_date_r1 <- as.Date("2020-05-01")
min_date_r2 <- as.Date("2020-06-19")
min_date_r3 <- as.Date("2020-07-24")
min_date_r4 <- as.Date("2020-08-22")

max_date_r1 <- as.Date("2020-06-01")
max_date_r2 <- as.Date("2020-07-07")
max_date_r3 <- as.Date("2020-08-11")
max_date_r4 <- as.Date("2020-09-08")

Fitting the exponential model

We fit the exponential model to subsets of the data corresponding to individual rounds and pairs of subsequent rounds

exp_mod_phe_r1 <- stan_exp_model_phe(X = phe[phe$date>=min_date_r1 & phe$date<= max_date_r1,]$date,
                                   Y= phe[phe$date>=min_date_r1 & phe$date<= max_date_r1,]$n_cases,
                                   iter = 20000,
                                   warmup = 500,
                                   cores = 1)
exp_mod_phe_r2 <- stan_exp_model_phe(X = phe[phe$date>=min_date_r2 & phe$date<= max_date_r2,]$date,
                                   Y= phe[phe$date>=min_date_r2 & phe$date<= max_date_r2,]$n_cases,
                                   iter = 20000,
                                   warmup = 500,
                                   cores = 1)
exp_mod_phe_r3 <- stan_exp_model_phe(X = phe[phe$date>=min_date_r3 & phe$date<= max_date_r3,]$date,
                                   Y= phe[phe$date>=min_date_r3 & phe$date<= max_date_r3,]$n_cases,
                                   iter = 20000,
                                   warmup = 500,
                                   cores = 1)
exp_mod_phe_r4 <- stan_exp_model_phe(X = phe[phe$date>=min_date_r4 & phe$date<= max_date_r4,]$date,
                                   Y= phe[phe$date>=min_date_r4 & phe$date<= max_date_r4,]$n_cases,
                                   iter = 20000,
                                   warmup = 500,
                                   cores = 1)

exp_mod_phe_r12 <- stan_exp_model_phe(X = phe[phe$date>=min_date_r1 & phe$date<= max_date_r2,]$date,
                                   Y= phe[phe$date>=min_date_r1 & phe$date<= max_date_r2,]$n_cases,
                                   iter = 20000,
                                   warmup = 500,
                                   cores = 1)
exp_mod_phe_r23 <- stan_exp_model_phe(X = phe[phe$date>=min_date_r2 & phe$date<= max_date_r3,]$date,
                                   Y= phe[phe$date>=min_date_r2 & phe$date<= max_date_r3,]$n_cases,
                                   iter = 20000,
                                   warmup = 500,
                                   cores = 1)
exp_mod_phe_r34 <- stan_exp_model_phe(X = phe[phe$date>=min_date_r3 & phe$date<= max_date_r4,]$date,
                                   Y= phe[phe$date>=min_date_r3 & phe$date<= max_date_r4,]$n_cases,
                                   iter = 20000,
                                   warmup = 500,
                                   cores = 1)

Using these model fits we can calculate the growth rate, R and the doubling/halving times for each model fit

R_estimates_phe_r1 <- exponential_estimate_R(exp_mod_phe_r1, n_mean = 2.29, b_mean =0.36, label ="PHE-Round1")
R_estimates_phe_r2 <- exponential_estimate_R(exp_mod_phe_r2, n_mean = 2.29, b_mean =0.36, label ="PHE-Round2")
R_estimates_phe_r3 <- exponential_estimate_R(exp_mod_phe_r3, n_mean = 2.29, b_mean =0.36, label ="PHE-Round3")
R_estimates_phe_r4 <- exponential_estimate_R(exp_mod_phe_r4, n_mean = 2.29, b_mean =0.36, label ="PHE-Round4")
R_estimates_phe_r12 <- exponential_estimate_R(exp_mod_phe_r12, n_mean = 2.29, b_mean =0.36, label ="PHE-Round1&2")
R_estimates_phe_r23 <- exponential_estimate_R(exp_mod_phe_r23, n_mean = 2.29, b_mean =0.36, label ="PHE-Round2&3")
R_estimates_phe_r34 <- exponential_estimate_R(exp_mod_phe_r34, n_mean = 2.29, b_mean =0.36, label ="PHE-Round3&4")

R_table <- rbind(R_estimates_phe_r1, R_estimates_phe_r2, R_estimates_phe_r3, R_estimates_phe_r4,
                 R_estimates_phe_r12, R_estimates_phe_r23, R_estimates_phe_r34)
print(R_table)

We can then plot the exponential model with 95% CI's for individual rounds and subsequent rounds.

individual_round_plots <- plot_exp_model_phe(X = phe[phe$date>=min_date_r1 & phe$date<= max_date_r4,]$date,
                                   Y= phe[phe$date>=min_date_r1 & phe$date<= max_date_r4,]$n_cases,
                                   fit_exp = list(exp_mod_phe_r1, exp_mod_phe_r2, exp_mod_phe_r3, exp_mod_phe_r4),
                                   X_model = list(rev(phe[phe$date>=min_date_r1 & phe$date<= max_date_r1,]$date),
                                                  rev(phe[phe$date>=min_date_r2 & phe$date<= max_date_r2,]$date),
                                                  rev(phe[phe$date>=min_date_r3 & phe$date<= max_date_r3,]$date),
                                                  rev(phe[phe$date>=min_date_r4 & phe$date<= max_date_r4,]$date)),
                                   color_list = list("red","red","red","red"),
                                   ylim = 5000.0)


print(individual_round_plots[[1]])

subsequent_round_plots <-  plot_exp_model_phe(X = phe[phe$date>=min_date_r1 & phe$date<= max_date_r4,]$date,
                                   Y= phe[phe$date>=min_date_r1 & phe$date<= max_date_r4,]$n_cases,
                                   fit_exp = list(exp_mod_phe_r12, exp_mod_phe_r23, exp_mod_phe_r34),
                                   X_model = list(rev(phe[phe$date>=min_date_r1 & phe$date<= max_date_r2,]$date),
                                                  rev(phe[phe$date>=min_date_r2 & phe$date<= max_date_r3,]$date),
                                                  rev(phe[phe$date>=min_date_r3 & phe$date<= max_date_r4,]$date)),
                                   color_list = list("red","blue", "dark green"),
                                   ylim = 5000.0)

print(subsequent_round_plots[[1]])

Fitting the Bayesian P-spline model

We fit the bayesian p-spline model to all the data corresponding to the period of time of the first 4 rounds of the REACT study

p_spline_mod_phe <- stan_p_spline_phe(X = phe[phe$date>=min_date_r1 & phe$date<= max_date_r4,]$date,
                                   Y= phe[phe$date>=min_date_r1 & phe$date<= max_date_r4,]$n_cases,
                                   target_dist_between_knots = 5,
                                   spline_degree = 3,
                                   iter = 5000,
                                   warmup = 1000,
                                   cores = 1)

We can then plot the model fit with 95%CI and 50%CI

p_spline_plot <- plot_p_spline_phe(X = phe[phe$date>=min_date_r1 & phe$date<= max_date_r4,]$date,
                                   Y= phe[phe$date>=min_date_r1 & phe$date<= max_date_r4,]$n_cases,
                                   p_spline_fit = p_spline_mod_phe, 
                                    target_dist_between_knots = 5,
                                  spline_degree = 3,
                                   ylim = 5000.0)


print(p_spline_plot[[1]])

From the p-spline model we can estimate the date of minimum prevalence and plot the posterior distribution

p_spline_min_date <- plot_p_spline_minimum_density(X = phe[phe$date>=min_date_r1 & phe$date<= max_date_r4,]$date,
                                   p_spline_fit = p_spline_mod_phe, 
                                    target_dist_between_knots = 5,
                                    spline_degree = 3)

print(p_spline_min_date[[1]])