Load long-term healthcare outbreak data

Data loaded includes number of cases since outbreak on a given day for each location, the number of outbreaks, the size of the facility, as well as the labels for the outbreaks.

n_outbreaks <- length(BC_LTHC_outbreaks_100Imputs[[100]]$capacity)
outbreak_sizes <- BC_LTHC_outbreaks_100Imputs[[100]]$capacity
outbreak_cases_series <- BC_LTHC_outbreaks_100Imputs[[100]]$time_series
ob_codes <- BC_LTHC_outbreaks_100Imputs[[100]]$Location


outbreak_cases <- BC_LTHC_outbreaks_100Imputs[[100]]$case_matrix

tmax <- 64

# plot cases from matrix
as_tibble(outbreak_cases,rownames="time") %>%
  mutate(time = as.double(time)) %>%
  pivot_longer(-time,names_to = "location",values_to="cases") %>%
  ggplot(aes(x=time,y=cases,color=location)) +
  geom_line() +
  theme(legend.position = "none")
#> Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if `.name_repair` is omitted as of tibble 2.0.0.
#> Using compatibility `.name_repair`.

Fit model

Fit model including estimating intervention. This code chunk is not run in the vignette, but provides the syntax for how to produce posterior samples. The bc_fit data object is including in the package.

stan_mod <- rstan::stan_model(system.file("stan", "hierarchical_SEIR_incidence_model.stan", package = "cr0eso"))

bc_fit <- seir_model_fit(
                       stan_model = stan_mod,
                       tmax,n_outbreaks,outbreak_cases,outbreak_sizes,
                       intervention_switch = TRUE,
                       multilevel_intervention = FALSE,
                       iter = 2000)

Model diagnostic checking

pairs(bc_fit$model, pars = c("params"))

rstan::check_hmc_diagnostics(bc_fit$model)
#> 
#> Divergences:
#> 0 of 4000 iterations ended with a divergence.
#> 
#> Tree depth:
#> 0 of 4000 iterations saturated the maximum tree depth of 10.
#> 
#> Energy:
#> E-BFMI indicated no pathological behavior.
# print(mod)

traceplot(bc_fit$model, pars = c("params[1]", "params[2]"), inc_warmup = TRUE, nrow = 2)


rhats <- bayesplot::rhat(bc_fit$model)
rhats[rhats>1.03]
#> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> 
#>   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA 
#> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> 
#>   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA 
#> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> 
#>   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA 
#> <NA> 
#>   NA
bayesplot::mcmc_rhat(rhats)
#> Warning: Dropped 49 NAs from 'new_rhat(rhat)'.


bayesplot::mcmc_nuts_divergence(bayesplot::nuts_params(bc_fit$model), 
                                bayesplot::log_posterior(bc_fit$model))

Re-label facilities

Extract the posterior sample and re-label the facilities by their estimated mean \(R_0^k\),

# Extract the posterior samples to a structured list:
posts <- rstan::extract(bc_fit$model)

location_labels <- spread_draws(bc_fit$model,r0k[location]) %>% 
  mutate(location = as.character(location)) %>%
  group_by(location) %>% 
  summarise(r0 = mean(r0k)) %>% 
  arrange(r0) %>% 
  mutate(site = LETTERS[1:n()]) %>%
  dplyr::select(location,site)

# get object of incidence and zeta
extracted_posts <- hom_extract_posterior_draws(posts,
                                               location_labels = location_labels) 
#> Warning in rpois(dplyr::n(), incidence): NAs produced

Plot model output


result <- hom_plot_r0_by_location(extracted_posts=extracted_posts)
#> Warning: Ignoring unknown parameters: fun.y

# plot results
show(result$plot + labs(y=TeX("$R_{0,k}$")))
#> No summary function supplied, defaulting to `mean_se()`


result$table %>%
  kableExtra::kbl() %>%
  kableExtra::kable_styling(bootstrap_options=bootstrap_options)

