Introduction"
In bennu: Bayesian Estimation of Naloxone Kit Number Under-Reporting

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(rstan)
library(bennu)
library(bayesplot)
library(ggplot2)

rstan_options(auto_write = TRUE)
options(mc.cores = 2)

## basic example code
d <- generate_model_data()
# note iter should be at least 2000 to generate a reasonable posterior sample
fit <- est_naloxone(d,iter=100)
mcmc_pairs(fit, pars = c("sigma","mu0","zeta"),
           off_diag_args = list(size = 1, alpha = 0.5))

Return median and 90th percentiles for posterior samples

library(posterior)

summarise_draws(fit, default_summary_measures())

We can compare two different model fits of the probability of prior kit use to simulated data using the plot_kit_use function

plot_kit_use(model = fit,data=d)

We can compare the posterior predictive check to the prior predictive check by re-running the model without the likelihood by setting the run_estimation parameter to FALSE

# note iter should be at least 2000 to generate a reasonable posterior sample
prior_fit <- est_naloxone(d,
                          run_estimation = FALSE,
                          iter=100)
mcmc_pairs(prior_fit, pars = c("sigma","mu0","zeta"),
           off_diag_args = list(size = 1, alpha = 0.5))

We can compare the prior and posterior predictive distribution of the probability of prior kit use using the plot_kit_use function

plot_kit_use(prior = prior_fit,posterior = fit, data=d)

Random walk of order 2 model

The second model uses a random walk of order 1 to characterize the time-dependence structure of the inverse logit probability of naloxone use,

$$\Delta^2c_i = c_i - c_{i+1} \sim N(0,\tau^{-1}) $$

The model can also be fit using a random walk of order 2 increments for the inverse logit probability of naloxone use using the following independent second order increments,

$$\Delta^2c_i = c_i - 2c_{i+1} + c_{i+2} \sim N(0,\tau^{-1}) $$

This version of the model will produce smoother estimates of probability of naloxone use in time.

# note iter should be at least 2000 to generate a reasonable posterior sample
# note use `rw_type` to specify order of random walk
rw2_fit <- est_naloxone(d,
                    rw_type = 2,
                    iter=100)
mcmc_pairs(rw2_fit, pars = c("sigma","mu0","zeta"),
           off_diag_args = list(size = 1, alpha = 0.5))

Summarize random walk order 2 model draws

summarise_draws(rw2_fit, default_convergence_measures())

We can compare two different model fits of the probability of prior kit use to simulated data using the plot_kit_use function

plot_kit_use(rw_1 = fit,rw_2 = rw2_fit,data=d)

fit missing distribution data

As distribution data are reported from sites some months may not be reported. This would lead to a difference in the frequency of the orders data and the distribution data. The model can account for these differences and infer the distribution during missing months. See the following example for generating distribution data once every three months

## basic example code
d_missing <- generate_model_data(reporting_freq = 3)
ggplot(aes(x=times,y=Reported_Used,color=as.factor(regions)),data=d_missing) +
  geom_point()

missing_fit <- est_naloxone(d_missing,
                    rw_type = 2,
                    iter=100)
mcmc_pairs(missing_fit, pars = c("sigma","mu0","zeta"),
           off_diag_args = list(size = 1, alpha = 0.5))

plot_kit_use(model = missing_fit, data=d_missing)

Change scale of data

the original model was designed for data aggregated into months. If the order and distribution data is aggregated differently then this can be incorporated into the model by updating the delay in reporting distribution and the delay in ordering to distributing distribution.

The delay in ordering to distributing distribution is composed of a gamma distribution with shape parameter $\alpha$ and inverse scale parameter $\beta = 1/\theta$. In order to update for example from monthly to weekly data, we can use the following property of the gamma distribution,

$$\Gamma(cx;\alpha,\beta) = \Gamma(x;\alpha,\beta/c)$$

To convert the delay distribution from month into week set $c = 1/4$ as one over the approximate number of weeks in a month. The reporting delay distribution is a simple vector of length 3 denoting the empirical delay distribution. This can be expanded into months through interpolation and normalization to retain its property as a probability distribution. Example code (not run) is shown below. Note that these distributions will need to be updated when considering another jurisdiction,

# monthly reporting delay distribution
psi_vec <- c(0.7, 0.2, 0.1)
# convert to weeks using interpolation
weekly_psi_vec <- rep(psi_vec,4) / sum(psi_vec)

# properties of order to distributed delay distribution in months
max_delays <- 3
delay_alpha <- 2
delay_beta <- 1

# convert to weeks
weekly_max_delays <- max_delays*4
weekly_delay_alpha <- delay_alpha
weekly_delay_beta <- 0.25 * delay_beta 

result <- est_naloxone(d,
 psi_vec = weekly_psi_vec,
 max_delays = weekly_max_delays,
 delay_alpha = weekly_delay_alpha,
 delay_beta = weekly_delay_beta)