This usage tutorial shows how to fit historical control-group
borrowing models using the historicalborrowlong
package.
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.
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.000128 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.28 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.063 seconds (Warm-up)
#> Chain 1: 0.554 seconds (Sampling)
#> Chain 1: 0.617 seconds (Total)
#> Chain 1:
#>
#> SAMPLING FOR MODEL 'historicalborrowlong' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 8.5e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.85 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.082 seconds (Warm-up)
#> Chain 1: 0.095 seconds (Sampling)
#> Chain 1: 0.177 seconds (Total)
#> Chain 1:
#>
#> SAMPLING FOR MODEL 'historicalborrowlong' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 9e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.9 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.066 seconds (Warm-up)
#> Chain 1: 0.071 seconds (Sampling)
#> Chain 1: 0.137 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>, …
The hierarchical model is the primary model of interest.^{1}
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 9.2e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.92 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: 20.577 seconds (Warm-up)
#> Chain 1: 15.529 seconds (Sampling)
#> Chain 1: 36.106 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 1.08 0.960 1.01 0.861 0.438 0.958 0.749
#> 2 1.04 0.981 1.01 0.869 0.159 0.979 0.953
#> 3 1.67 0.605 0.856 0.868 0.156 0.990 1.00
#> 4 1.03 1.00 1.00 0.696 0.647 0.985 0.838
#> 5 1.08 0.922 1.01 0.830 0.626 0.993 1.20
#> 6 1.37 0.925 0.749 0.852 0.416 1.00 1.11
#> 7 1.33 0.780 0.982 1.19 0.114 0.987 1.02
#> 8 1.28 0.907 0.886 0.465 0.685 0.961 0.922
#> 9 1.29 0.952 0.985 1.10 0.495 0.997 1.27
#> 10 1.14 0.864 0.950 0.781 0.555 0.999 0.866
#> # ℹ 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 8.6e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.86 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.915 seconds (Warm-up)
#> Chain 1: 1.229 seconds (Sampling)
#> Chain 1: 4.144 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.06 0.974 1.07 1.13 0.0230 -0.0255 0.00654
#> 2 0.942 0.995 1.01 1.03 -0.340 0.0419 -0.0797
#> 3 1.03 0.995 1.03 0.882 -0.455 -0.0371 -0.424
#> 4 1.11 0.977 1.02 0.894 0.794 1.22 0.911
#> 5 1.07 0.982 1.01 0.917 0.746 1.24 0.957
#> 6 0.946 0.980 1.00 0.981 -0.0591 0.298 0.152
#> 7 0.986 0.978 1.04 0.931 -0.453 0.0943 -0.294
#> 8 0.960 0.977 1.00 0.948 0.332 0.500 0.230
#> 9 1.05 0.974 1.03 0.951 0.933 0.765 0.884
#> 10 1.07 0.986 1.03 0.852 0.419 1.01 0.602
#> # ℹ 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 8.5e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.85 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: 3.202 seconds (Warm-up)
#> Chain 1: 1.196 seconds (Sampling)
#> Chain 1: 4.398 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.35 0.832 0.961 0.960 0.581 0.982 0.983
#> 2 0.817 0.911 0.995 0.616 0.393 0.974 1.01
#> 3 1.02 0.822 1.09 0.891 0.148 0.970 0.956
#> 4 1.18 0.767 0.963 0.860 0.207 0.986 1.29
#> 5 0.941 1.05 0.916 0.679 0.106 0.989 0.810
#> 6 0.958 0.965 1.02 0.639 0.837 0.988 1.40
#> 7 1.40 1.06 0.899 0.621 0.250 0.981 0.983
#> 8 0.779 1.00 1.05 0.752 0.507 0.991 1.27
#> 9 0.792 0.932 1.06 0.815 0.586 0.990 1.10
#> 10 1.33 0.913 1.02 0.881 0.419 0.985 0.871
#> # ℹ 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>, …
A typical workflow will run all three models. Since each model can take several hours to run, it is strongly recommended to:
cores
argument equal to chains
.cores
will usually be greater than 1, it is
strongly recommended that each cluster job have as many cores/slots as
the cores
argument.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.06 18.7 21.5
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.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.31 0.685
#> 2 2 visit2 0.806 0.746
#> 3 3 visit3 2.72 0.755
#> 4 4 visit4 0.661 0.773
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)
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.↩︎