Informative Prior Construction based on Historical Data

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)
map_ASAS20 <- gMAP(cbind(r, n-r) ~ 1 | study,
                   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.

map_automix <- automixfit(map_ASAS20)
map_automix

## EM for Beta Mixture Model
## Log-Likelihood = 600.5567
## 
## Univariate beta mixture
## Mixture Components:
##   comp1       comp2      
## w   0.5811855   0.4188145
## a   8.8668236  84.2578948
## b  14.8278270 156.7240118

plot(map_automix)$mix

The resulting MAP prior is approximated by a mixture of conjugate priors, given by π₁(θ) = 0.63Beta(42.5, 77.2) + 0.37Beta(7.2, 12.4), with θ̂_h ≈ 0.36.

SAM Weight Determination

Let θ and θ_h denote the treatment effects associated with the current arm data D and historical D_h, respectively. Let δ denote the clinically significant difference such that is |θ_h − θ| ≥ δ, then θ_h is regarded as clinically distinct from θ, and it is therefore inappropriate to borrow any information from D_h. Consider two hypotheses:

H₀ : θ = θ_h,  H₁ : θ = θ_h + δ or θ = θ_h − δ. H₀ represents that D_h and D are consistent (i.e., no prior-data conflict) and thus information borrowing is desirable, whereas H₁ 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|⋅) 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₀) and P(H₁) is the prior probabilities of H₀ and H₁ being true. BF is the Bayes Factor that in this case is the same as LRT.

The SAM prior, denoted as π_sam(θ), is then defined as a mixture of an informative prior π₁(θ), constructed based on D_h, with a non-informative prior π₀(θ): π_sam(θ) = wπ₁(θ) + (1 − w)π₀(θ) 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, π_sam(θ) 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:

n <- 35; r = 10 
wSAM <- SAM_weight(if.prior = map_automix, 
                   delta = 0.2,
                   n = n, r = r)
cat('SAM weight: ', wSAM)

## SAM weight:  0.7577612

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:

wSAM <- SAM_weight(if.prior = map_automix, 
                   delta = 0.2,
                   method.w = 'PPR',
                   prior.odds = 3/7,
                   n = n, r = r)
cat('SAM weight: ', wSAM)

## SAM weight:  0.5727662

The prior.odds indicates the prior probability of H₀ over the prior probability of H₁. 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.

SAM Prior Construction

To construct the SAM prior, we mix the derived informative prior π₁(θ) with a vague prior π₀(θ) using pre-determined mixture weight by SAM_prior function in SAMprior as follows:

SAM.prior <- SAM_prior(if.prior = map_automix, 
                       nf.prior = mixbeta(nf.prior = c(1,1,1)),
                       weight = wSAM)
SAM.prior

## Univariate beta mixture
## Mixture Components:
##   comp1       comp2       nf.prior   
## w   0.3328834   0.2398828   0.4272338
## a   8.8668236  84.2578948   1.0000000
## b  14.8278270 156.7240118   1.0000000

where the non-informative prior π₀(θ) follows a uniform distribution.

Operating Characteristics

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 (θ_t − θ > 0) > 0.95.

In SAMprior, the operating characteristics can be considered in following steps:

1. Specify priors: This step involves constructing informative prior based on historical data and non-informative prior.

2. 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.

3. 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.

set.seed(123)
TypeI <- get_OC(if.prior = map_automix,       ## MAP prior from historical data
                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)

set.seed(123)
Power <- get_OC(if.prior = map_automix,       ## MAP prior from historical data
                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)

References

[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.

R Session Info

sessionInfo()

## R version 4.4.2 (2024-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] knitr_1.49       SAMprior_1.1.1   ggplot2_3.5.1    Metrics_0.1.4   
## [5] checkmate_2.3.2  assertthat_0.2.1 RBesT_1.7-4      rmarkdown_2.29  
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6         tensorA_0.36.2.1     xfun_0.49           
##  [4] bslib_0.8.0          QuickJSR_1.4.0       inline_0.3.20       
##  [7] vctrs_0.6.5          tools_4.4.2          generics_0.1.3      
## [10] stats4_4.4.2         curl_6.0.1           parallel_4.4.2      
## [13] tibble_3.2.1         fansi_1.0.6          pkgconfig_2.0.3     
## [16] distributional_0.5.0 RcppParallel_5.1.9   lifecycle_1.0.4     
## [19] farver_2.1.2         compiler_4.4.2       stringr_1.5.1       
## [22] munsell_0.5.1        codetools_0.2-20     htmltools_0.5.8.1   
## [25] sys_3.4.3            buildtools_1.0.0     sass_0.4.9          
## [28] bayesplot_1.11.1     yaml_2.3.10          Formula_1.2-5       
## [31] crayon_1.5.3         pillar_1.9.0         jquerylib_0.1.4     
## [34] cachem_1.1.0         StanHeaders_2.32.10  abind_1.4-8         
## [37] posterior_1.6.0      rstan_2.32.6         tidyselect_1.2.1    
## [40] digest_0.6.37        mvtnorm_1.3-2        stringi_1.8.4       
## [43] dplyr_1.1.4          reshape2_1.4.4       maketools_1.3.1     
## [46] labeling_0.4.3       fastmap_1.2.0        grid_4.4.2          
## [49] colorspace_2.1-1     cli_3.6.3            magrittr_2.0.3      
## [52] loo_2.8.0            pkgbuild_1.4.5       utf8_1.2.4          
## [55] withr_3.0.2          scales_1.3.0         backports_1.5.0     
## [58] matrixStats_1.4.1    gridExtra_2.3        evaluate_1.0.1      
## [61] V8_6.0.0             rstantools_2.4.0     rlang_1.1.4         
## [64] Rcpp_1.0.13-1        glue_1.8.0           jsonlite_1.8.9      
## [67] R6_2.5.1             plyr_1.8.9