The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

Analysis Example of Bayesian MCPMod for Continuous Data

library(BayesianMCPMod)
library(clinDR)
library(dplyr)

set.seed(7015)

Background and data

In this vignette we will show the use of the Bayesian MCPMod package for continuous distributed data. The focus lies on the utilization of an informative prior and the Bayesian MCPMod evaluation of a single trial. We will use data that is included in the clinDR package. More specifically, trial results of BRINTELLIX will be used to illustrate the specification of an informative prior and the usage of such a prior for the Bayesian evaluation of a (hypothetical) new trial. BRINTELLIX is a medication used to treat major depressive disorder. Various clinical trials with different dose groups, including control groups, were conducted.

Calculation of a MAP prior

In a first step, a meta analytic prior will be calculated using historical data from 4 trials with main endpoint Change from baseline in MADRS score after 8 weeks. Please note that only information from the control group will be integrated leading to an informative mixture prior for the control group, while for the active groups a non-informative prior will be specified.

data("metaData")
dataset     <- filter(as.data.frame(metaData), bname == "BRINTELLIX")
histcontrol <- filter(
  dataset,
  dose       == 0,
  primtime   == 8,
  indication == "MAJOR DEPRESSIVE DISORDER",
  protid     != 5)

hist_data   <- data.frame(
  trial = histcontrol$nctno,
  est   = histcontrol$rslt,
  se    = histcontrol$se,
  sd    = histcontrol$sd,
  n     = histcontrol$sampsize)

Here, we suggest a function to construct a list of prior distributions for the different dose groups. This function is adapted to the needs of this example. Other applications may need a different way to construct prior distributions.

getPriorList <- function (
  
  hist_data,
  dose_levels,
  dose_names    = NULL,
  robust_weight = 0.5
  
) {
  
  sd_tot <- with(hist_data, sum(sd * n) / sum(n))
  
  gmap <- RBesT::gMAP(
    formula    = cbind(est, se) ~ 1 | trial,
    weights    = hist_data$n,
    data       = hist_data,
    family     = gaussian,
    beta.prior = cbind(0, 100 * sd_tot),
    tau.dist   = "HalfNormal",
    tau.prior  = cbind(0, sd_tot / 4))
  
  prior_ctr <- RBesT::automixfit(gmap)
  
  if (!is.null(robust_weight)) {
    
    prior_ctr <- suppressMessages(RBesT::robustify(
      priormix = prior_ctr,
      weight   = robust_weight,
      sigma    = sd_tot))
    
  }
  
  prior_trt <- RBesT::mixnorm(
    comp1 = c(w = 1, m = summary(prior_ctr)[1], n = 1),
    sigma = sd_tot,
    param = "mn")
  
  prior_list <- c(list(prior_ctr),
                  rep(x     = list(prior_trt),
                      times = length(dose_levels[-1])))
  
  if (is.null(dose_names)) {
    
    dose_names <- c("Ctr", paste0("DG_", seq_along(dose_levels[-1])))
    
  }
  
  names(prior_list) <- dose_names
  
  return (prior_list)
  
}

With the dose levels to be investigated, the prior distribution can be constructed.

dose_levels <- c(0, 2.5, 5, 10)

prior_list  <- getPriorList(
  hist_data     = hist_data,
  dose_levels   = dose_levels,
  robust_weight = 0.3)

getESS(prior_list)
#>  Ctr DG_1 DG_2 DG_3 
#> 20.1  1.0  1.0  1.0

Specifications for the new trial

To be able to apply the Bayesian MCPMod approach, candidate models need to be specified using functions from the R package DoseFinding. Since there are only 3 active dose levels we will limit the set of candidate models to a linear, an exponential, and an emax model. Note that the linear candidate model does not require a guesstimate.

exp_guesst  <- DoseFinding::guesst(
  d     = 5,
  p     = c(0.2),
  model = "exponential",
  Maxd  = max(dose_levels))

emax_guesst <- DoseFinding::guesst(
  d     = 2.5,
  p     = c(0.9),
  model = "emax")

mods <- DoseFinding::Mods(
  linear      = NULL,
  emax        = emax_guesst,
  exponential = exp_guesst,
  doses       = dose_levels,
  maxEff      = -1,
  placEff     = -12.8)

We will use the trial with ct.gov number NCT00735709 as exemplary new trial.

