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

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

In order to subset the data we define the dates for rounds 14 to 18 of the study

min_date_r1 <- as.Date("2020-05-01")
min_date_r8 <- as.Date("2020-12-30")
max_date_r13 <- as.Date("2021-07-12")

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 rounds 1 to 13 of the study. Note that the dataframes pos and tot also include daily weighted counts by region as well as nationally and so all analysis can easily be extended to look at each region of England.

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

The p-spline model can then be plotted with estimated 95%CI for prevalence (all the raw data used to produce the plots is also avialble in the list returned). In the paper this analysis was perfomed for we only present the model for rounds 8-13.

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

print(p_spline_plot[[1]])

We can then calculate and plot the rolling two week average Reproduction number using a generation time suitable for Alpha

p_spline_Rt_A <- plot_p_spline_R(X = pos[pos$X>=min_date_r1 & pos$X<= max_date_r13,]$X,
                                   p_spline_fit = p_spline_mod_react, 
                                   target_dist_between_knots = 5,
                                   spline_degree = 3,
                                   ylim = 2.,
                                   n=2.29, b=0.36,
                                   tau_max = 14)

print(p_spline_Rt_A[[1]])

Or as a sensitivity analysis using a generation time approporate for the Delta variant.

p_spline_Rt_D <- plot_p_spline_R(X = pos[pos$X>=min_date_r1 & pos$X<= max_date_r13,]$X,
                                   p_spline_fit = p_spline_mod_react, 
                                   target_dist_between_knots = 5,
                                   spline_degree = 3,
                                   ylim = 2.,
                                   n=2.20, b=0.48,
                                   tau_max = 14)

print(p_spline_Rt_D[[1]])

Fitting the Segmented model to the data

We can fit the Segmented model to all of the REACT data for rounds 8 to 13 of the study.

#Setting the key dtaes of the changes in restrictions
res0 <- as.Date("2021-01-06")
res1 <- as.Date("2021-03-08")
res2 <- as.Date("2021-03-29")
res3 <- as.Date("2021-04-12")
res4 <- as.Date("2021-05-17")

rest_dates <- c(res0,res1,res2,res3,res4)

breakpoint_mod_react <- stan_breakpoint_model_weighted(pos[pos$X>=min_date_r8 & pos$X<= max_date_r13,]$X,
                                   pos[pos$X>=min_date_r8 & pos$X<= max_date_r13,]$England,
                                   tot[tot$X>=min_date_r8 & tot$X<= max_date_r13,]$England,
                                   restriction_dates = rest_dates,
                                   iter = 2000,
                                   warmup = 500,
                                   cores = 1)

We can then plot the estimated R's and multiplicative changes in R as inferred from the Segmented model

breakpoint_plot <- plot_breakpoint_R(pos[pos$X>=min_date_r8 & pos$X<= max_date_r13,]$X,
                                     breakpoint_mod_react,
                                     rest_dates)

breakpoint_plot[[1]]

Fitting the Segmented model with vaccine + variant effects

We can fit the Segmented model including effects due to vaccination and the Delta variant to all of the REACT data for rounds 8 to 13 of the study.

Vax1 <- readRDS(system.file("extdata", "Vax1.rds", package = "reactidd"))
Vax2 <- readRDS(system.file("extdata", "Vax2.rds", package = "reactidd"))
Delta <- readRDS(system.file("extdata", "Delta.rds", package = "reactidd"))

breakpoint_mod_complex_react <- stan_breakpoint_model_complex_weighted(pos[pos$X>=min_date_r8 & pos$X<= max_date_r13,]$X,
                                   pos[pos$X>=min_date_r8 & pos$X<= max_date_r13,]$England,
                                   tot[tot$X>=min_date_r8 & tot$X<= max_date_r13,]$England,
                                   restriction_dates = rest_dates,
                                   Vax1,
                                   Vax2,
                                   Delta,
                                   iter = 2000,
                                   warmup = 500,
                                   cores = 1)

We can then plot the estimated R's inferred from the Segmented model including vaccination and Delta proportion as additional effects.

breakpoint_complex_plot <- plot_breakpoint_complex_R(pos[pos$X>=min_date_r8 & pos$X<= max_date_r13,]$X,
                                                     breakpoint_mod_complex_react,
                                                     rest_dates,
                                                     ylim=3.0)

breakpoint_complex_plot[[1]]