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 REACT data available using reactidd::load_example_data()

pos <- load_example_data()[[1]]
tot <- load_example_data()[[2]]

In order to subset the data we define the dates of each of the first 7 rounds 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-20")
min_date_r5 <- as.Date("2020-09-18")
min_date_r6 <- as.Date("2020-10-16")
min_date_r7 <- as.Date("2020-11-13")

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")
max_date_r5 <- as.Date("2020-10-05")
max_date_r6 <- as.Date("2020-11-02")
max_date_r7 <- as.Date("2020-12-03")

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_react_r1 <- stan_exp_model(pos[pos$X>=min_date_r1 & pos$X<= max_date_r1,]$X,
                                   pos[pos$X>=min_date_r1 & pos$X<= max_date_r1,]$England,
                                   tot[tot$X>=min_date_r1 & tot$X<= max_date_r1,]$England,
                                   iter = 5000,
                                   warmup = 500,
                                   cores = 1)
exp_mod_react_r2 <- stan_exp_model(pos[pos$X>=min_date_r2 & pos$X<= max_date_r2,]$X,
                                   pos[pos$X>=min_date_r2 & pos$X<= max_date_r2,]$England,
                                   tot[tot$X>=min_date_r2 & tot$X<= max_date_r2,]$England,
                                   iter = 5000,
                                   warmup = 500,
                                   cores = 1)
exp_mod_react_r3 <- stan_exp_model(pos[pos$X>=min_date_r3 & pos$X<= max_date_r3,]$X,
                                   pos[pos$X>=min_date_r3 & pos$X<= max_date_r3,]$England,
                                   tot[tot$X>=min_date_r3 & tot$X<= max_date_r3,]$England,
                                   iter = 5000,
                                   warmup = 500,
                                   cores = 1)
exp_mod_react_r4 <- stan_exp_model(pos[pos$X>=min_date_r4 & pos$X<= max_date_r4,]$X,
                                   pos[pos$X>=min_date_r4 & pos$X<= max_date_r4,]$England,
                                   tot[tot$X>=min_date_r4 & tot$X<= max_date_r4,]$England,
                                   iter = 5000,
                                   warmup = 500,
                                   cores = 1)
exp_mod_react_r5 <- stan_exp_model(pos[pos$X>=min_date_r5 & pos$X<= max_date_r5,]$X,
                                   pos[pos$X>=min_date_r5 & pos$X<= max_date_r5,]$England,
                                   tot[tot$X>=min_date_r5 & tot$X<= max_date_r5,]$England,
                                   iter = 5000,
                                   warmup = 500,
                                   cores = 1)
exp_mod_react_r6 <- stan_exp_model(pos[pos$X>=min_date_r6 & pos$X<= max_date_r6,]$X,
                                   pos[pos$X>=min_date_r6 & pos$X<= max_date_r6,]$England,
                                   tot[tot$X>=min_date_r6 & tot$X<= max_date_r6,]$England,
                                   iter = 5000,
                                   warmup = 500,
                                   cores = 1)
exp_mod_react_r7 <- stan_exp_model(pos[pos$X>=min_date_r7 & pos$X<= max_date_r7,]$X,
                                   pos[pos$X>=min_date_r7 & pos$X<= max_date_r7,]$England,
                                   tot[tot$X>=min_date_r7 & tot$X<= max_date_r7,]$England,
                                   iter = 5000,
                                   warmup = 500,
                                   cores = 1)
exp_mod_react_r12 <- stan_exp_model(pos[pos$X>=min_date_r1 & pos$X<= max_date_r2,]$X,
                                   pos[pos$X>=min_date_r1 & pos$X<= max_date_r2,]$England,
                                   tot[tot$X>=min_date_r1 & tot$X<= max_date_r2,]$England,
                                   iter = 5000,
                                   warmup = 500,
                                   cores = 1)
