Usage

This usage tutorial shows how to fit historical control-group borrowing models using the historicalborrowlong package.

Data

historicalborrowlong expects data on multiple patients partitioned into studies, groups, and repeated measures (“reps”). Here is a simulated example. There are functions to simulate from the prior predictive distribution of each of the hierarchical, independent, and pooled models.

library(historicalborrowlong)
library(dplyr)
set.seed(0)
data <- hbl_sim_independent(
  n_continuous = 1,
  n_study = 3,
  n_group = 2,
  n_rep = 4,
  alpha = rep(1, 12),
  delta = rep(0.5, 4),
  n_patient = 10
)$data %>%
  rename(
    outcome = response,
    trial = study,
    arm = group,
    subject = patient,
    visit = rep,
    factor1 = covariate_study1_continuous1,
    factor2 = covariate_study2_continuous1
  ) %>%
  mutate(
    trial = paste0("trial", trial),
    arm = paste0("arm", arm),
    subject = paste0("subject", subject),
    visit = paste0("visit", visit),
  )
data
#> # A tibble: 160 × 8
#>    trial  arm   subject  visit  factor1 factor2 covariate_study3_conti…¹ outcome
#>    <chr>  <chr> <chr>    <chr>    <dbl>   <dbl>                    <dbl>   <dbl>
#>  1 trial1 arm1  subject1 visit1   0.764       0                        0  2.82  
#>  2 trial1 arm1  subject1 visit2   0.764       0                        0  3.40  
#>  3 trial1 arm1  subject1 visit3   0.764       0                        0  3.30  
#>  4 trial1 arm1  subject1 visit4   0.764       0                        0  3.02  
#>  5 trial1 arm1  subject2 visit1  -0.799       0                        0 -0.504 
#>  6 trial1 arm1  subject2 visit2  -0.799       0                        0  0.0952
#>  7 trial1 arm1  subject2 visit3  -0.799       0                        0  0.441 
#>  8 trial1 arm1  subject2 visit4  -0.799       0                        0 -0.0921
#>  9 trial1 arm1  subject3 visit1  -1.15        0                        0 -0.0738
#> 10 trial1 arm1  subject3 visit2  -1.15        0                        0 -1.13  
#> # ℹ 150 more rows
#> # ℹ abbreviated name: ¹covariate_study3_continuous1

You as the user will choose a reference level of the study (“trial”) column to indicate which study is the current one (the other are historical). Likewise, you will choose a level of the group (“arm”) column to indicate which group is the control group and a level of the rep (“visit”) column to indicate the first measurement of each patient (baseline). To see how historicalborrowlong assigns numeric indexes to the study and group levels, use hbl_data(). Viewing this output may assist with interpreting the results later on.

standardized_data <- hbl_data(
  data = data,
  response = "outcome",
  study = "trial",
  study_reference = "trial3",
  group = "arm",
  group_reference = "arm1",
  patient = "subject",
  rep = "visit",
  rep_reference = "visit1",
   # Can be continuous, categorical, or binary columns:
  covariates = c("factor1", "factor2")
)
standardized_data
#> # A tibble: 160 × 11
#>    study patient patient_label   rep rep_label response study_label group_label
#>    <int>   <int> <chr>         <int> <chr>        <dbl> <chr>       <chr>      
#>  1     1       1 subject1          1 visit1      2.82   trial1      arm1       
#>  2     1       1 subject1          2 visit2      3.40   trial1      arm1       
#>  3     1       1 subject1          3 visit3      3.30   trial1      arm1       
#>  4     1       1 subject1          4 visit4      3.02   trial1      arm1       
#>  5     1       2 subject10         1 visit1      0.343  trial1      arm1       
#>  6     1       2 subject10         2 visit2     -1.46   trial1      arm1       
#>  7     1       2 subject10         3 visit3     -0.747  trial1      arm1       
#>  8     1       2 subject10         4 visit4      0.155  trial1      arm1       
#>  9     1      12 subject2          1 visit1     -0.504  trial1      arm1       
#> 10     1      12 subject2          2 visit2      0.0952 trial1      arm1       
#> # ℹ 150 more rows
#> # ℹ 3 more variables: group <int>, covariate_factor1 <dbl>,
#> #   covariate_factor2 <dbl>

distinct(
  standardized_data,
  study,
  study_label,
  group,
  group_label,
  rep,
  rep_label
) %>%
  select(
    study,
    study_label,
    group,
    group_label,
    rep,
    rep_label
  )
#> # A tibble: 16 × 6
#>    study study_label group group_label   rep rep_label
#>    <int> <chr>       <int> <chr>       <int> <chr>    
#>  1     1 trial1          1 arm1            1 visit1   
#>  2     1 trial1          1 arm1            2 visit2   
#>  3     1 trial1          1 arm1            3 visit3   
#>  4     1 trial1          1 arm1            4 visit4   
#>  5     2 trial2          1 arm1            1 visit1   
#>  6     2 trial2          1 arm1            2 visit2   
#>  7     2 trial2          1 arm1            3 visit3   
#>  8     2 trial2          1 arm1            4 visit4   
#>  9     3 trial3          1 arm1            1 visit1   
#> 10     3 trial3          1 arm1            2 visit2   
#> 11     3 trial3          1 arm1            3 visit3   
#> 12     3 trial3          1 arm1            4 visit4   
#> 13     3 trial3          2 arm2            1 visit1   
#> 14     3 trial3          2 arm2            2 visit2   
#> 15     3 trial3          2 arm2            3 visit3   
#> 16     3 trial3          2 arm2            4 visit4

