| Title: | Bayesian Dynamic Borrowing with Propensity Score |
|---|---|
| Description: | A tool which aims to help evaluate the effect of external borrowing using an integrated approach described in Lewis et al., (2019) <doi:10.1080/19466315.2018.1497533> that combines propensity score and Bayesian dynamic borrowing methods. |
| Authors: | Isaac Gravestock [cre, ctb], Craig Gower-Page [aut], Matt Secrest [ctb], Yichen Lu [aut], Aijing Lin [aut], F. Hoffmann-La Roche AG [cph, fnd] |
| Maintainer: | Isaac Gravestock <[email protected]> |
| License: | Apache License (>= 2) |
| Version: | 0.2.1 |
| Built: | 2024-10-10 04:15:08 UTC |
| Source: | https://github.com/Genentech/psborrow |
S4 Class for specifying parameters for enrollment time, drop-out pattern and analysis start time
S4 Class for setting parameters for time-to-events
S4 Class for specifying prior distributions and predictors for MCMC methods
Fit a dynamic borrowing Weibull survival model to the given dataset and extract the posterior
samples using MCMC.
See the user guide for more information on the model formulation.
See run_mcmc() for more information on the available parameters for tuning the MCMC sampling
process
apply_mcmc(dt, formula_cov, ...) extract_samples(object) ## S3 method for class 'apply_mcmc' summary(object, ...)apply_mcmc(dt, formula_cov, ...) extract_samples(object) ## S3 method for class 'apply_mcmc' summary(object, ...)
dt |
A data.frame containing data required for modelling. See details |
formula_cov |
A one sided formula specifying which non-treatment covariates should be included into the model. See details |
... |
Additional arguments passed onto |
object |
A |
apply_mcmc()The dt data.frame must contain 1 row per subject with the following variables:
time - A continuous non-zero number specifying the time that the subject had an event at
cnsr - A column of 0/1's where 1 indicates that the event was censored/right truncated
ext - A column of 0/1's where 1 indicates that the subject was part of the external control
trt - A column of 0/1's where 1 indicates that the subject was receiving the experimental treatment
The dt data.frame may also contain any additional covariates to be used in the Weibull model
as specified by formula_cov. In order to fit a valid model formula_cov must contain
the intercept term. The formula will be automatically adjusted to include the treatment term
and as such should not be included here, if you want to include a treatment interaction term
this should be done by using ~ trt:covariate and NOT via ~ trt*covariate.
extract_samples()This function can be used to extract the samples generated by apply_mcmc()
summary()This function provides summary statistics about the samples generated by apply_mcmc()
The extracted samples can be roughly defined as follows (see the user guide for full details):
HR_cc_hc - The hazard ratio between the concurrent control arm and the historical
control arm. This can be
be thought of as the ratio of the scale parameter between the baseline trial distribution and
the baseline
external control distribution. This is equivalent to exp(alpha[2] - alpha[1])
HR_trt_cc - The hazard ratio between the treatment arm and the concurrent control arm.
This is equivalent
to exp(beta_trt)
alpha[1] - The shape parameter for the trial's baseline distribution
alpha[2] - The shape parameter for the historical control's baseline distribution
beta_trt - The log-hazard ratio for the treatment effect. This is equivalent to
log(HR_trt_cc)
beta_<var> - The log-hazard ratio for any other covariate provided to the model via
formula_cov
r0 - The scale parameter for the baseline distribution of both the trial and the
historical control
tau/sigma - The precision/variance for alpha[1] i.e. controls how much
information is borrowed from the
historical control arm
.covClasss classesConcatenate multiple .covClasss classes
## S4 method for signature '.covClass' c(x, ...)## S4 method for signature '.covClass' c(x, ...)
x |
A |
... |
Other |
A vector of .covClasss classes
# combine two sets of covariates covset1 = set_cov(n_cat = 2, n_cont = 0, mu_int = 0, mu_ext = 0, var = 1) covset2 = set_cov(n_cat = 0, n_cont = 1, mu_int = 62, mu_ext = 65, var = 11) cov_list = c(covset1, covset2)# combine two sets of covariates covset1 = set_cov(n_cat = 2, n_cont = 0, mu_int = 0, mu_ext = 0, var = 1) covset2 = set_cov(n_cat = 0, n_cont = 1, mu_int = 62, mu_ext = 65, var = 11) cov_list = c(covset1, covset2)
.priorClasss classConcatenate multiple .priorClasss class
## S4 method for signature '.priorClass' c(x, ...)## S4 method for signature '.priorClass' c(x, ...)
x |
A |
... |
A |
A vector of .priorClasss classes
Utility function to make the mcmc column names more human friendly
fix_col_names(x, column_names)fix_col_names(x, column_names)
x |
a mcmc results object created by |
column_names |
The names to change the beta columns to |
Generate summary statistics of a simulation scenario
get_summary(dt)get_summary(dt)
dt |
a |
a data.frame containing the mean and sd of posterior HR between treatment and control arm, the posterior mean and sd of HR between internal control and external control arm, reject rate, variance, bias and mse of the simulation set
Simple function which leverages the DESCRIPTION file to check if the user is in a development environment for psborrow.
is_psborrow_dev()is_psborrow_dev()
TRUE/FALSE flag (TRUE = in development environment)
Match
match_cov(dt, match)match_cov(dt, match)
dt |
a list of |
match |
A vector of covariates name to match on |
a list of matrix containing matched cohort information
# match internal and external trial data using different covariates smp = set_n(ssC = 140, ssE = 275, ssExt = 100) covset1 = set_cov(n_cat = 2, n_cont = 0, mu_int = 0, mu_ext = 0, var = 1) covset2 = set_cov(n_cat = 0, n_cont = 1, mu_int = 62, mu_ext = 65, var = 11) cObj = c(covset1, covset2) sample_cov <- simu_cov(ssObj = smp, covObj = cObj, HR = 1, driftHR = 1.2, nsim = 2) # match on covariates 1 and 2 match_cov(dt = sample_cov, match = c("cov1", "cov2")) # match on all 3 covariates match_cov(dt = sample_cov, match = c("cov1", "cov2", "cov3"))# match internal and external trial data using different covariates smp = set_n(ssC = 140, ssE = 275, ssExt = 100) covset1 = set_cov(n_cat = 2, n_cont = 0, mu_int = 0, mu_ext = 0, var = 1) covset2 = set_cov(n_cat = 0, n_cont = 1, mu_int = 62, mu_ext = 65, var = 11) cObj = c(covset1, covset2) sample_cov <- simu_cov(ssObj = smp, covObj = cObj, HR = 1, driftHR = 1.2, nsim = 2) # match on covariates 1 and 2 match_cov(dt = sample_cov, match = c("cov1", "cov2")) # match on all 3 covariates match_cov(dt = sample_cov, match = c("cov1", "cov2", "cov3"))
Plot bias for each prior distribution according to selected simulation parameters
plot_bias(dt, HR = 1, driftHR = 1, pred = "none")plot_bias(dt, HR = 1, driftHR = 1, pred = "none")
dt |
a |
HR |
pre-specified HR between treatment and control arm in the internal trial for which the bias should be plotted. Must be within |
driftHR |
pre-specified HR between external control arm and internal control arm for which the bias should be plotted. Must be within |
pred |
predictors to use when fitting exponential distribution in MCMC for which the bias should be plotted. Must be within |
a bar plot of class ggplot containing the bias for each prior distribution.
Plot mean posterior hazard ratio between treatment and control
plot_hr(dt, HR = 0.67, driftHR = 1, pred = "none")plot_hr(dt, HR = 0.67, driftHR = 1, pred = "none")
dt |
a |
HR |
pre-specified HR between treatment and control arm in the internal trial for which the mean posterior hazard ratio should be plotted. Must be within |
driftHR |
pre-specified HR between external control arm and internal control arm for which the mean posterior hazard ratio should be plotted. Must be within |
pred |
predictors to use when fitting exponential distribution in MCMC for which the mean posterior hazard ratio should be plotted. Must be within |
a plot of class ggplot containing the mean posterior hazard ratio for each prior distribution.
Plot mean squared error (MSE) for each prior distribution according to selected simulation parameters
plot_mse(dt, HR = 1, driftHR = 1, pred = "none")plot_mse(dt, HR = 1, driftHR = 1, pred = "none")
dt |
a |
HR |
pre-specified HR between treatment and control arm in the internal trial for which the MSE should be plotted. Must be within |
driftHR |
pre-specified HR between external control arm and internal control arm for which the MSE should be plotted. Must be within |
pred |
predictors to use when fitting exponential distribution in MCMC for which the MSE should be plotted. Must be within |
a bar plot of class ggplot containing the MSE for each prior distribution.
Plot power for each prior distribution according to selected simulation parameters
plot_power(dt, HR = 0.67, driftHR = 1, pred = "none")plot_power(dt, HR = 0.67, driftHR = 1, pred = "none")
dt |
a |
HR |
pre-specified HR between treatment and control arm in the internal trial for which the power should be plotted. Must be within |
driftHR |
pre-specified HR between external control arm and internal control arm for which the power should be plotted. Must be within |
pred |
predictors to use when fitting exponential distribution in MCMC for which the power should be plotted. Must be within |
a bar plot of class ggplot containing the power for each prior distribution.
Plot type 1 error for each prior distribution according to selected simulation parameters
plot_type1error(dt, driftHR = 1, pred = "none")plot_type1error(dt, driftHR = 1, pred = "none")
dt |
a |
driftHR |
the driftHR between the external and internal control arms for which the type 1 error should be plotted. Must be within |
pred |
the predictors used when fitting the exponential distribution in MCMC for which the type 1 error should be plotted. Must be within |
a bar plot of class ggplot containing type 1 error proportions for each prior distribution.
Simple wrapper function around message() that will supress
printing messages if the option psborrow.quiet is set to TRUE
i.e.
options("psborrow.quiet" = TRUE)
ps_message(...)ps_message(...)
... |
Values passed onto |
Generate summary statistics for the MCMC chains
rej_est(samples)rej_est(samples)
samples |
an object of class |
a vector containing the mean, median, sd, reject rate for the MCMC chains
Run MCMC for multiple scenarios with provided data
run_mcmc(dt, priorObj, n.chains, n.adapt, n.burn, n.iter, seed, path)run_mcmc(dt, priorObj, n.chains, n.adapt, n.burn, n.iter, seed, path)
dt |
a list of |
priorObj |
an object of class |
n.chains |
number of parallel chains for the model |
n.adapt |
number of iterations for adaptation |
n.burn |
number of iterations discarded as burn-in |
n.iter |
number of iterations to monitor |
seed |
the seed of random number generator. Default is the first element of .Random.seed |
path |
file name for saving the output including folder path |
a data.frame containing summary statistics of the posterior distribution for each simulation
# examples in vignette# examples in vignette
Run MCMC for multiple scenarios with provided data with parallel processing
run_mcmc_p( dt, priorObj, n.chains, n.adapt, n.burn, n.iter, seed, path, n.cores = 2 )run_mcmc_p( dt, priorObj, n.chains, n.adapt, n.burn, n.iter, seed, path, n.cores = 2 )
dt |
a list of |
priorObj |
an object of class |
n.chains |
number of parallel chains for the model |
n.adapt |
number of iterations for adaptation |
n.burn |
number of iterations discarded as burn-in |
n.iter |
number of iterations to monitor |
seed |
the seed of random number generator. Default is the first element of .Random.seed |
path |
file name for saving the output including folder path |
n.cores |
number of processes to parallelize over (default = 2) |
a data.frame containing summary statistics of the posterior distribution for each simulation
# similar to run_mcmc# similar to run_mcmc
This function allows user to specify the enrollment and drop-out rate, and the type of clinical cut-off Date. Both enrollment times and drop-out times follow piece-wise exponential distribution.
set_clin(gamma, e_itv, CCOD, CCOD_t, etaC, etaE, d_itv)set_clin(gamma, e_itv, CCOD, CCOD_t, etaC, etaE, d_itv)
gamma |
A vector of rate of enrollment per unit of time |
e_itv |
A vector of duration of time periods for recruitment with rates specified in |
CCOD |
Type of analysis start time. Analysis starts at |
CCOD_t |
Time difference between analysis start and first patient's enrollment if |
etaC |
A vector for dropout rate per unit time for control arm |
etaE |
A vector for dropout rate per unit time for experimental arm. If left |
d_itv |
A vector of duration of time periods for dropping out the study with rates specified in |
A .clinClass class containing information on enrollment time, drop-out pattern and analysis start time
# set the operational parameter values for the trial # analysis starts at64 time units after first patient in set_clin(gamma = 10, e_itv = 4, etaC = 0.003, CCOD = "fixed-first", CCOD_t = 64) # analysis starts at 12 time units after last patient in set_clin(gamma = 2, e_itv = 18, etaC = 0.005, CCOD = "fixed-last", CCOD_t = 12)# set the operational parameter values for the trial # analysis starts at64 time units after first patient in set_clin(gamma = 10, e_itv = 4, etaC = 0.003, CCOD = "fixed-first", CCOD_t = 64) # analysis starts at 12 time units after last patient in set_clin(gamma = 2, e_itv = 18, etaC = 0.005, CCOD = "fixed-last", CCOD_t = 12)
This function saves the mean, variance and covariance among covariates. For technical details, see the vignette.
set_cov(n_cat, n_cont, mu_int, mu_ext, var, cov, prob_int, prob_ext)set_cov(n_cat, n_cont, mu_int, mu_ext, var, cov, prob_int, prob_ext)
n_cat |
Number of binary variable. See details |
n_cont |
Number of continuous variable |
mu_int |
Mean of covariates in the internal trial. All the covariates are simulated from a
multivariate normal distribution. If left |
mu_ext |
Mean of covariates in the external trial. If left |
var |
Variance of covariates. If left |
cov |
Covariance between each pair of covariates. Covariance needs to be provided in
a certain order and users are encouraged to read the example provided in the vignette. If
left |
prob_int |
Probability of binary covariate equalling 1 in the internal trial. If left
|
prob_ext |
Probability of binary covariate equalling 1 in the external trial. If
left |
Categorical variables are created by sampling a continuous variable from the multivariate
normal
distribution (thus respecting the correlation to other covariates specified by cov)
and then applying a cut point derived from the prob_int or prob_ext quantile
of said distribution i.e. for a univariate variable it would be derived as:
binvar <- as.numeric(rnorm(n, mu, sqrt(var)) < qnorm(prob, mu, sqrt(var)))
Please note that this means that the value of mu_int & mu_ext has no impact on categorical
covariates and thus can be set to any value.
As an example of how this process works assume n_cat=3 and n_cont=2. First 5 variables are
sampled from the multivariate normal distribution as specified by mu_int/mu_ext, var &
cov. Then, the first 3 of these variables are converted to binary based on the probabilities
specified by prob_int and prob_ext. This means that that the 2 continuous variables will
take their mean and sd from the last 2 entries in the vectors mu_int/mu_ext and var.
A .covClass class containing covariate information
Defines the model formula and distribution to be used when simulating time-to-events. Please see the user-guide for the model formulations
set_event(event, lambdaC, beta, shape, t_itv, change, keep)set_event(event, lambdaC, beta, shape, t_itv, change, keep)
event |
Distribution of time-to-events: |
lambdaC |
Baseline hazard rate of internal control arm. Specify a vector for piece-wise
hazard with duration specified in |
beta |
covariates' coefficients (i.e. log hazard ratios). Must be equal in length to the number of covariates
created by |
shape |
the shape parameter of Weibull distribution if |
t_itv |
a vector indicating interval lengths where the exponential rates provided in
|
change |
A list of additional derivered covariates to be used in simulating time-to-events. See details |
keep |
A character vector specifying which of the original covariates (i.e. those not
derived via the |
The change argument is used to specify additional derived covariates to be used when
simulating time-to-events. For example, let’s say have 3 covariates cov1, cov2 & cov3
but that we also wish to include a new covariate that is an interaction
between cov1 and cov2 as well as another covariate that is equal to the sum of
cov2 and cov3; we could implement this as follows:
set_event(
event = "weibull",
shape = 0.9,
lambdaC = 0.0135,
beta = c(5, 3, 1, 7, 9),
change = list(
c("cov1", "*", "cov2"),
c("cov2", "+", "cov3")
)
)
Note that in the above example 5 values have been specified to beta,
3 for the original three covariates
and 2 for the two additional derived covariates included via change.
Variables derived via change are automatically included in the model regardless
of whether they are listed in keep or not. Likewise, these covariates are derived
separately and not via a standard R formula, that is to say including an interaction
term does not automatically include the individual fixed effects.
a .eventClass class containing time-to-events information
a matrix containing simulated time-to-events information
# time-to-event follows a Weibull distribution set_event(event = "weibull", shape = 0.9, lambdaC = 0.0135) # time-to-event follows a piece-wise exponential distribution set_event(event = "pwexp", t_itv = 1, lambdaC = c(0.1, 0.02))# time-to-event follows a Weibull distribution set_event(event = "weibull", shape = 0.9, lambdaC = 0.0135) # time-to-event follows a piece-wise exponential distribution set_event(event = "pwexp", t_itv = 1, lambdaC = c(0.1, 0.02))
This function conducts validity check and generates a matrix with two binary variables indicating
if the observation belongs to the external trial
if the observation belongs to the treatment arm.
set_n(ssC, ssE, ssExt)set_n(ssC, ssE, ssExt)
ssC |
Number of observations in the internal control arm. Default is 100 |
ssE |
Number of observations in the internal experiment arm. Default is the same number
of observations as |
ssExt |
Number of observations in the external control arm. Default is the same number
of observations as |
A matrix containing external trial indicator and treatment indicator
Specify prior distributions and predictors for MCMC methods
set_prior(pred, prior, r0, alpha, sigma)set_prior(pred, prior, r0, alpha, sigma)
pred |
Predictors to include in the weibull distribution. No covariates except for
treatment indicator is included if |
prior |
Prior distribution for the precision parameter that controls the degree of
borrowing. Half-cauchy distribution if |
r0 |
Initial values for the shape of the weibull distribution for time-to-events |
alpha |
Initial values for log of baseline hazard rate for external and internal control
arms. Length of |
sigma |
Initial values for precision parameter if |
a .priorClass class containing survival data and prior information
# hierachical Bayesian model with precision parameter follows a half-cauchy distribution set_prior(pred = "none", prior = "cauchy", r0 = 1, alpha = c(0, 0), sigma = 0.03) # hierachical Bayesian model with precision parameter follows a gamma distribution set_prior(pred = "none", prior = "gamma", r0 = 1, alpha = c(0, 0)) # conventional Bayesian model to not borrow from external control arm set_prior(pred = "none", prior = "no_ext", alpha = 0) # conventional Bayesian model to fully borrow from external control arm set_prior(pred = "none", prior = "full_ext", alpha = 0)# hierachical Bayesian model with precision parameter follows a half-cauchy distribution set_prior(pred = "none", prior = "cauchy", r0 = 1, alpha = c(0, 0), sigma = 0.03) # hierachical Bayesian model with precision parameter follows a gamma distribution set_prior(pred = "none", prior = "gamma", r0 = 1, alpha = c(0, 0)) # conventional Bayesian model to not borrow from external control arm set_prior(pred = "none", prior = "no_ext", alpha = 0) # conventional Bayesian model to fully borrow from external control arm set_prior(pred = "none", prior = "full_ext", alpha = 0)
This function generates continuous and binary covariates through simulating from a multivariate normal distribution. Outcomes are further converted to binary variables using quantiles of the normal distribution calculated from the probability provided. Then the covariates are added to the external trial and treatment arm indicators.
simu_cov(ssObj, covObj, driftHR, HR, nsim, seed, path)simu_cov(ssObj, covObj, driftHR, HR, nsim, seed, path)
ssObj |
an object of class |
covObj |
an object of class |
driftHR |
hazard ratio of external control and internal control arms |
HR |
a list of hazard ratio of treatment and control arms |
nsim |
number of simulation. Default is 5 |
seed |
the seed of R‘s random number generator. Default is the first element of .Random.seed |
path |
file name for saving the output including folder path |
a list of matrix containing simulated covariates information
# simulate patient-level data with 1 continuous covariate sample = set_n(ssC = 10, ssE = 20, ssExt = 40) cov1 = set_cov(n_cat = 0, n_cont = 1, mu_int = 0, mu_ext = 0, var = 1) simu_cov(ssObj = sample, covObj = cov1, HR = 0.5, driftHR = 1, nsim = 2) # simulate patient-level data with 1 binary and 2 continuous covariate cov2 = set_cov(n_cat = 1, n_cont = 2, mu_int = 0, mu_ext = 0, var = 1, cov = 0.3, prob_int = 0.2, prob_ext = 0.3) simu_cov(ssObj = sample, covObj = cov2, HR = 0.5, driftHR = 1, nsim = 2)# simulate patient-level data with 1 continuous covariate sample = set_n(ssC = 10, ssE = 20, ssExt = 40) cov1 = set_cov(n_cat = 0, n_cont = 1, mu_int = 0, mu_ext = 0, var = 1) simu_cov(ssObj = sample, covObj = cov1, HR = 0.5, driftHR = 1, nsim = 2) # simulate patient-level data with 1 binary and 2 continuous covariate cov2 = set_cov(n_cat = 1, n_cont = 2, mu_int = 0, mu_ext = 0, var = 1, cov = 0.3, prob_int = 0.2, prob_ext = 0.3) simu_cov(ssObj = sample, covObj = cov2, HR = 0.5, driftHR = 1, nsim = 2)
Simulate time-to-events for multiple scenarios
simu_time(dt, eventObj, clinInt, clinExt, seed, path)simu_time(dt, eventObj, clinInt, clinExt, seed, path)
dt |
a list of |
eventObj |
an object of class |
clinInt |
an object of class |
clinExt |
an object of class |
seed |
the seed of R‘s random number generator. Default is the first element of .Random.seed |
path |
file name for saving the output including folder path |
a list of matrix containing simulated time-to-events information
# simulate patient-level data without covariates # simulate survival time following weibull distribution # simulate trial indicator and set hazard ratios sample = set_n(ssC = 10, ssE = 20, ssExt = 40) sample_hr <- simu_cov(ssObj = sample, HR = 1, driftHR=c(1,1.2), nsim = 10) # enrollment pattern, drop-out, analysis start time c_int = set_clin(gamma = 2, e_itv = 10, etaC = 0.5, CCOD = "fixed-first", CCOD_t = 64) c_ext = c_int # simulate time-to-event with a weibull distribution evt1 <- set_event(event = "weibull", shape = 0.8, lambdaC = 0.01) simu_time(dt = sample_hr, eventObj = evt1, clinInt = c_int, clinExt = c_ext) # simulate time-to-event with an exponential distribution evt2 <- set_event(event = "pwexp", t_itv = 1, lambdaC = c(0.1, 0.02)) simu_time(dt = sample_hr, eventObj = evt2, clinInt = c_int, clinExt = c_int)# simulate patient-level data without covariates # simulate survival time following weibull distribution # simulate trial indicator and set hazard ratios sample = set_n(ssC = 10, ssE = 20, ssExt = 40) sample_hr <- simu_cov(ssObj = sample, HR = 1, driftHR=c(1,1.2), nsim = 10) # enrollment pattern, drop-out, analysis start time c_int = set_clin(gamma = 2, e_itv = 10, etaC = 0.5, CCOD = "fixed-first", CCOD_t = 64) c_ext = c_int # simulate time-to-event with a weibull distribution evt1 <- set_event(event = "weibull", shape = 0.8, lambdaC = 0.01) simu_time(dt = sample_hr, eventObj = evt1, clinInt = c_int, clinExt = c_ext) # simulate time-to-event with an exponential distribution evt2 <- set_event(event = "pwexp", t_itv = 1, lambdaC = c(0.1, 0.02)) simu_time(dt = sample_hr, eventObj = evt2, clinInt = c_int, clinExt = c_int)