exp_mod_react_r23 <- stan_exp_model(pos[pos$X>=min_date_r2 & pos$X<= max_date_r3,]$X,
                                   pos[pos$X>=min_date_r2 & pos$X<= max_date_r3,]$England,
                                   tot[tot$X>=min_date_r2 & tot$X<= max_date_r3,]$England,
                                   iter = 5000,
                                   warmup = 500,
                                   cores = 1)
exp_mod_react_r34 <- stan_exp_model(pos[pos$X>=min_date_r3 & pos$X<= max_date_r4,]$X,
                                   pos[pos$X>=min_date_r3 & pos$X<= max_date_r4,]$England,
                                   tot[tot$X>=min_date_r3 & tot$X<= max_date_r4,]$England,
                                   iter = 5000,
                                   warmup = 500,
                                   cores = 1)
exp_mod_react_r45 <- stan_exp_model(pos[pos$X>=min_date_r4 & pos$X<= max_date_r5,]$X,
                                   pos[pos$X>=min_date_r4 & pos$X<= max_date_r5,]$England,
                                   tot[tot$X>=min_date_r4 & tot$X<= max_date_r5,]$England,
                                   iter = 5000,
                                   warmup = 500,
                                   cores = 1)
exp_mod_react_r56 <- stan_exp_model(pos[pos$X>=min_date_r5 & pos$X<= max_date_r6,]$X,
                                   pos[pos$X>=min_date_r5 & pos$X<= max_date_r6,]$England,
                                   tot[tot$X>=min_date_r5 & tot$X<= max_date_r6,]$England,
                                   iter = 5000,
                                   warmup = 500,
                                   cores = 1)