As explained in the hbl_data() and hb_mcmc_*() help files, before running the MCMC, dataset is pre-processed. This includes expanding the rows of the data so every rep of every patient gets an explicit row. So if your original data has irregular rep IDs, e.g. unscheduled visits in a clinical trial that few patients attend, please remove them before the analysis. Only the most common rep IDs should be added.

After expanding the rows, the function fills in missing values for every column except the response. That includes covariates. Missing covariate values are filled in, first with last observation carried forward, then with last observation carried backward. If there are still missing values after this process, the program throws an informative error.

Models

Recommended: run all models using `hbl_mcmc_sge()`

The hbl_mcmc_sge() function runs all 3 models of interest on a Sun Grid Engine (SGE) computing cluster and returns a list of results from all three models. For standardized analyses of real studies, this is highly recommended over the alternative of running each model separately (hbl_mcmc_hierarchical(), hbl_mcmc_pool(), and hbl_mcmc_independent(), as in the next subsection).

In hbl_mcmc_sge(), each model runs in its own SGE job, and chains run in parallel across the cores available to each job. The return value is a list of results from each model (hierarchical, independent, and pooled) that you would get by running hbl_mcmc_hierarchical(), hbl_mcmc_pool(), and hbl_mcmc_independent() separately. Each of these data frames each be supplied separately to hbl_summary() as explained later in this vignette.

mcmc <- hbl_mcmc_sge(
  data = data,
  response = "outcome",
  study = "trial",
  study_reference = "trial3",
  group = "arm",
  group_reference = "arm1",
  patient = "subject",
  rep = "visit",
  rep_reference = "visit1",
   # Can be continuous, categorical, or binary columns:
  covariates = c("factor1", "factor2"),
  # Raise these arguments for serious analyses:
  chains = 1, # Increase to about 3 or 4 in real-life use cases.
  cores = 1, # *HIGHLY* recommended to have cores = chains
  iter = 20, # Increase to several thousand in real-life use cases.
  warmup = 10, # Increase to several thousand in real-life use cases.
  log = "/dev/null", # optional SGE log file, /dev/null to disregard the log
  scheduler = "local" # Set to "sge" for serious analysis.
)

#> 
#> SAMPLING FOR MODEL 'historicalborrowlong' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.000103 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.03 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: WARNING: No variance estimation is
#> Chain 1:          performed for num_warmup < 20
#> Chain 1: 
#> Chain 1: Iteration:  1 / 20 [  5%]  (Warmup)
#> Chain 1: Iteration:  2 / 20 [ 10%]  (Warmup)
#> Chain 1: Iteration:  4 / 20 [ 20%]  (Warmup)
#> Chain 1: Iteration:  6 / 20 [ 30%]  (Warmup)
#> Chain 1: Iteration:  8 / 20 [ 40%]  (Warmup)
#> Chain 1: Iteration: 10 / 20 [ 50%]  (Warmup)
#> Chain 1: Iteration: 11 / 20 [ 55%]  (Sampling)
#> Chain 1: Iteration: 12 / 20 [ 60%]  (Sampling)
#> Chain 1: Iteration: 14 / 20 [ 70%]  (Sampling)
#> Chain 1: Iteration: 16 / 20 [ 80%]  (Sampling)
#> Chain 1: Iteration: 18 / 20 [ 90%]  (Sampling)
#> Chain 1: Iteration: 20 / 20 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 0.052 seconds (Warm-up)
#> Chain 1:                0.462 seconds (Sampling)
#> Chain 1:                0.514 seconds (Total)
#> Chain 1: 
#> 
#> SAMPLING FOR MODEL 'historicalborrowlong' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 7e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.7 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: WARNING: No variance estimation is
#> Chain 1:          performed for num_warmup < 20
#> Chain 1: 
#> Chain 1: Iteration:  1 / 20 [  5%]  (Warmup)
#> Chain 1: Iteration:  2 / 20 [ 10%]  (Warmup)
#> Chain 1: Iteration:  4 / 20 [ 20%]  (Warmup)
#> Chain 1: Iteration:  6 / 20 [ 30%]  (Warmup)
#> Chain 1: Iteration:  8 / 20 [ 40%]  (Warmup)
#> Chain 1: Iteration: 10 / 20 [ 50%]  (Warmup)
#> Chain 1: Iteration: 11 / 20 [ 55%]  (Sampling)
#> Chain 1: Iteration: 12 / 20 [ 60%]  (Sampling)
#> Chain 1: Iteration: 14 / 20 [ 70%]  (Sampling)
#> Chain 1: Iteration: 16 / 20 [ 80%]  (Sampling)
#> Chain 1: Iteration: 18 / 20 [ 90%]  (Sampling)
#> Chain 1: Iteration: 20 / 20 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 0.069 seconds (Warm-up)
#> Chain 1:                0.08 seconds (Sampling)
#> Chain 1:                0.149 seconds (Total)
#> Chain 1: 
#> 
#> SAMPLING FOR MODEL 'historicalborrowlong' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 6.6e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.66 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: WARNING: No variance estimation is
#> Chain 1:          performed for num_warmup < 20
#> Chain 1: 
#> Chain 1: Iteration:  1 / 20 [  5%]  (Warmup)
#> Chain 1: Iteration:  2 / 20 [ 10%]  (Warmup)
#> Chain 1: Iteration:  4 / 20 [ 20%]  (Warmup)
#> Chain 1: Iteration:  6 / 20 [ 30%]  (Warmup)
#> Chain 1: Iteration:  8 / 20 [ 40%]  (Warmup)
#> Chain 1: Iteration: 10 / 20 [ 50%]  (Warmup)
#> Chain 1: Iteration: 11 / 20 [ 55%]  (Sampling)
#> Chain 1: Iteration: 12 / 20 [ 60%]  (Sampling)
#> Chain 1: Iteration: 14 / 20 [ 70%]  (Sampling)
#> Chain 1: Iteration: 16 / 20 [ 80%]  (Sampling)
#> Chain 1: Iteration: 18 / 20 [ 90%]  (Sampling)
#> Chain 1: Iteration: 20 / 20 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 0.054 seconds (Warm-up)
#> Chain 1:                0.059 seconds (Sampling)
#> Chain 1:                0.113 seconds (Total)
#> Chain 1:

mcmc
#> $hierarchical
#> # A tibble: 10 × 69
#>    `alpha[1]` `alpha[2]` `alpha[3]` `alpha[4]` `alpha[5]` `alpha[6]` `alpha[7]`
#>         <dbl>      <dbl>      <dbl>      <dbl>      <dbl>      <dbl>      <dbl>
#>  1      2.37       1.32       0.679      1.42       3.06      -0.732      0.187
#>  2      0.778      1.75       1.28       0.870      0.748      0.603      1.30 
#>  3      0.770      1.14       1.02       0.866      0.722      0.655      1.02 
#>  4      1.04       1.58       1.18       1.23       0.945      1.08       1.18 
#>  5      0.566      2.07       1.10       0.506      0.497      0.911      1.10 
#>  6      0.521      1.61       1.00       0.896      0.431      0.970      1.01 
#>  7      0.810      1.14       1.08       0.724      0.868      0.986      1.09 
#>  8      1.03       0.963      1.03       0.752      0.619      0.966      0.529
#>  9      1.27       0.857      0.930      0.752      0.657      1.01       1.22 
#> 10      0.560      1.75       0.975      0.256      0.604      0.980      0.634
#> # ℹ 62 more variables: `alpha[8]` <dbl>, `alpha[9]` <dbl>, `alpha[10]` <dbl>,
#> #   `alpha[11]` <dbl>, `alpha[12]` <dbl>, `delta[1]` <dbl>, `delta[2]` <dbl>,
#> #   `delta[3]` <dbl>, `delta[4]` <dbl>, `beta[1]` <dbl>, `beta[2]` <dbl>,
#> #   `sigma[1,1]` <dbl>, `sigma[2,1]` <dbl>, `sigma[3,1]` <dbl>,
#> #   `sigma[1,2]` <dbl>, `sigma[2,2]` <dbl>, `sigma[3,2]` <dbl>,
#> #   `sigma[1,3]` <dbl>, `sigma[2,3]` <dbl>, `sigma[3,3]` <dbl>,
#> #   `sigma[1,4]` <dbl>, `sigma[2,4]` <dbl>, `sigma[3,4]` <dbl>, …
#> 
#> $independent
#> # A tibble: 10 × 61
#>    `alpha[1]` `alpha[2]` `alpha[3]` `alpha[4]` `alpha[5]` `alpha[6]` `alpha[7]`
#>         <dbl>      <dbl>      <dbl>      <dbl>      <dbl>      <dbl>      <dbl>
#>  1      1.00       1.45       1.01       0.902      0.408     -1.74        1.10
#>  2      1.37       0.852      0.948      0.839      0.923     -0.784       1.70
#>  3      0.913      0.914      1.07       1.20       0.284      0.935       1.26
#>  4      0.868      0.969      1.02       1.48       0.661      0.986       2.32
#>  5      1.05       1.01       1.03       1.52       0.264      0.992       2.24
#>  6      0.501      1.04       1.05       0.808      1.05       0.994       2.10
#>  7      0.965      1.16       1.00       0.559     -0.256      0.998       1.31
#>  8      0.830      1.10       1.02       0.612     -0.457      1.01        1.35
#>  9      1.49       0.927      0.891      0.522      1.24       0.981       1.07
#> 10      1.18       0.840      0.972      0.894      0.982      0.979       1.25
#> # ℹ 54 more variables: `alpha[8]` <dbl>, `alpha[9]` <dbl>, `alpha[10]` <dbl>,
#> #   `alpha[11]` <dbl>, `alpha[12]` <dbl>, `delta[1]` <dbl>, `delta[2]` <dbl>,
#> #   `delta[3]` <dbl>, `delta[4]` <dbl>, `beta[1]` <dbl>, `beta[2]` <dbl>,
#> #   `sigma[1,1]` <dbl>, `sigma[2,1]` <dbl>, `sigma[3,1]` <dbl>,
#> #   `sigma[1,2]` <dbl>, `sigma[2,2]` <dbl>, `sigma[3,2]` <dbl>,
#> #   `sigma[1,3]` <dbl>, `sigma[2,3]` <dbl>, `sigma[3,3]` <dbl>,
#> #   `sigma[1,4]` <dbl>, `sigma[2,4]` <dbl>, `sigma[3,4]` <dbl>, lp__ <dbl>, …
#> 
#> $pool
#> # A tibble: 10 × 53
#>    `alpha[1]` `alpha[2]` `alpha[3]` `alpha[4]` `delta[1]` `delta[2]` `delta[3]`
#>         <dbl>      <dbl>      <dbl>      <dbl>      <dbl>      <dbl>      <dbl>
#>  1      1.02       0.942      1.01       1.06     -0.797     -0.276   -0.665   
#>  2      0.941      1.02       1.04       0.927     0.528      0.0151   0.263   
#>  3      1.07       0.839      0.998      0.854    -0.311      0.310    0.187   
#>  4      0.741      1.17       1.04       1.00     -0.0223    -0.102   -0.272   
#>  5      0.994      1.04       1.04       0.964    -0.0103     0.275   -0.165   
#>  6      0.954      1.06       0.984      0.951    -0.493      0.765    0.278   
#>  7      1.02       0.958      1.02       0.874     0.207      0.681    0.453   
#>  8      0.974      0.962      1.02       0.908    -0.192      0.322    0.0158  
#>  9      0.867      0.955      0.873      1.03     -0.0356     0.0230  -0.000104
#> 10      0.936      0.970      0.969      1.01     -0.731      0.170   -0.384   
#> # ℹ 46 more variables: `delta[4]` <dbl>, `beta[1]` <dbl>, `beta[2]` <dbl>,
#> #   `sigma[1,1]` <dbl>, `sigma[2,1]` <dbl>, `sigma[3,1]` <dbl>,
#> #   `sigma[1,2]` <dbl>, `sigma[2,2]` <dbl>, `sigma[3,2]` <dbl>,
#> #   `sigma[1,3]` <dbl>, `sigma[2,3]` <dbl>, `sigma[3,3]` <dbl>,
#> #   `sigma[1,4]` <dbl>, `sigma[2,4]` <dbl>, `sigma[3,4]` <dbl>, lp__ <dbl>,
#> #   .chain <int>, .iteration <int>, .draw <int>, `lambda_current[1,2,1]` <dbl>,
#> #   `lambda_current[1,3,1]` <dbl>, `lambda_current[1,4,1]` <dbl>, …

