2. Conduct a hybrid control analysis on a dataset using BDB
Matt Secrest and Isaac Gravestock
Source:vignettes/dataset.Rmd
dataset.Rmd
In this vignette, you’ll learn how to conduct Bayesian dynamic
borrowing (BDB) analyses using psborrow2
.
Before you start
The functionality in this article relies on Stan
for model fitting,
specifically via the CmdStan
and cmdstanr
tools.
If you haven’t used CmdStan
before you’ll need to
install the R package and the external program. More information can be
found in the cmdstanr
installation guide.
The short version is:
# Install the cmdstanr package
install.packages("cmdstanr", repos = c("https://standev.runiverse.dev", getOption("repos")))
library(cmdstanr)
# Install the external CmdStan program
check_cmdstan_toolchain()
install_cmdstan(cores = 2)
Now you’re ready to start with psborrow2
.
Creating an analysis object
For a BDB analysis in psborrow2
, we need to create an
object of class Analysis
which contains all the information
needed to build a model and compile an MCMC sampler using Stan. To
create an Analysis
object, we will call the function
create_analysis_obj()
. Let’s look at the four required
arguments to this function and evaluate them oneatatime.
create_analysis_obj(
data_matrix,
outcome,
borrowing,
treatment
)
1. data_matrix
Required elements
data_matrix
is where we input the onerowperpatient
numeric
matrix for our analysis. The column names of the
matrix are not fixed, so the names of columns will be specified in the
outcome, treatment, and borrowing sections.
There are two columns required for all analyses:
 A flag denoting receipt of the experimental intervention
(
1
) or not (0
)  A flag denoting whether the patient was part of the external data
source (
1
) or the internal trial (0
)
Timetoevent
If the outcome is timetoevent, then two additional columns are needed:
 The duration of followup for each patient
 A flag denoting whether the patient was censored (
1
) or not (0
)
Binary endpoints
If the outcome is binary, one additional column is needed:
 A flag denoting whether a patient had the event of interest
(
1
) or not (0
)
Covariates
Covariates may also be included in BDB analyses. These should be included in the data matrix if the plan is to adjust for them.
Note Only numeric
matrices are
supported. See Example data for creating
such a matrix from a data.frame
.
Note No missing data is currently allowed, all values must be nonmissing.
Example data
We will be using an example dataset stored in psborrow2
(example_matrix
). If you are starting from a data frame or
tibble, you can easily create a suitable matrix with the
psborrow2
helper function
create_data_matrix()
.
Creating a data matrix with create_data_matrix()
# Start with data.frame
diabetic_df < survival::diabetic
# For demonstration purposes, let some patients be external controls
diabetic_df$external < ifelse(diabetic_df$trt == 0 & diabetic_df$id > 1000, 1, 0)
# Create the censor flag
diabetic_df$cens < ifelse(diabetic_df$status == 0, 1, 0)
diabetes_matrix < create_data_matrix(
diabetic_df,
outcome = c("time", "cens"),
trt_flag_col = "trt",
ext_flag_col = "external",
covariates = ~ age + laser + risk
)
# Call `add_covariates()` with `covariates = c("age", "laserargon", "risk") `
head(diabetes_matrix)
# time cens trt external age laserargon risk
# 1 46.23 1 0 0 28 1 9
# 2 46.23 1 1 0 28 1 9
# 3 42.50 1 1 0 12 0 8
# 4 31.30 0 0 0 12 0 6
# 5 42.27 1 1 0 9 0 11
# 6 42.27 1 0 0 9 0 11
psborrow2
example matrix
Let’s look at the first few rows of the example matrix:
head(example_matrix)
# id ext trt cov4 cov3 cov2 cov1 time status cnsr resp
# [1,] 1 0 0 1 1 1 0 2.4226411 1 0 1
# [2,] 2 0 0 1 1 0 1 50.0000000 0 1 1
# [3,] 3 0 0 0 0 0 1 0.9674372 1 0 1
# [4,] 4 0 0 1 1 0 1 14.5774738 1 0 1
# [5,] 5 0 0 1 1 0 0 50.0000000 0 1 0
# [6,] 6 0 0 1 1 0 1 50.0000000 0 1 0
The column definitions are below:

ext
, 0/1, flag for external controls 
trt
, 0/1, flag for treatment arm 
cov1
, 0/1, a baseline covariate 
cov2
, 0/1, a baseline covariate 
time
, positive numeric, survival time 
cnsr
, 0/1, censoring indicator 
resp
, 0/1, indicator for binary response outcome
2. outcome
psborrow2
currently supports four outcomes:
 Timetoevent with exponential distribution (constant hazard),
created with
outcome_surv_exponential()
 Timetoevent 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 timetoevent 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] );