exp_mod_react_r67 <- stan_exp_model(pos[pos$X>=min_date_r6 & pos$X<= max_date_r7,]$X,
                                   pos[pos$X>=min_date_r6 & pos$X<= max_date_r7,]$England,
                                   tot[tot$X>=min_date_r6 & tot$X<= max_date_r7,]$England,
                                   iter = 5000,
                                   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_react_r1 <- exponential_estimate_R(exp_mod_react_r1, n_mean = 2.29, b_mean =0.36, label ="React-Round1")
R_estimates_react_r2 <- exponential_estimate_R(exp_mod_react_r2, n_mean = 2.29, b_mean =0.36, label ="React-Round2")
R_estimates_react_r3 <- exponential_estimate_R(exp_mod_react_r3, n_mean = 2.29, b_mean =0.36, label ="React-Round3")
R_estimates_react_r4 <- exponential_estimate_R(exp_mod_react_r4, n_mean = 2.29, b_mean =0.36, label ="React-Round4")
R_estimates_react_r5 <- exponential_estimate_R(exp_mod_react_r5, n_mean = 2.29, b_mean =0.36, label ="React-Round5")
R_estimates_react_r6 <- exponential_estimate_R(exp_mod_react_r6, n_mean = 2.29, b_mean =0.36, label ="React-Round6")
R_estimates_react_r7 <- exponential_estimate_R(exp_mod_react_r7, n_mean = 2.29, b_mean =0.36, label ="React-Round7")
R_estimates_react_r12 <- exponential_estimate_R(exp_mod_react_r12, n_mean = 2.29, b_mean =0.36, label ="React-Round1&2")
R_estimates_react_r23 <- exponential_estimate_R(exp_mod_react_r23, n_mean = 2.29, b_mean =0.36, label ="React-Round2&3")
R_estimates_react_r34 <- exponential_estimate_R(exp_mod_react_r34, n_mean = 2.29, b_mean =0.36, label ="React-Round3&4")
R_estimates_react_r45 <- exponential_estimate_R(exp_mod_react_r45, n_mean = 2.29, b_mean =0.36, label ="React-Round3&4")
R_estimates_react_r56 <- exponential_estimate_R(exp_mod_react_r56, n_mean = 2.29, b_mean =0.36, label ="React-Round3&4")
R_estimates_react_r67 <- exponential_estimate_R(exp_mod_react_r67, n_mean = 2.29, b_mean =0.36, label ="React-Round3&4")

R_table <- rbind(R_estimates_react_r1, R_estimates_react_r2, R_estimates_react_r3, R_estimates_react_r4,R_estimates_react_r5,R_estimates_react_r6,R_estimates_react_r7,
                 R_estimates_react_r12, R_estimates_react_r23, R_estimates_react_r34, R_estimates_react_r45, R_estimates_react_r56, R_estimates_react_r67)
print(R_table)

Plotting the exponential model fits

The model fits can then be plotted. First the individual round fits:

individual_round_plots <- plot_exp_model(X = pos[pos$X>=min_date_r1 & pos$X<= max_date_r7,]$X,
                                   Y= pos[pos$X>=min_date_r1 & pos$X<= max_date_r7,]$England,
                                   N = tot[tot$X>=min_date_r1 & tot$X<= max_date_r7,]$England,
                                   fit_exp = list(exp_mod_react_r1, exp_mod_react_r2, exp_mod_react_r3, exp_mod_react_r4, exp_mod_react_r5, exp_mod_react_r6, exp_mod_react_r7),
                                   X_model = list(pos[pos$X>=min_date_r1 & pos$X<= max_date_r1,]$X,
                                                  pos[pos$X>=min_date_r2 & pos$X<= max_date_r2,]$X,
                                                  pos[pos$X>=min_date_r3 & pos$X<= max_date_r3,]$X,
                                                  pos[pos$X>=min_date_r4 & pos$X<= max_date_r4,]$X,
                                                  pos[pos$X>=min_date_r5 & pos$X<= max_date_r5,]$X,
                                                  pos[pos$X>=min_date_r6 & pos$X<= max_date_r6,]$X,
                                                  pos[pos$X>=min_date_r7 & pos$X<= max_date_r7,]$X),
                                   color_list = list("red","red","red","red","red","red","red"),
                                   ylim = 3.0)

print(individual_round_plots[[1]])

Then plot the models fit to subsequent rounds

subsequent_round_plots <- plot_exp_model(X = pos[pos$X>=min_date_r1 & pos$X<= max_date_r7,]$X,
                                   Y= pos[pos$X>=min_date_r1 & pos$X<= max_date_r7,]$England,
                                   N = tot[tot$X>=min_date_r1 & tot$X<= max_date_r7,]$England,
                                   fit_exp = list(exp_mod_react_r12, exp_mod_react_r23, exp_mod_react_r34, exp_mod_react_r45, exp_mod_react_r56, exp_mod_react_r67),
                                   X_model = list(pos[pos$X>=min_date_r1 & pos$X<= max_date_r2,]$X,
                                                  pos[pos$X>=min_date_r2 & pos$X<= max_date_r3,]$X,
                                                  pos[pos$X>=min_date_r3 & pos$X<= max_date_r4,]$X,
                                                  pos[pos$X>=min_date_r4 & pos$X<= max_date_r5,]$X,
                                                  pos[pos$X>=min_date_r5 & pos$X<= max_date_r6,]$X,
                                                  pos[pos$X>=min_date_r6 & pos$X<= max_date_r7,]$X),
                                   color_list = list("red","blue","dark green","cyan","purple","green"),
                                   ylim = 3.0)

print(subsequent_round_plots[[1]])

Fitting the Bayesian P-Spline Model to the data

We can then fit the Bayesian P-spline model to all of the REACT data for the first 7 rounds of the study (iterations and warmup will need to be bigger for accurate estimates)

p_spline_mod_react <- stan_p_spline(pos[pos$X>=min_date_r1 & pos$X<= max_date_r7,]$X,
                                   pos[pos$X>=min_date_r1 & pos$X<= max_date_r7,]$England,
                                   tot[tot$X>=min_date_r1 & tot$X<= max_date_r7,]$England,
                                   target_dist_between_knots = 5,
                                   spline_degree = 3,
                                   iter = 5000,
                                   warmup = 500,
                                   cores = 1)

The p-spline model can then be plotted with estimated 95%CI and 50%CI (the raw data used to produced the plots are avialble by using [[2]] and [[3]] instead of [[1]])

p_spline_plot <- plot_p_spline_prev(pos[pos$X>=min_date_r1 & pos$X<= max_date_r7,]$X,
                                   pos[pos$X>=min_date_r1 & pos$X<= max_date_r7,]$England,
                                   tot[tot$X>=min_date_r1 & tot$X<= max_date_r7,]$England,
                                   p_spline_fit = p_spline_mod_react, 
                                    target_dist_between_knots = 5,
                                    spline_degree = 3,
                                   ylim = 2.0)

print(p_spline_plot[[1]])

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

p_spline_min_date <- plot_p_spline_minimum_density(X = pos[pos$X>=min_date_r1 & pos$X<= max_date_r7,]$X,
                                   p_spline_fit = p_spline_mod_react, 
                                    target_dist_between_knots = 5,
                                    spline_degree = 3)

print(p_spline_min_date[[1]])

We can then calculate and plot the instanteous growth rate over the study period

p_spline_igr <- plot_p_spline_igr(X = pos[pos$X>=min_date_r1 & pos$X<= max_date_r7,]$X,
                                   p_spline_fit = p_spline_mod_react, 
                                    target_dist_between_knots = 5,
                                    spline_degree = 3,
                                  ylim = 0.2)

print(p_spline_igr[[1]])

and also the rolling two week average Reproduction number

p_spline_Rt <- plot_p_spline_R(X = pos[pos$X>=min_date_r1 & pos$X<= max_date_r7,]$X,
                                   p_spline_fit = p_spline_mod_react, 
                                    target_dist_between_knots = 5,
                                    spline_degree = 3)

print(p_spline_Rt[[1]])

From the Bayeisan P-spline posterior we can also calculate the average growth rate over periods of time covering individual and subsequent rounds.

R_estimate_pspline_round1 <- estimate_p_spline_average_R(X = pos[pos$X>=min_date_r1 & pos$X<= max_date_r7,]$X,
                                   p_spline_fit = p_spline_mod_react, 
                                    target_dist_between_knots = 5,
                                    spline_degree = 3,
                                   min_date_num = min_date_r1,
                                   max_date_num = max_date_r1,
                                   label = "Round1")
R_estimate_pspline_round2 <- estimate_p_spline_average_R(X = pos[pos$X>=min_date_r1 & pos$X<= max_date_r7,]$X,
                                   p_spline_fit = p_spline_mod_react, 
                                    target_dist_between_knots = 5,
                                    spline_degree = 3,
                                   min_date_num = min_date_r2,
                                   max_date_num = max_date_r2,
                                   label = "Round2")
R_estimate_pspline_round3 <- estimate_p_spline_average_R(X = pos[pos$X>=min_date_r1 & pos$X<= max_date_r7,]$X,
                                   p_spline_fit = p_spline_mod_react, 
                                    target_dist_between_knots = 5,
                                    spline_degree = 3,
                                   min_date_num = min_date_r3,
                                   max_date_num = max_date_r3,
                                   label = "Round3")
R_estimate_pspline_round4 <- estimate_p_spline_average_R(X = pos[pos$X>=min_date_r1 & pos$X<= max_date_r7,]$X,
                                   p_spline_fit = p_spline_mod_react, 
                                    target_dist_between_knots = 5,
                                    spline_degree = 3,
                                   min_date_num = min_date_r4,
                                   max_date_num = max_date_r4,
                                   label = "Round4")
R_estimate_pspline_round5 <- estimate_p_spline_average_R(X = pos[pos$X>=min_date_r1 & pos$X<= max_date_r7,]$X,
                                   p_spline_fit = p_spline_mod_react, 
                                    target_dist_between_knots = 5,
                                    spline_degree = 3,
                                   min_date_num = min_date_r5,
                                   max_date_num = max_date_r5,
                                   label = "Round5")
R_estimate_pspline_round6 <- estimate_p_spline_average_R(X = pos[pos$X>=min_date_r1 & pos$X<= max_date_r7,]$X,
                                   p_spline_fit = p_spline_mod_react, 
                                    target_dist_between_knots = 5,
                                    spline_degree = 3,
                                   min_date_num = min_date_r6,
                                   max_date_num = max_date_r6,
                                   label = "Round6")
R_estimate_pspline_round7 <- estimate_p_spline_average_R(X = pos[pos$X>=min_date_r1 & pos$X<= max_date_r7,]$X,
                                   p_spline_fit = p_spline_mod_react, 
                                    target_dist_between_knots = 5,
                                    spline_degree = 3,
                                   min_date_num = min_date_r7,
                                   max_date_num = max_date_r7,
                                   label = "Round7")

R_estimate_pspline_round12 <- estimate_p_spline_average_R(X = pos[pos$X>=min_date_r1 & pos$X<= max_date_r7,]$X,
                                   p_spline_fit = p_spline_mod_react, 
                                    target_dist_between_knots = 5,
                                    spline_degree = 3,
                                   min_date_num = min_date_r1,
                                   max_date_num = max_date_r2,
                                   label = "Round12")
R_estimate_pspline_round23 <- estimate_p_spline_average_R(X = pos[pos$X>=min_date_r1 & pos$X<= max_date_r7,]$X,
                                   p_spline_fit = p_spline_mod_react, 
                                    target_dist_between_knots = 5,
                                    spline_degree = 3,
                                   min_date_num = min_date_r2,
                                   max_date_num = max_date_r3,
                                   label = "Round23")
R_estimate_pspline_round34 <- estimate_p_spline_average_R(X = pos[pos$X>=min_date_r1 & pos$X<= max_date_r7,]$X,
                                   p_spline_fit = p_spline_mod_react, 
                                    target_dist_between_knots = 5,
                                    spline_degree = 3,
                                   min_date_num = min_date_r3,
                                   max_date_num = max_date_r4,
                                   label = "Round34")
R_estimate_pspline_round45 <- estimate_p_spline_average_R(X = pos[pos$X>=min_date_r1 & pos$X<= max_date_r7,]$X,
                                   p_spline_fit = p_spline_mod_react, 
                                    target_dist_between_knots = 5,
                                    spline_degree = 3,
                                   min_date_num = min_date_r4,
                                   max_date_num = max_date_r5,
                                   label = "Round45")
R_estimate_pspline_round56 <- estimate_p_spline_average_R(X = pos[pos$X>=min_date_r1 & pos$X<= max_date_r7,]$X,
                                   p_spline_fit = p_spline_mod_react, 
                                    target_dist_between_knots = 5,
                                    spline_degree = 3,
                                   min_date_num = min_date_r5,
                                   max_date_num = max_date_r6,
                                   label = "Round56")
R_estimate_pspline_round67 <- estimate_p_spline_average_R(X = pos[pos$X>=min_date_r1 & pos$X<= max_date_r7,]$X,
                                   p_spline_fit = p_spline_mod_react, 
                                    target_dist_between_knots = 5,
                                    spline_degree = 3,
                                   min_date_num = min_date_r6,
                                   max_date_num = max_date_r7,
                                   label = "Round67")


R_table <- rbind(R_estimate_pspline_round1,R_estimate_pspline_round2,R_estimate_pspline_round3,
                 R_estimate_pspline_round4,R_estimate_pspline_round5,R_estimate_pspline_round6,
                 R_estimate_pspline_round7,
                 R_estimate_pspline_round12,R_estimate_pspline_round23,R_estimate_pspline_round34,
                 R_estimate_pspline_round45,R_estimate_pspline_round56,R_estimate_pspline_round67)
print(R_table)