location	r0
Total	2.29 (0.22 - 6.75)
R	9.17 (7.16 - 11.97)
Q	8.17 (5.29 - 11.38)
P	6.09 (4.7 - 8.04)
O	5.73 (4.37 - 7.76)
N	4.55 (3.05 - 6.98)
M	3.32 (2.43 - 4.58)
L	3.27 (2.42 - 4.51)
K	3.2 (2.41 - 4.35)
J	2.9 (2.13 - 3.98)
I	2.1 (1.29 - 3.23)
H	1.57 (0.7 - 2.69)
G	1.52 (0.91 - 2.33)
F	1.47 (0.7 - 2.5)
E	1.2 (0.46 - 2.21)
D	0.93 (0.3 - 1.89)
C	0.86 (0.35 - 1.54)
B	0.79 (0.22 - 1.79)
A	0.56 (0.16 - 1.17)

Plot model fit to incidence


result <- hom_plot_incidence_by_location(extracted_posts=extracted_posts,
                                         outbreak_cases = outbreak_cases, 
                                         end_time = tmax+5,
                                         location_labels = location_labels)
# plot results
result$plot
#> Warning: Removed 18 row(s) containing missing values (geom_path).

Plot counterfactual scenario

The code below extracts and plots the counterfactual scenario and also provides a summary table of cases by location in the baseline, scenario, where there was no intervention and the difference and proportional difference representing the cases averted. The final row provides a total summary of cases.

result <- hom_plot_counterfactual_by_location(bc_fit,
                                      outbreak_cases = outbreak_cases,
                                      location_labels = location_labels)

# plot results
show(result$plot)



# show table of results
result$table %>%
  kableExtra::kable() %>%
  kableExtra::kable_styling(bootstrap_options = c("striped","responsive"))

location	baseline	counterfactual	averted	proportion_averted
A	10 (5 - 17)	17 (7 - 52)	7 (-4 - 38)	0.42 (-0.45 - 0.8)
B	3 (0 - 6)	6 (1 - 31)	3 (-2 - 27)	0.62 (-1.5 - 1)
C	10 (4 - 17)	23 (8 - 79.05)	13 (-1 - 67)	0.59 (-0.14 - 0.87)
D	4 (1 - 9)	10 (2 - 48)	6 (-2 - 42)	0.67 (-0.67 - 0.94)
E	5 (1 - 10)	18 (3 - 64)	13 (-1 - 56.05)	0.75 (-0.17 - 0.94)
F	6 (2 - 12)	30 (6 - 74)	24 (2 - 65)	0.81 (0.29 - 0.94)
G	14 (7 - 23)	76 (24 - 160)	62 (12 - 143)	0.81 (0.5 - 0.92)
H	6 (2 - 12)	31 (6 - 70)	25 (2 - 62)	0.81 (0.25 - 0.94)
I	9 (3 - 17)	60 (22 - 90)	50 (16 - 79)	0.84 (0.67 - 0.93)
J	30 (20 - 44)	164 (129 - 195)	132 (98 - 164)	0.81 (0.73 - 0.88)
K	46 (32 - 63)	215 (182 - 248)	168 (134 - 203)	0.78 (0.7 - 0.85)
L	23 (14 - 35)	105 (82 - 130)	81 (57 - 106)	0.78 (0.66 - 0.87)
M	32 (21 - 46)	142 (116 - 169)	109 (82 - 137)	0.77 (0.67 - 0.85)
N	16 (8 - 27)	42 (24 - 61)	26 (4.95 - 46)	0.61 (0.16 - 0.82)
O	83 (63 - 104)	142 (117 - 168)	59 (26 - 91)	0.42 (0.21 - 0.57)
P	94 (72 - 116)	148 (123 - 174)	55 (20 - 88)	0.37 (0.15 - 0.53)
Q	93 (62 - 118)	114 (89 - 146)	18.5 (-12.05 - 77)	0.16 (-0.13 - 0.55)
R	80 (63 - 100)	89 (70 - 109)	8 (-14 - 31)	0.09 (-0.18 - 0.31)
total	569 (510 - 627)	1459 (1244.95 - 1735)	890 (658.95 - 1184.05)	0.61 (0.52 - 0.69)

Provide posterior distribution for the proportion of cases averted by location,

show(result$violin_plot)
#> Warning: Removed 151 rows containing non-finite values (stat_ydensity).
#> Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
#> collapsing to unique 'x' values

#> Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
#> collapsing to unique 'x' values

#> Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
#> collapsing to unique 'x' values

#> Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
#> collapsing to unique 'x' values

#> Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
#> collapsing to unique 'x' values

#> Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
#> collapsing to unique 'x' values

Plot critical times and uncertainty


locations <- BC_LTHC_outbreaks_100Imputs[[100]]$Location
num_cases <- BC_LTHC_outbreaks_100Imputs[[100]]$num_cases
capacities <- BC_LTHC_outbreaks_100Imputs[[100]]$capacity

posts <- rstan::extract(bc_fit$model)


hom_plot_critical_times_by_location(posts,
                                    locations,
                                    num_cases,
                                    capacities,
                                    location_labels = location_labels)

Hierarchichal intervention model

Fit model including estimating intervention by each location in a hierarchical design,

bc_fit_zeta <- seir_model_fit(
  stan_model = stan_mod,
  tmax,n_outbreaks,outbreak_cases,outbreak_sizes,
  intervention_switch = TRUE,
  multilevel_intervention = TRUE,
                       iter = 2000)

Plot model output - zeta model

# Extract the posterior samples to a structured list:
posts <- rstan::extract(bc_fit_zeta$model)

extracted_posts <- hom_extract_posterior_draws(posts,
                                  location_labels = location_labels) 
#> Warning in rpois(dplyr::n(), incidence): NAs produced

result <- hom_plot_r0_by_location(extracted_posts=extracted_posts, sort=FALSE)
#> Warning: Ignoring unknown parameters: fun.y

# plot results
show(result$plot + labs(y=TeX("$R_{0,k}$")))
#> No summary function supplied, defaulting to `mean_se()`


result$table %>%
  kableExtra::kbl() %>%
  kableExtra::kable_styling(bootstrap_options=bootstrap_options)

location	r0
Total	2.96 (0.33 - 7.79)
R	9.13 (7.06 - 12.35)
Q	9.95 (7.44 - 13.21)
P	6.91 (5.31 - 9.01)
O	6.7 (5.06 - 8.9)
N	6.65 (4.1 - 10.32)
M	4.73 (3.08 - 7.23)
L	2.37 (1.6 - 3.43)
K	3.82 (2.67 - 5.55)
J	2.71 (1.83 - 4.03)
I	2.87 (1.61 - 5.14)
H	5.39 (2.45 - 9.52)
G	3.62 (1.85 - 6.54)
F	3.68 (1.8 - 6.9)
E	3.73 (1.7 - 7.46)
D	3.28 (1.59 - 6.81)
C	2.62 (1.35 - 4.87)
B	2.97 (1.38 - 5.99)
A	2.6 (1.26 - 4.96)

Plot hierarchical intervention strength by location (zeta)


result <- hom_plot_zeta_by_location(extracted_posts=extracted_posts,sort=FALSE)
#> Warning: Ignoring unknown parameters: fun.y

# plot results
show(result$plot)
#> No summary function supplied, defaulting to `mean_se()`


result$table %>%
  kableExtra::kbl() %>%
  kableExtra::kable_styling(bootstrap_options=bootstrap_options)

location	zeta
Total	0.1 (0.03 - 0.32)
R	0.03 (0.01 - 0.09)
Q	0.11 (0.09 - 0.14)
P	0.06 (0.05 - 0.08)
O	0.06 (0.05 - 0.08)
N	0.09 (0.05 - 0.13)
M	0.08 (0.05 - 0.11)
L	0.03 (0.01 - 0.05)
K	0.06 (0.04 - 0.09)
J	0.05 (0.03 - 0.07)
I	0.07 (0.03 - 0.13)
H	0.19 (0.1 - 0.33)
G	0.13 (0.07 - 0.22)
F	0.15 (0.07 - 0.27)
E	0.18 (0.09 - 0.35)
D	0.21 (0.09 - 0.42)
C	0.16 (0.08 - 0.32)
B	0.22 (0.09 - 0.46)
A	0.27 (0.13 - 0.51)

Plot model fit to incidence - zeta model

result <- hom_plot_incidence_by_location(extracted_posts=extracted_posts,
                                         outbreak_cases = outbreak_cases,
                                         end_time = tmax+5, 
                                         location_labels = location_labels,
                                         sort=FALSE)
# plot results
result$plot
#> Warning: Removed 18 row(s) containing missing values (geom_path).

Plot counterfactual scenario - zeta model