Alternative: run specific models

The hierarchical model is the primary model of interest.¹

mcmc_hierarchical <- hbl_mcmc_hierarchical(
  data = data,
  response = "outcome",
  study = "trial",
  study_reference = "trial3",
  group = "arm",
  group_reference = "arm1",
  patient = "subject",
  rep = "visit",
  rep_reference = "visit1",
  # Can be continuous, categorical, or binary columns.
  covariates = c("factor1", "factor2"),
  # Raise these arguments for serious analyses:
  chains = 1, # Increase to about 3 or 4 in real-life use cases.
  iter = 400, # Increase to several thousand in real-life use cases.
  warmup = 200, # Increase to several thousand in real-life use cases.
  cores = 1 # Optionally run different chains in different processes.
)

#> 
#> SAMPLING FOR MODEL 'historicalborrowlong' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 7.2e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.72 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:   1 / 400 [  0%]  (Warmup)
#> Chain 1: Iteration:  40 / 400 [ 10%]  (Warmup)
#> Chain 1: Iteration:  80 / 400 [ 20%]  (Warmup)
#> Chain 1: Iteration: 120 / 400 [ 30%]  (Warmup)
#> Chain 1: Iteration: 160 / 400 [ 40%]  (Warmup)
#> Chain 1: Iteration: 200 / 400 [ 50%]  (Warmup)
#> Chain 1: Iteration: 201 / 400 [ 50%]  (Sampling)
#> Chain 1: Iteration: 240 / 400 [ 60%]  (Sampling)
#> Chain 1: Iteration: 280 / 400 [ 70%]  (Sampling)
#> Chain 1: Iteration: 320 / 400 [ 80%]  (Sampling)
#> Chain 1: Iteration: 360 / 400 [ 90%]  (Sampling)
#> Chain 1: Iteration: 400 / 400 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 8.99 seconds (Warm-up)
#> Chain 1:                9.457 seconds (Sampling)
#> Chain 1:                18.447 seconds (Total)
#> Chain 1:

mcmc_hierarchical
#> # A tibble: 200 × 69
#>    `alpha[1]` `alpha[2]` `alpha[3]` `alpha[4]` `alpha[5]` `alpha[6]` `alpha[7]`
#>         <dbl>      <dbl>      <dbl>      <dbl>      <dbl>      <dbl>      <dbl>
#>  1      0.867      1.00       1.05       0.813      0.571      0.984      0.927
#>  2      1.15       1.01       0.982      0.783      0.481      0.990      0.982
#>  3      0.824      1.13       1.01       0.535      0.609      0.975      1.03 
#>  4      1.30       0.923      0.984      0.979      0.252      0.978      1.01 
#>  5      0.956      0.902      0.953      0.878      1.13       0.957      0.898
#>  6      1.36       0.892      0.924      0.856      0.542      0.990      1.73 
#>  7      0.676      1.12       1.01       0.498      0.301      0.973      0.897
#>  8      1.18       0.875      0.929      1.02       0.531      0.996      1.07 
#>  9      1.20       0.844      1.00       0.393      0.569      0.994      1.20 
#> 10      0.942      1.01       1.06       0.954      0.462      0.984      0.977
#> # ℹ 190 more rows
#> # ℹ 62 more variables: `alpha[8]` <dbl>, `alpha[9]` <dbl>, `alpha[10]` <dbl>,
#> #   `alpha[11]` <dbl>, `alpha[12]` <dbl>, `delta[1]` <dbl>, `delta[2]` <dbl>,
#> #   `delta[3]` <dbl>, `delta[4]` <dbl>, `beta[1]` <dbl>, `beta[2]` <dbl>,
#> #   `sigma[1,1]` <dbl>, `sigma[2,1]` <dbl>, `sigma[3,1]` <dbl>,
#> #   `sigma[1,2]` <dbl>, `sigma[2,2]` <dbl>, `sigma[3,2]` <dbl>,
#> #   `sigma[1,3]` <dbl>, `sigma[2,3]` <dbl>, `sigma[3,3]` <dbl>, …

The pooled and independent models are benchmarks used to quantify the borrowing strength of the hierarchical model. To run these benchmark models, run the functions below. Each function returns a data frame with one column per parameter and one row per posterior sample.

mcmc_pool <- hbl_mcmc_pool(
  data = data,
  response = "outcome",
  study = "trial",
  study_reference = "trial3",
  group = "arm",
  group_reference = "arm1",
  patient = "subject",
  rep = "visit",
  rep_reference = "visit1",
  # Can be continuous, categorical, or binary columns:
  covariates = c("factor1", "factor2"),
  # Raise these arguments for serious analyses:
  chains = 1, # Increase to about 3 or 4 in real-life use cases.
  iter = 400, # Increase to several thousand in real-life use cases.
  warmup = 200, # Increase to several thousand in real-life use cases.
  cores = 1 # Optionally run different chains in different processes.
)

#> 
#> SAMPLING FOR MODEL 'historicalborrowlong' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 6.2e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.62 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:   1 / 400 [  0%]  (Warmup)
#> Chain 1: Iteration:  40 / 400 [ 10%]  (Warmup)
#> Chain 1: Iteration:  80 / 400 [ 20%]  (Warmup)
#> Chain 1: Iteration: 120 / 400 [ 30%]  (Warmup)
#> Chain 1: Iteration: 160 / 400 [ 40%]  (Warmup)
#> Chain 1: Iteration: 200 / 400 [ 50%]  (Warmup)
#> Chain 1: Iteration: 201 / 400 [ 50%]  (Sampling)
#> Chain 1: Iteration: 240 / 400 [ 60%]  (Sampling)
#> Chain 1: Iteration: 280 / 400 [ 70%]  (Sampling)
#> Chain 1: Iteration: 320 / 400 [ 80%]  (Sampling)
#> Chain 1: Iteration: 360 / 400 [ 90%]  (Sampling)
#> Chain 1: Iteration: 400 / 400 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 2.449 seconds (Warm-up)
#> Chain 1:                0.829 seconds (Sampling)
#> Chain 1:                3.278 seconds (Total)
#> Chain 1:

mcmc_pool
#> # A tibble: 200 × 53
#>    `alpha[1]` `alpha[2]` `alpha[3]` `alpha[4]` `delta[1]` `delta[2]` `delta[3]`
#>         <dbl>      <dbl>      <dbl>      <dbl>      <dbl>      <dbl>      <dbl>
#>  1      1.10       1.00       1.02       0.961     0.0574    -0.0813     0.127 
#>  2      0.971      1.01       0.982      1.02     -0.205      0.629      0.172 
#>  3      0.947      0.982      1.06       0.992     0.422      0.125      0.243 
#>  4      0.892      0.976      1.01       0.887    -0.293      0.413      0.0658
#>  5      0.794      0.974      0.972      0.860    -0.250      0.115      0.0172
#>  6      1.13       0.981      0.993      0.985     0.900      1.38       0.952 
#>  7      0.853      0.993      0.984      0.877     1.11       0.989      1.10  
#>  8      1.13       0.970      0.968      1.13      0.0338    -0.187      0.0172
#>  9      1.03       0.976      1.03       1.01     -0.0903     0.327      0.140 
#> 10      0.987      0.984      0.998      0.996     0.237      0.581      0.282 
#> # ℹ 190 more rows
#> # ℹ 46 more variables: `delta[4]` <dbl>, `beta[1]` <dbl>, `beta[2]` <dbl>,
#> #   `sigma[1,1]` <dbl>, `sigma[2,1]` <dbl>, `sigma[3,1]` <dbl>,
#> #   `sigma[1,2]` <dbl>, `sigma[2,2]` <dbl>, `sigma[3,2]` <dbl>,
#> #   `sigma[1,3]` <dbl>, `sigma[2,3]` <dbl>, `sigma[3,3]` <dbl>,
#> #   `sigma[1,4]` <dbl>, `sigma[2,4]` <dbl>, `sigma[3,4]` <dbl>, lp__ <dbl>,
#> #   .chain <int>, .iteration <int>, .draw <int>, …

mcmc_independent <- hbl_mcmc_independent(
  data = data,
  response = "outcome",
  study = "trial",
  study_reference = "trial3",
  group = "arm",
  group_reference = "arm1",
  patient = "subject",
  rep = "visit",
  rep_reference = "visit1",
  covariates = c("factor1", "factor2"),
  # Raise these arguments for serious analyses:
  chains = 1, # Increase to about 3 or 4 in real-life use cases.
  iter = 400, # Increase to several thousand in real-life use cases.
  warmup = 200, # Increase to several thousand in real-life use cases.
  cores = 1 # Optionally run different chains in different processes.
)

#> 
#> SAMPLING FOR MODEL 'historicalborrowlong' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 6.6e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.66 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:   1 / 400 [  0%]  (Warmup)
#> Chain 1: Iteration:  40 / 400 [ 10%]  (Warmup)
#> Chain 1: Iteration:  80 / 400 [ 20%]  (Warmup)
#> Chain 1: Iteration: 120 / 400 [ 30%]  (Warmup)
#> Chain 1: Iteration: 160 / 400 [ 40%]  (Warmup)
#> Chain 1: Iteration: 200 / 400 [ 50%]  (Warmup)
#> Chain 1: Iteration: 201 / 400 [ 50%]  (Sampling)
#> Chain 1: Iteration: 240 / 400 [ 60%]  (Sampling)
#> Chain 1: Iteration: 280 / 400 [ 70%]  (Sampling)
#> Chain 1: Iteration: 320 / 400 [ 80%]  (Sampling)
#> Chain 1: Iteration: 360 / 400 [ 90%]  (Sampling)
#> Chain 1: Iteration: 400 / 400 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 1.964 seconds (Warm-up)
#> Chain 1:                0.897 seconds (Sampling)
#> Chain 1:                2.861 seconds (Total)
#> Chain 1:

mcmc_independent
#> # A tibble: 200 × 61
#>    `alpha[1]` `alpha[2]` `alpha[3]` `alpha[4]` `alpha[5]` `alpha[6]` `alpha[7]`
#>         <dbl>      <dbl>      <dbl>      <dbl>      <dbl>      <dbl>      <dbl>
#>  1      1.71       1.06       0.931      0.693     0.112       0.973      0.903
#>  2      1.17       0.770      0.942      0.593     0.804       0.979      1.23 
#>  3      1.55       1.05       0.916      0.963     0.673       0.985      1.35 
#>  4      1.27       0.524      0.998      0.599     0.344       0.993      1.35 
#>  5      1.43       1.01       1.02       0.691     0.569       0.991      1.38 
#>  6      1.25       1.03       0.910      1.03      0.629       0.983      1.16 
#>  7      1.18       0.996      0.991      0.267     0.346       0.983      1.15 
#>  8      1.02       0.929      1.06       0.772     0.481       0.985      1.37 
#>  9      1.07       1.10       0.979      0.645     0.252       0.992      1.44 
#> 10      0.837      1.16       0.918      0.745    -0.0183      0.987      1.66 
#> # ℹ 190 more rows
#> # ℹ 54 more variables: `alpha[8]` <dbl>, `alpha[9]` <dbl>, `alpha[10]` <dbl>,
#> #   `alpha[11]` <dbl>, `alpha[12]` <dbl>, `delta[1]` <dbl>, `delta[2]` <dbl>,
#> #   `delta[3]` <dbl>, `delta[4]` <dbl>, `beta[1]` <dbl>, `beta[2]` <dbl>,
#> #   `sigma[1,1]` <dbl>, `sigma[2,1]` <dbl>, `sigma[3,1]` <dbl>,
#> #   `sigma[1,2]` <dbl>, `sigma[2,2]` <dbl>, `sigma[3,2]` <dbl>,
#> #   `sigma[1,3]` <dbl>, `sigma[2,3]` <dbl>, `sigma[3,3]` <dbl>, …

Model performance

A typical workflow will run all three models. Since each model can take several hours to run, it is strongly recommended to:

For each model, run all chains in parallel by setting the cores argument equal to chains.
Run all models concurrently on different jobs on a computing cluster. Since cores will usually be greater than 1, it is strongly recommended that each cluster job have as many cores/slots as the cores argument.

Convergence

It is important to check convergence diagnostics on each model. The hbl_convergence() function returns data frame of summarized convergence diagnostics. max_rhat is the maximum univariate Gelman/Rubin potential scale reduction factor over all the parameters of the model, min_ess_bulk is the minimum bulk effective sample size over the parameters, and min_ess_tail is the minimum tail effective sample size. max_rhat should be below 1.01, and the ESS metrics should both be above 100 times the number of MCMC chains. If any of these conditions are not true, the MCMC did not converge, and it is recommended to try running the model for more saved iterations (and if max_rhat is high, possibly more warmup iterations). You could also try increasing adapt_delta and max_treedepth in the control argument, or adjust other MCMC/HMC settings through control or the informal arguments (ellipsis “...”). See ?rstan::stan for details.

hbl_convergence(mcmc_hierarchical)

#> # A tibble: 1 × 3
#>   max_rhat min_ess_bulk min_ess_tail
#>      <dbl>        <dbl>        <dbl>
#> 1     1.11         8.53         11.9

Results

Each model can be summarized with the hbl_summary() function. The output is a table with few rows and many columns.

summary_hierarchical <- hbl_summary(
  mcmc = mcmc_hierarchical,
  data = data,
  response = "outcome",
  study = "trial",
  study_reference = "trial3",
  group = "arm",
  group_reference = "arm1",
  patient = "subject",
  rep = "visit",
  rep_reference = "visit1",
  covariates = c("factor1", "factor2"),
  eoi = c(0, 1),
  direction = c(">", "<")
)
summary_hierarchical
#> # A tibble: 8 × 54
#>   group group_label   rep rep_label data_n data_N data_n_study_1 data_n_study_2
#>   <dbl> <chr>       <int> <chr>      <int>  <int>          <int>          <int>
#> 1     1 arm1            1 visit1        30     30             10             10
#> 2     1 arm1            2 visit2        30     30             10             10
#> 3     1 arm1            3 visit3        30     30             10             10
#> 4     1 arm1            4 visit4        30     30             10             10
#> 5     2 arm2            1 visit1        10     10              0              0
#> 6     2 arm2            2 visit2        10     10              0              0
#> 7     2 arm2            3 visit3        10     10              0              0
#> 8     2 arm2            4 visit4        10     10              0              0
#> # ℹ 46 more variables: data_n_study_3 <int>, data_N_study_1 <int>,
#> #   data_N_study_2 <int>, data_N_study_3 <int>, data_mean <dbl>, data_sd <dbl>,
#> #   data_lower <dbl>, data_upper <dbl>, response_mean <dbl>,
#> #   response_variance <dbl>, response_sd <dbl>, response_lower <dbl>,
#> #   response_upper <dbl>, response_mean_mcse <dbl>, response_sd_mcse <dbl>,
#> #   response_lower_mcse <dbl>, response_upper_mcse <dbl>, change_mean <dbl>,
#> #   change_lower <dbl>, change_upper <dbl>, change_mean_mcse <dbl>, …

hbl_summary() returns a tidy data frame with one row per group (e.g. treatment arm) and the columns in the following list. Unless otherwise specified, the quantities are calculated at the group-by-rep level. Some are calculated for the current (non-historical) study only, while others pertain to the combined dataset which includes all historical studies.

group: group index.
group_label: original group label in the data.
rep: rep index.
rep_label: original rep label in the data.
data_mean: observed mean of the response specific to the current study.
data_sd: observed standard deviation of the response specific to the current study.
data_lower: lower bound of a simple frequentist 95% confidence interval of the observed data mean specific to the current study.
data_upper: upper bound of a simple frequentist 95% confidence interval of the observed data mean specific to the current study.
data_n: number of non-missing observations in the combined dataset (all studies).
data_N: total number of observations (missing and non-missing) in the combined dataset (all studies).
data_n_study_*: number of non-missing observations in each study. The suffixes of these column names are integer study indexes. Call dplyr::distinct(hbl_data(your_data), study, study_label) to see which study labels correspond to these integer indexes.
data_N_study_*: total number of observations (missing and non-missing) within each study. The suffixes of these column names are integer study indexes. Call dplyr::distinct(hbl_data(your_data), study, study_label) to see which study labels correspond to these integer indexes.
response_mean: Estimated posterior mean of the response from the model. (Here, the response variable in the data should be a change from baseline outcome.) Specific to the current study.
response_sd: Estimated posterior standard deviation of the mean response from the model. Specific to the current study.
response_variance: Estimated posterior variance of the mean response from the model. Specific to the current study.
response_lower: Lower bound of a 95% posterior interval on the mean response from the model. Specific to the current study.
response_upper: Upper bound of a 95% posterior interval on the mean response from the model. Specific to the current study.
response_mean_mcse: Monte Carlo standard error of response_mean.
response_sd_mcse: Monte Carlo standard error of response_sd.
response_lower_mcse: Monte Carlo standard error of response_lower.
response_upper_mcse: Monte Carlo standard error of response_upper.
change_*: same as the response_* columns, but for change from baseline instead of the response. Not included if response_type is "change" because in that case the response is already change from baseline. Specific to the current study.
change_percent_*: same as the change_* columns, but for the percent change from baseline (from 0% to 100%). Not included if response_type is "change" because in that case the response is already change from baseline. Specific to the current study.
diff_*: same as the response_* columns, but for treatment effect.
P(diff > EOI), P(diff < EOI): CSF probabilities on the treatment effect specified with the eoi and direction arguments. Specific to the current study.
effect_mean: same as the response_* columns, but for the effect size (diff / residual standard deviation). Specific to the current study.
precision_ratio*: same as the response_* columns, but for the precision ratio, which compares within-study variance to among-study variance. Only returned for the hierarchical model. Specific to the current study.