new_trial  <- filter(
  dataset,
  primtime   == 8,
  indication == "MAJOR DEPRESSIVE DISORDER",
  protid     == 5)

n_patients <- c(128, 124, 129, 122)

Combination of prior information and trial results

As outlined in Fleischer et al. [Bayesian MCPMod. Pharmaceutical Statistics. 2022; 21(3): 654-670.], in a first step the posterior is calculated combining the prior information with the estimated results of the new trial. Via the summary function it is possible to print out the summary information of the posterior distributions.

posterior <- getPosterior(
  prior    = prior_list,
  mu_hat   = new_trial$rslt,
  se_hat   = new_trial$se,
  calc_ess = TRUE)

summary(posterior)
#>           mean        sd      2.5%     50.0%      97.5%
#> Ctr  -11.17511 0.7149157 -12.57452 -11.16897  -9.789677
#> DG_1 -14.88115 0.7130817 -16.27877 -14.88115 -13.483539
#> DG_2 -15.08016 0.7101057 -16.47195 -15.08016 -13.688383
#> DG_3 -15.63661 0.7259755 -17.05950 -15.63661 -14.213724

Execution of Bayesian MCPMod Test step

For the execution of the testing step of Bayesian MCPMod a critical value on the probability scale will be determined for a given alpha level. This critical value is calculated using the re-estimated contrasts for the frequentist MCPMod to ensure that, when using non-informative priors, the actual error level for falsely declaring a significant trial in the Bayesian MCPMod is controlled by the specified alpha level.

A pseudo-optimal contrast matrix is generated based on the variability of the posterior distribution (see Fleischer et al. 2022 for more details).

crit_pval <- getCritProb(
  mods           = mods,
  dose_levels    = dose_levels,
  se_new_trial   = new_trial$se,
  alpha_crit_val = 0.05)

contr_mat <- getContr(
  mods         = mods,
  dose_levels  = dose_levels,
  sd_posterior = summary(posterior)[, 2])

Please note that there are different ways to derive the contrasts. The following code shows the implementation of some of these ways but it is not executed and the contrast specification above is used.

# i) the frequentist contrast
contr_mat_prior <- getContr(
  mods           = mods,
  dose_levels    = dose_levels,
  dose_weights   = n_patients,
  prior_list     = prior_list)
# ii) re-estimated frequentist contrasts
contr_mat_prior <- getContr(
  mods           = mods,
  dose_levels    = dose_levels,
  se_new_trial   = new_trial$se)
# iii)  Bayesian approach using number of patients for new trial and prior distribution
contr_mat_prior <- getContr(
  mods           = mods,
  dose_levels    = dose_levels,
  dose_weights   = n_patients,
  prior_list     = prior_list)

The Bayesian MCP testing step is then executed based on the posterior information, the provided contrasts and the multiplicity adjusted critical value.

BMCP_result <- performBayesianMCP(
  posterior_list = posterior,
  contr          = contr_mat, 
  crit_prob_adj  = crit_pval)

BMCP_result
#> Bayesian Multiple Comparison Procedure
#>      sign crit_prob_adj max_post_prob post_probs.linear post_probs.emax
#> [1,]    1     0.9740992     0.9999998         0.9999529       0.9999998
#>      post_probs.exponential
#> [1,]              0.9980234
#> attr(,"successRate")
#> [1] 1
#> attr(,"ess_avg")
#>   Ctr  DG_1  DG_2  DG_3 
#> 186.1 186.6 188.2 180.0

The testing step is significant indicating a non-flat dose-response shape. In detail, all 3 models are significant and the p-value for the emax model indicates deviation from the null hypothesis the most.

Model fitting and visualization of results

In the model fitting step the posterior distribution is used as basis. Both simplified and full fitting are performed. For the simplified fit, the multivariate normal distribution of the control group is approximated and reduced by a one-dimensional normal distribution. The actual fit (on this approximated posterior distribution) is then performed using generalized least squares criterion. In contrast, for the full fit, the non-linear optimization problem is addressed via the Nelder Mead algorithm.

The output of the fit includes information about the predicted effects for the included dose levels, the generalized AIC, and the corresponding weights. For the considered case, the simplified and the full fit are very similar.

# Option a) Simplified approach by using approximated posterior distribution
fit_simple <- getModelFits(
  models      = names(mods),
  dose_levels = dose_levels,
  posterior   = posterior,
  simple      = TRUE)