# }
# }
# }
Sampling from an analysis object
We can take draws from the posterior distribution using the function
mcmc_sample()
. This function takes as input our
Analysis
object and any arguments (other than the
data
argument) that are passed to CmdStanModel
objects. Note that running this may take a few minutes.
results < mcmc_sample(anls_obj,
iter_warmup = 2000,
iter_sampling = 50000,
chains = 4,
seed = 112233
)
Summarizing results
As a CmdStanMCMC
object, results
has
several methods which are outlined on the cmdstanr
website. For instance, we can see a see a summary of the posterior
distribution samples with results$summary()
:
results$summary()
# # A tibble: 6 × 10
# variable mean median sd mad q5 q95 rhat ess_bulk
# <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 lp__ 1618. 1618. 1.50 1.29 1621. 1.62e+3 1.00 71314.
# 2 beta_trt 0.155 0.158 0.199 0.198 0.477 1.74e1 1.00 85216.
# 3 tau 1.21 0.506 1.92 0.693 0.00418 4.74e+0 1.00 81024.
# 4 alpha[1] 3.36 3.35 0.161 0.161 3.63 3.10e+0 1.00 86029.
# 5 alpha[2] 2.40 2.40 0.0557 0.0557 2.49 2.31e+0 1.00 132062.
# 6 HR_trt 0.873 0.854 0.176 0.168 0.620 1.19e+0 1.00 85216.
# # ℹ 1 more variable: ess_tail <dbl>
The summary includes information for several parameter estimates from
our BDB model. Because it may not be immediately clear what the
parameters from the Stan model refer to, psborrow2
has a
function which returns a variable dictionary from the analysis object to
help interpret these parameters:
variable_dictionary(anls_obj)
# Stan_variable Description
# 1 tau commensurability parameter
# 2 alpha[1] baseline log hazard rate, internal
# 3 alpha[2] baseline log hazard rate, external
# 4 beta_trt treatment log HR
# 5 HR_trt treatment HR
We can also capture all of the draws by calling
results$draws()
, which returns an object of class
draws
. draws
objects are common in many MCMC
sampling software packages and allow us to leverage packages such as
posterior
and bayesplot
.
draws < results$draws()
print(draws)
# # A draws_array: 50000 iterations, 4 chains, and 6 variables
# , , variable = lp__
#
# chain
# iteration 1 2 3 4
# 1 1616 1616 1616 1616
# 2 1616 1616 1618 1618
# 3 1616 1617 1619 1617
# 4 1616 1618 1619 1618
# 5 1617 1617 1617 1617
#
# , , variable = beta_trt
#
# chain
# iteration 1 2 3 4
# 1 0.20 0.009 0.022 0.109
# 2 0.31 0.092 0.167 0.305
# 3 0.16 0.082 0.027 0.092
# 4 0.24 0.319 0.092 0.226
# 5 0.30 0.081 0.012 0.439
#
# , , variable = tau
#
# chain
# iteration 1 2 3 4
# 1 1.62 1.14 1.94 2.58
# 2 0.84 1.03 1.32 0.23
# 3 0.48 1.90 0.61 1.12
# 4 1.46 0.48 0.18 0.48
# 5 0.33 1.64 0.32 2.15
#
# , , variable = alpha[1]
#
# chain
# iteration 1 2 3 4
# 1 3.4 3.4 3.5 3.4
# 2 3.3 3.4 3.5 3.4
# 3 3.4 3.4 3.7 3.5
# 4 3.3 3.4 3.6 3.6
# 5 3.3 3.5 3.5 3.2
#
# # ... with 49995 more iterations, and 2 more variables
psborrow2
also has a function to rename variables in
draws
objects to be more interpretable,
rename_draws_covariates()
. This function uses the
variable_dictionary
labels. Let’s use it here to make the
results easier to interpret:
draws < rename_draws_covariates(draws, anls_obj)
summary(draws)
# # A tibble: 6 × 10
# variable mean median sd mad q5 q95 rhat ess_bulk
# <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 lp__ 1.62e+3 1.62e+3 1.50 1.29 1.62e+3 1.62e+3 1.00 71314.
# 2 treatment lo… 1.55e1 1.58e1 0.199 0.198 4.77e1 1.74e1 1.00 85216.
# 3 commensurabi… 1.21e+0 5.06e1 1.92 0.693 4.18e3 4.74e+0 1.00 81024.
# 4 baseline log… 3.36e+0 3.35e+0 0.161 0.161 3.63e+0 3.10e+0 1.00 86029.
# 5 baseline log… 2.40e+0 2.40e+0 0.0557 0.0557 2.49e+0 2.31e+0 1.00 132062.
# 6 treatment HR 8.73e1 8.54e1 0.176 0.168 6.20e1 1.19e+0 1.00 85216.
# # ℹ 1 more variable: ess_tail <dbl>
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.
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("mixbluepink")
bayesplot::mcmc_trace(
draws[1:5000, 1:2, ], # Using a subset of draws only
pars = c("treatment HR", "commensurability parameter"),
n_warmup = 1000
)
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.