Borrowing metrics

The hbl_metrics() function shows borrowing metrics like the mean shift ratio and variance shift ratio which require input from benchmark models. The mean shift ratio and variance shift ratio are defined in the methods vignette. Below, we compute these metrics on the hierarchical model.

summary_pool <- hbl_summary(
  mcmc = mcmc_pool,
  data = data,
  response = "outcome",
  study = "trial",
  study_reference = "trial3",
  group = "arm",
  group_reference = "arm1",
  patient = "subject",
  rep = "visit",
  rep_reference = "visit1",
  covariates = c("factor1", "factor2")
)

summary_independent <- hbl_summary(
  mcmc = mcmc_independent,
  data = data,
  response = "outcome",
  study = "trial",
  study_reference = "trial3",
  group = "arm",
  group_reference = "arm1",
  patient = "subject",
  rep = "visit",
  rep_reference = "visit1",
  covariates = c("factor1", "factor2")
)

hbl_metrics(
  borrow = summary_hierarchical,
  pool = summary_pool,
  independent = summary_independent
)
#> # A tibble: 4 × 4
#>     rep rep_label mean_shift_ratio variance_shift_ratio
#>   <int> <chr>                <dbl>                <dbl>
#> 1     1 visit1               1.70                 0.719
#> 2     2 visit2               0.629                0.665
#> 3     3 visit3              43.2                  0.726
#> 4     4 visit4               0.567                0.676

Plots

The hbl_plot_borrow() function visualizes the results from the hierarchical model against the benchmark models (independent and pooled) to gain intuition about the overall effect of borrowing on estimation.

hbl_plot_borrow(
  borrow = summary_hierarchical,
  pool = summary_pool,
  independent = summary_independent
)

hbl_plot_group() shows the same information but grouped by the group designations in the data (e.g. treatment arm). The results below are not actually correct because the MCMCs ran for so few iterations. For serious analyses, increase the iter and warmup arguments to several thousand and increase the number of chains to about 3 or 4.

hbl_plot_group(
  borrow = summary_hierarchical,
  pool = summary_pool,
  independent = summary_independent
)

For the hierarchical model, we can plot the marginal posterior of the \(\tau\) parameter to see the influence of the uniform prior upper bound. If the probability mass of \(\tau\) builds up against the uniform upper bound, then the estimated hierarchical variance may be larger than assumed, and should indicate that the algorithm wants to borrow less. If the probability mass of \(\tau\) is concentrated on the left, towards 0, then this may indicate that the data suggests borrowing and the uniform prior upper bound has little effect on the model estimates.

hbl_plot_tau(mcmc_hierarchical)

Ongoing problems with borrowing metrics

For some analyses, the mean shift ratio can fall below 0 or above 1. This may occur if pooled model estimates are very close to independent model estimates for some time points, or when the uniform upper bound prior is set to an uninformative large value and pooled model estimates are very far away from independent model estimates, due to random variability in the estimates of the independent, pooled, and hierarchical model estimates. Mean shift ratio issues may also occur when enforcing nontrivial correlation structures among time points because of differences in longitudinal smoothing behavior among models. (In the latter case, specifying diagonal covariance matrices may resolve the issues.)

Similar issues with the variance shift ratio may occur for time points where the pooled model estimate is far away from the independent model, since large shifts in the mean may increase the magnitude of the estimate of residual error for the current study.

The precision ratio metric is based on the conditional estimation of treatment effect, rather than the unconditional estimation, so it will usually be an underestimate of the mean shift ratio or the variance shift ratio. (One benefit of unconditional estimation over conditional estimation is reduced variability in the estimate of treatment effect.)

It is almost always best to set a diffuse prior on \(\tau_t\) so that historical borrowing is fully dynamic. But in extreme cases with small numbers of historical studies, it may be necessary to choose a value of s_tau to approximate the prior amount of borrowing given a desired precision ratio and supposed residual variance. See the function hbl_s_tau() to suggest a rough value to begin exploring. For a large number of historical studies, the computation may take an extremely long time to run. It may be necessary to change arguments covariance_current and covariance_historical (especially the latter) for computational efficiency, but not even that may be enough if the data is enormous.↩︎