2. outcome

psborrow2 currently supports four outcomes:

• Time-to-event with exponential distribution (constant hazard), created with outcome_surv_exponential()
• Time-to-event with Weibull distribution and proportional hazards parametrization, created with outcome_surv_weibull_ph()
• Binary endpoints with a Bernoulli distribution and using logistic regression, created with outcome_bin_logistic()
• Continuous endpoints with a normal distribution, created with outcome_cont_normal()

After we select which outcome and distribution we want, we need to specify a prior distribution for the baseline event rate, baseline_prior. In this case, baseline_prior is a log hazard rate. Let’s assume we have no prior knowledge on this event rate, so we’ll specify an uninformative prior: prior_normal(0, 1000).

For our example, let’s conduct a time-to-event analysis using the exponential distribution.

outcome <- outcome_surv_exponential(
  time_var = "time",
  cens_var = "cnsr",
  baseline_prior = prior_normal(0, 1000)
)
outcome
# Outcome object with class OutcomeSurvExponential 
# 
# Outcome variables:
# time_var cens_var 
#   "time"   "cnsr" 
# 
# Baseline prior:
# Normal Distribution
# Parameters:
#  Stan  R    Value
#  mu    mean    0 
#  sigma sd   1000

3. borrowing

psborrow2 supports three different borrowing methods, each of which has its own class:

• No borrowing: This is the internal trial comparison without any external data. Use borrowing_none() to specify this.
• Full borrowing: This is pooling of the external and internal control arms. Use borrowing_full() to specify this.
• Bayesian dynamic borrowing with the hierarchical commensurate prior: This borrowing is as described in Hobbs et al. (2011) and uses the hierarchical commensurate prior. Use borrowing_hierarchical_commensurate() to specify this.

The column name for the external control column flag in our matrix is also required and passed to ext_flag_col.

Finally, for dynamic borrowing only, the hyperprior distribution on the commensurability parameter must be specified. This hyperprior determines (along with the comparability of the outcomes between internal and external controls) how much borrowing of the external control group will be performed. Example hyperpriors include largely uninformative inverse gamma distributions e.g., prior_gamma(alpha = .001, beta = .001) as well as more informative distributions e.g., prior_gamma(alpha = 1, beta = .001), though any distribution on the positive real line can be used. Distributions with more density at higher values (i.e., higher precision) will lead to more borrowing. We’ll choose an uninformative gamma prior in this example.

Note: Prior distributions are outlined in greater detail in a separate vignette, see vignette('prior_distributions', package = 'psborrow2').

borrowing <- borrowing_hierarchical_commensurate(
  ext_flag_col = "ext",
  tau_prior = prior_gamma(0.001, 0.001)
)
borrowing
# Borrowing object using the  Bayesian dynamic borrowing with the hierarchical commensurate prior  approach
# 
# External control flag: ext 
# 
# Commensurability parameter prior:
# Gamma Distribution
# Parameters:
#  Stan  R     Value
#  alpha shape 0.001
#  beta  rate  0.001
# Constraints: <lower=0>

3. treatment

Finally, treatment details are outlined in treatment_details(). Here, we first specify the column for the treatment flag in trt_flag_col. In addition, we need to specify the prior on the effect estimate, trt_prior. We’ll use another uninformative normal distribution for the prior on the treatment effect:

treatment <- treatment_details(
  trt_flag_col = "trt",
  trt_prior = prior_normal(0, 1000)
)
treatment
# Treatment object
# 
# Treatment flag column: trt 
# 
# Treatment effect prior:
# Normal Distribution
# Parameters:
#  Stan  R    Value
#  mu    mean    0 
#  sigma sd   1000

Application

Now that we have thought through each of the inputs to create_analysis_obj(), let’s create an analysis object:

anls_obj <- create_analysis_obj(
  data_matrix = example_matrix,
  outcome = outcome,
  borrowing = borrowing,
  treatment = treatment,
  quiet = TRUE
)

The Stan model compiled successfully and informed us that we are ready to begin sampling.

Note that if you are interested in seeing the Stan code that was generated, you can use the get_stan_code() function to see the full Stan code that will be compiled.

get_stan_code(anls_obj)
# 
#   functions {
#   
# }
# 
#   data {
#   int<lower=0> N;
#   vector[N] trt;
#   vector[N] time;
# vector[N] cens;
# 
#   matrix[N,2] Z;
#   
# }
# 
#   parameters {
# real beta_trt;
# 
# real<lower=0> tau;
# vector[2] alpha;
# 
# }
# 
#   transformed parameters {
#   real HR_trt = exp(beta_trt);
# }
# 
#   model {
#   vector[N] lp;
#   vector[N] elp;
#   beta_trt ~ normal(0, 1000);
#   lp = Z * alpha + trt * beta_trt;
# elp = exp(lp) ;
#   
#   
#   tau ~ gamma(0.001, 0.001) ;
# real sigma;
# sigma = 1 / tau;
# alpha[2] ~ normal(0, 1000) ;
# alpha[1] ~ normal(alpha[2], sqrt(sigma)) ;
#   for (i in 1:N) {
#    if (cens[i] == 1) {
#       target += exponential_lccdf(time[i] | elp[i] );
#    } else {
#       target += exponential_lpdf(time[i] | elp[i] );
#    }
# }
# }

Using bayesplot

With draws objects and the bayesplot package, we can create many useful visual summary plots. We can visualize the marginal posterior distribution of a variable we are interested in by plotting histograms of the draws with the function mcmc_hist(). Let’s do that for the Hazard ratio for the treatment effect and for our commensurability parameter, tau.

bayesplot::mcmc_hist(draws, c("treatment HR", "commensurability parameter"))

plot of chunk mcmc-hist

We can see other plots that help us understand and diagnose problems with the MCMC sampler, such as trace and rank plots:

bayesplot::color_scheme_set("mix-blue-pink")

bayesplot::mcmc_trace(
  draws[1:5000, 1:2, ], # Using a subset of draws only
  pars = c("treatment HR", "commensurability parameter"),
  n_warmup = 1000
)

plot of chunk mcmc-trace

Many other functions are outlined in the bayesplot vignettes.

Using posterior

draws objects are also supported by the posterior package, which provides many other tools for analyzing MCMC draw data. For instance, we can use the summarize_draws() function to derive 80% credible intervals for all parameters:

library(posterior)
summarize_draws(draws, ~ quantile(.x, probs = c(0.1, 0.9)))
# # A tibble: 6 × 3
#   variable                                `10%`     `90%`
#   <chr>                                   <dbl>     <dbl>
# 1 lp__                               -1620.     -1616.   
# 2 treatment log HR                      -0.408      0.100
# 3 commensurability parameter             0.0172     3.22 
# 4 baseline log hazard rate, internal    -3.57      -3.15 
# 5 baseline log hazard rate, external    -2.47      -2.33 
# 6 treatment HR                           0.665      1.11

Another useful application of the posterior package is the evaluation of the Monte Carlo standard error for quantiles of a variable of interest:

vm <- extract_variable_matrix(draws, "treatment HR")
mcse_quantile(x = vm, probs = c(0.1, 0.5, 0.9))
#  mcse_q10  mcse_q50  mcse_q90 
# 0.0006465 0.0006845 0.0012300

posterior contains many other helpful functions, as outlined in their vignettes.