# Option b) Making use of the complete posterior distribution
fit <- getModelFits(
  models      = names(mods),
  dose_levels = dose_levels,
  posterior   = posterior,
  simple      = FALSE)

Via the predict() function, one can also receive estimates for dose levels that were not included in the trial.

predict(fit, doses = c(0, 2.5, 4, 5, 7, 10))
#> $linear
#> [1] -12.34765 -13.36060 -13.96837 -14.37355 -15.18391 -16.39945
#> 
#> $emax
#> [1] -11.21518 -14.80555 -15.13389 -15.25711 -15.40776 -15.52834
#> 
#> $exponential
#> [1] -12.52970 -13.31742 -13.83954 -14.21002 -15.00894 -16.36759
#> 
#> attr(,"doses")
#> [1]  0.0  2.5  4.0  5.0  7.0 10.0

It is possible to plot the fitted dose response models and an AIC based average model (black lines).

plot(fit_simple)

plot(fit)

To assess the uncertainty, one can additionally visualize credible bands (orange shaded areas, default levels are 50% and 95%). These credible bands are calculated with a bootstrap method as follows: Samples from the posterior distribution are drawn and for every sample the simplified fitting step and a prediction is performed. These predictions are then used to identify and visualize the specified quantiles.

plot(fit, cr_bands = TRUE)

The bootstrap based quantiles can also be directly calculated via the getBootstrapQuantiles() function. For this example, only 6 quantiles are bootstrapped for each model fit.

bs_sample <- getBootstrapSamples(
  model_fits = fit,
  doses      = c(0, 2.5, 4, 5, 7, 10),
  n_samples  = 6)

getBootstrapQuantiles(
  bs_sample = bs_sample,
  quantiles = c(0.025, 0.5, 0.975))
#>          model dose      2.5%       50%     97.5%
#> 1       linear  0.0 -13.68887 -12.76192 -11.90959
#> 2       linear  2.5 -14.09215 -13.64127 -12.86736
#> 3       linear  4.0 -14.33412 -14.16889 -13.44203
#> 4       linear  5.0 -14.53047 -14.46881 -13.82514
#> 5       linear  7.0 -15.25634 -14.91039 -14.56200
#> 6       linear 10.0 -16.34513 -15.91349 -15.25640
#> 7         emax  0.0 -12.19232 -11.25725 -10.30861
#> 8         emax  2.5 -15.50039 -14.96948 -14.43386
#> 9         emax  4.0 -15.50758 -15.09123 -14.61432
#> 10        emax  5.0 -15.50998 -15.14971 -14.67588
#> 11        emax  7.0 -15.52479 -15.25516 -14.72865
#> 12        emax 10.0 -15.53704 -15.33852 -14.76919
#> 13 exponential  0.0 -13.82903 -12.97193 -12.11521
#> 14 exponential  2.5 -14.11136 -13.63587 -12.85005
#> 15 exponential  4.0 -14.29850 -14.07595 -13.33711
#> 16 exponential  5.0 -14.43446 -14.37550 -13.68272
#> 17 exponential  7.0 -15.07893 -14.77887 -14.40640
#> 18 exponential 10.0 -16.24617 -15.86902 -15.15993
#> 19      avgFit  0.0 -12.19232 -11.25725 -10.30861
#> 20      avgFit  2.5 -15.50039 -14.96948 -14.43386
#> 21      avgFit  4.0 -15.50758 -15.09123 -14.61432
#> 22      avgFit  5.0 -15.50998 -15.14971 -14.67588
#> 23      avgFit  7.0 -15.52479 -15.25516 -14.72865
#> 24      avgFit 10.0 -15.53704 -15.33852 -14.76919

Technical note: The median quantile of the bootstrap based procedure is not necessary similar to the main model fit, as they are derived via different procedures. The main fit, i.e. the black lines in the plot, show the best fit of a certain model based on minimizing the residuals for the posterior distribution, while the bootstrap based 50% quantile shows the median fit of the random sampling and fitting procedure.

Additional note

It is also possible to perform the testing and modeling step in a combined fashion via the performBayesianMCPMod() function. This code serves merely as an example and is not run in this vignette.

performBayesianMCPMod(
      posterior_list   = posterior,
      contr            = contr_mat, 
      crit_prob_adj    = crit_pval,
      simple           = FALSE)

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.