The self-adapting mixture prior (SAMprior) package is designed to enhance the effectiveness and practicality of clinical trials by leveraging historical information or real-world data [1]. The package incorporate historical data into a new trial using an informative prior constructed based on historical data while mixing a non-informative prior to enhance the robustness of information borrowing. It utilizes a data-driven way to determine a self-adapting mixture weight that dynamically favors the informative (non-informative) prior component when there is little (substantial) evidence of prior-data conflict. Operating characteristics are evaluated and compared to the robust Meta-Analytic-Predictive (rMAP) prior [2], which assigns a fixed weight of 0.5.
Consider a randomized clinical trial to compare a treatment with a control in patients with ankylosing spondylitis. The primary efficacy endpoint is binary, indicating whether a patient achieves 20% improvement at week six according to the Assessment of SpondyloArthritis International Society criteria [3]. Nine historical data available to the control were used to construct the MAP prior:
study | n | r |
---|---|---|
Baeten (2013) | 6 | 1 |
Deodhar (2016) | 122 | 35 |
Deodhar (2019) | 104 | 31 |
Erdes (2019) | 23 | 10 |
Huang (2019) | 153 | 56 |
Kivitz (2018) | 117 | 55 |
Pavelka (2017) | 76 | 28 |
Sieper (2017) | 74 | 21 |
Van der Heijde (2018) | 87 | 35 |
SAM prior is constructed by mixing an informative prior \(\pi_1(\theta)\), constructed based on
historical data, with a non-informative prior \(\pi_0(\theta)\) using the mixture weight
\(w\) determined by
SAM_weight
function to achieve the degree
of prior-data conflict [1]. The following sections describe how to
construct SAM prior in details.
To construct informative priors based on the aforementioned nine
historical data, we apply gMAP
function
from RBesT to perform meta-analysis. This informative prior results in a
representative form from a large MCMC samples, and it can be converted
to a parametric representation with the
automixfit
function using
expectation-maximization (EM) algorithm [4]. This informative prior is
also called MAP prior.
# load R packages
library(ggplot2)
theme_set(theme_bw()) # sets up plotting theme
set.seed(22)
<- gMAP(cbind(r, n-r) ~ 1 | study,
map_ASAS20 family = binomial,
data = ASAS20,
tau.dist = "HalfNormal",
tau.prior = 1,
beta.prior = 2)
## Assuming default prior location for beta: 0
## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess
## Final MCMC sample equivalent to less than 1000 independent draws.
## Please consider increasing the MCMC simulation size.
<- automixfit(map_ASAS20)
map_automix map_automix
## EM for Beta Mixture Model
## Log-Likelihood = 629.3091
##
## Univariate beta mixture
## Mixture Components:
## comp1 comp2
## w 0.530831 0.469169
## a 50.769450 9.059985
## b 89.281035 15.747092
plot(map_automix)$mix
The resulting MAP prior is approximated by a mixture of conjugate priors, given by \(\pi_1(\theta) = 0.63 Beta(42.5, 77.2) + 0.37 Beta(7.2, 12.4)\), with \(\hat{\theta}_h \approx 0.36\).
Let \(\theta\) and \(\theta_h\) denote the treatment effects associated with the current arm data \(D\) and historical \(D_h\), respectively. Let \(\delta\) denote the clinically significant difference such that is \(|\theta_h - \theta| \ge \delta\), then \(\theta_h\) is regarded as clinically distinct from \(\theta\), and it is therefore inappropriate to borrow any information from \(D_h\). Consider two hypotheses:
\[ H_0: \theta = \theta_h, ~~ H_1: \theta = \theta_h + \delta ~ \text{or} ~ \theta = \theta_h - \delta. \] \(H_0\) represents that \(D_h\) and \(D\) are consistent (i.e., no prior-data conflict) and thus information borrowing is desirable, whereas \(H_1\) represents that the treatment effect of \(D\) differs from \(D_h\) to such a degree that no information should be borrowed.
The SAM prior uses the likelihood ratio test (LRT) statistics \(R\) to quantify the degree of prior-data conflict and determine the extent of information borrowing. \[ R = \frac{P(D | H_0, \theta_h)}{P(D | H_1, \theta_h)} = \frac{P(D | \theta = \theta_h)}{\max \{ P(D | \theta = \theta_h + \delta), P(D | \theta = \theta_h - \delta) \}} , \] where \(P(D | \cdot)\) denotes the likelihood function. An alternative Bayesian choice is the posterior probability ratio (PPR): \[ R = \frac{P(D | H_0, \theta_h)}{P(D | H_1, \theta_h)} = \frac{P(H_0)}{P(H_1)} \times BF , \] where \(P(H_0)\) and \(P(H_1)\) is the prior probabilities of \(H_0\) and \(H_1\) being true. \(BF\) is the Bayes Factor that in this case is the same as LRT.
The SAM prior, denoted as \(\pi_{sam}(\theta)\), is then defined as a mixture of an informative prior \(\pi_1(\theta)\), constructed based on \(D_h\), with a non-informative prior \(\pi_0(\theta)\): \[\pi_{sam}(\theta) = w \pi_1(\theta) + (1 - w) \pi_0(\theta)\] where the mixture weight \(w\) is calculated as: \[w = \frac{R}{1 + R}.\] As the level of prior-data conflict increases, the likelihood ratio \(R\) decreases, resulting in a decrease in the weight \(w\) assigned to the informative prior and a decrease in information borrowing. As a result, \(\pi_{sam}(\theta)\) is data-driven and has the ability to self-adapt the information borrowing based on the degree of prior-data conflict.
To calculate mixture weight \(w\) of
the SAM prior, we assume the sample size enrolled in the control arm is
\(n = 35\), with \(r = 10\) responses, then we can apply
function SAM_weight
in SAMprior as
follows:
<- 35; r = 10
n <- SAM_weight(if.prior = map_automix,
wSAM delta = 0.2,
n = n, r = r)
cat('SAM weight: ', wSAM)
## SAM weight: 0.7588881
The default method to calculate \(w\) is using LRT, which is fully data-driven. However, if investigators want to incorporate prior information on prior-data conflict to determine the mixture weight \(w\), this can be achieved by using PPR method as follows:
<- SAM_weight(if.prior = map_automix,
wSAM delta = 0.2,
method.w = 'PPR',
prior.odds = 3/7,
n = n, r = r)
cat('SAM weight: ', wSAM)
## SAM weight: 0.5742702
The prior.odds
indicates the prior
probability of \(H_0\) over the prior
probability of \(H_1\). In this case
(e.g., prior.odds = 3/7
), the prior
information favors the presence prior-data conflict and it results in a
decreased mixture weight.
When historical information is congruent with the current control arm, SAM weight reaches to the highest peak. As the level of prior-data conflict increases, SAM weight decreases. This demonstrates that SAM prior is data-driven and self-adapting, favoring the informative (non-informative) prior component when there is little (substantial) evidence of prior-data conflict.
## Warning: `qplot()` was deprecated in ggplot2 3.4.0.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
To construct the SAM prior, we mix the derived informative prior
\(\pi_1(\theta)\) with a vague prior
\(\pi_0(\theta)\) using pre-determined
mixture weight by SAM_prior
function in
SAMprior as follows:
<- SAM_prior(if.prior = map_automix,
SAM.prior nf.prior = mixbeta(nf.prior = c(1,1,1)),
weight = wSAM)
SAM.prior
## Univariate beta mixture
## Mixture Components:
## comp1 comp2 nf.prior
## w 0.3048404 0.2694298 0.4257298
## a 50.7694503 9.0599849 1.0000000
## b 89.2810355 15.7470916 1.0000000
where the non-informative prior \(\pi_0(\theta)\) follows a uniform distribution.
In this section, we aim to investigate the operating characteristics of the SAM prior, constructed based on the historical data in the context of ankylosing spondylitis trial, via simulation. The incorporation of historical information is expected to be beneficial in reducing the required sample size for the current arms. To achieve this, we assume a 1:2 ratio between the control and treatment arms.
We compare the operating characteristics of the SAM prior and rMAP prior with pre-specified fixed weight under various scenarios. Specifically, we will evaluate the relative bias and relative mean squared error (MSE) of these methods. The relative bias and relative MSE are defined as the differences between the bias/MSE of a given method and the bias/MSE obtained when using a non-informative prior.
Additionally, we investigate the type I error and power of the methods under different degrees of prior-data conflicts. The decision regarding whether a treatment is superior or inferior to a standard control will be based on the condition: \[\Pr(\theta_t - \theta > 0) > 0.95.\]
In SAMprior, the operating characteristics can be considered in following steps:
Specify priors: This step involves constructing informative prior based on historical data and non-informative prior.
Specify the decision criterion: The
decision2S
function is used to initialize
the decision criterion. This criterion determines whether a treatment is
considered superior or inferior to a standard control.
Specify design parameters for the
get_OC
function: This step involves
defining the design parameters for evaluating the operating
characteristics. These parameters include the clinically significant
difference (CSD) used in SAM prior calculation, the method used to
determine the mixture weight for the SAM prior, the sample sizes for the
control and treatment arms, the number of trials used for simulation,
the choice of weight for the robust MAP prior used as a benchmark, and
the vector of response rates for both the control and treatment
arms.
To compute the type I error, we consider four scenarios, with the first and last two scenarios representing minimal and substantial prior-data conflicts, respectively. In general, the results show that both methods effectively control the type I error.
set.seed(123)
<- get_OC(if.prior = map_automix, ## MAP prior from historical data
TypeI nf.prior = mixbeta(c(1,1,1)), ## Non-informative prior for treatment arm
delta = 0.2, ## CSD for SAM prior
## Method to determine the mixture weight for the SAM prior
method.w = 'LRT',
n = 35, n.t = 70, ## Sample size for control and treatment arms
## Decisions
decision = decision2S(0.95, 0, lower.tail=FALSE),
ntrial = 1000, ## Number of trials simulated
if.MAP = TRUE, ## Output robust MAP prior for comparison
weight = 0.5, ## Weight for robust MAP prior
## Response rates for control and treatment arms
theta = c(0.36, 0.36, 0.11, 0.55),
theta.t = c(0.34, 0.33, 0.11, 0.55)
)kable(TypeI)
scenarios | Bias_SAM | Bias_rMAP | RMSE_SAM | RMSE_rMAP | wSAM | res_SAM | res_rMAP | res_NP |
---|---|---|---|---|---|---|---|---|
1 | -0.0040642 | -0.0040324 | -0.0026221 | -0.0034043 | 0.7245092 | 0.020 | 0.015 | 0.026 |
2 | -0.0052884 | -0.0050513 | -0.0028415 | -0.0035849 | 0.7117479 | 0.018 | 0.012 | 0.026 |
3 | 0.0017143 | 0.0172913 | 0.0004427 | 0.0023177 | 0.0245039 | 0.031 | 0.026 | 0.031 |
4 | -0.0132404 | -0.0356093 | 0.0021340 | 0.0025261 | 0.1875810 | 0.098 | 0.107 | 0.046 |
For power evaluation, we also consider four scenarios, with the first and last two scenarios representing minimal and substantial prior-data conflicts, respectively. In general, it is observed that the SAM prior achieves better power compared to rMAP, particularly when there is strong evidence of prior-data conflicts.
set.seed(123)
<- get_OC(if.prior = map_automix, ## MAP prior from historical data
Power nf.prior = mixbeta(c(1,1,1)), ## Non-informative prior for treatment arm
delta = 0.2, ## CSD for SAM prior
n = 35, n.t = 70, ## Sample size for control and treatment arms
## Decisions
decision = decision2S(0.95, 0, lower.tail=FALSE),
ntrial = 1000, ## Number of trials simulated
if.MAP = TRUE, ## Output robust MAP prior for comparison
weight = 0.5, ## Weight for robust MAP prior
## Response rates for control and treatment arms
theta = c(0.37, 0.34, 0.16, 0.11),
theta.t = c(0.57, 0.54, 0.36, 0.31)
)kable(Power)
scenarios | Bias_SAM | Bias_rMAP | RMSE_SAM | RMSE_rMAP | wSAM | res_SAM | res_rMAP | res_NP |
---|---|---|---|---|---|---|---|---|
1 | -0.0067439 | -0.0064999 | -0.0026717 | -0.0033999 | 0.7181020 | 0.799 | 0.771 | 0.622 |
2 | 0.0029869 | 0.0042857 | -0.0023390 | -0.0028376 | 0.7379218 | 0.780 | 0.750 | 0.626 |
3 | 0.0074003 | 0.0278075 | 0.0014730 | 0.0030087 | 0.1098638 | 0.665 | 0.524 | 0.708 |
4 | 0.0020211 | 0.0181728 | 0.0005283 | 0.0026173 | 0.0286888 | 0.748 | 0.587 | 0.766 |
Finally, we present an example of how to make a final decision on
whether the treatment is superior or inferior to a standard control once
the trial has been completed and data has been collected. This step can
be accomplished using the postmix
function
from RBesT, as shown below:
## Sample size and number of responses for treatment arm
<- 70; x_t <- 22
n_t
## first obtain posterior distributions...
<- postmix(priormix = SAM.prior, ## SAM Prior
post_SAM r = r, n = n)
<- postmix(priormix = mixbeta(c(1,1,1)), ## Non-informative prior
post_trt r = x_t, n = n_t)
## Define the decision function
= decision2S(0.95, 0, lower.tail=FALSE)
decision
## Decision-making
decision(post_trt, post_SAM)
## [1] 0
[1] Yang P. et al., Biometrics, 2023; 00, 1–12. https://doi.org/10.1111/biom.13927
[2] Schmidli H. et al., Biometrics 2014; 70(4):1023-1032.
[3] Baeten D. et al., The Lancet, 2013; (382), 9906, p
1705.
[4] Neuenschwander B. et al., Clin Trials. 2010; 7(1):5-18.
sessionInfo()
## R version 4.2.1 (2022-06-23)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur ... 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
##
## locale:
## [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] knitr_1.43 SAMprior_1.1.1 ggplot2_3.4.3 Metrics_0.1.4
## [5] checkmate_2.2.0 assertthat_0.2.1 RBesT_1.7-2
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.11 mvtnorm_1.2-3 prettyunits_1.1.1
## [4] ps_1.7.5 digest_0.6.33 utf8_1.2.3
## [7] plyr_1.8.8 R6_2.5.1 backports_1.4.1
## [10] stats4_4.2.1 evaluate_0.21 highr_0.10
## [13] pillar_1.9.0 rlang_1.1.1 rstudioapi_0.15.0
## [16] callr_3.7.3 jquerylib_0.1.4 rmarkdown_2.24
## [19] labeling_0.4.3 stringr_1.5.0 loo_2.6.0
## [22] munsell_0.5.0 compiler_4.2.1 xfun_0.40
## [25] rstan_2.21.8 pkgconfig_2.0.3 pkgbuild_1.4.2
## [28] rstantools_2.3.1.1 htmltools_0.5.6 tidyselect_1.2.0
## [31] tibble_3.2.1 tensorA_0.36.2 gridExtra_2.3
## [34] codetools_0.2-19 matrixStats_1.0.0 fansi_1.0.4
## [37] crayon_1.5.2 dplyr_1.1.3 withr_2.5.0
## [40] grid_4.2.1 distributional_0.3.2 jsonlite_1.8.7
## [43] gtable_0.3.4 lifecycle_1.0.3 DBI_1.1.3
## [46] magrittr_2.0.3 posterior_1.4.1 StanHeaders_2.26.27
## [49] scales_1.2.1 RcppParallel_5.1.7 stringi_1.7.12
## [52] cli_3.6.1 cachem_1.0.8 reshape2_1.4.4
## [55] farver_2.1.1 bslib_0.5.1 generics_0.1.3
## [58] vctrs_0.6.3 Formula_1.2-5 tools_4.2.1
## [61] glue_1.6.2 processx_3.8.2 abind_1.4-5
## [64] parallel_4.2.1 fastmap_1.1.1 yaml_2.3.7
## [67] inline_0.3.19 colorspace_2.1-0 bayesplot_1.10.0
## [70] sass_0.4.7