# Introduction

The purpose of this vignette is to introduce how to simulate and analyze an adaptive Bayesian clinical trial for continuous-valued, i.e., normal, outcomes. The simulation section compromises the design of the trial itself which provides type I error rate and power at each interim look. We use the normal conjugate prior for the estimation of posterior values. Available historical data can be used as an informative prior; we use the bayesDP package as the engine for incorporating the historical data. Instead of using traditional R function, we use pipes to input our parameters.

# Estimation of Treatment Difference

Let $$\bar{y}$$, $$s$$, and $$N$$ denote the sample mean, sample standard deviation, and sample size, respectively. Then, the posterior distribution of the mean under vague (flat) priors is

$\begin{array}{rcl} \tilde{\sigma}^2\mid\bar{y},s,N & \sim & InverseGamma\left(\frac{N-1}{2},\,\frac{N-1}{2}s^2 \right),\\ \\ \tilde{\mu}\mid\bar{y},N,\tilde{\sigma}^2 & \sim & \mathcal{N}ormal\left(\bar{y},\, \frac{1}{N}\tilde{\sigma}^2 \right). \end{array}$

When historical data is present, let $$\bar{y}_0$$, $$s_0$$, and $$N_0$$ denote the sample mean, sample standard deviation, and sample size of the historical data, respectively. The posterior distribution of the mean for the historical data is $\begin{array}{rcl} \sigma^2_0\mid\bar{y_0},s_0,N_0 & \sim & InverseGamma\left(\frac{N_0-1}{2},\,\frac{N_0-1}{2}s_0^2 \right),\\ \\ \mu_0 \mid \bar{y}_0, N_0, \sigma^2_0 & \sim & \mathcal{N}ormal\left(\bar{y}_0,\, \frac{1}{N_0}\sigma^2_0 \right). \end{array}$ The weight of the historical data included in the study design and analysis is denoted by $$\hat\alpha$$. For more details on computation of $$\hat{\alpha}$$, please refer to the vignette of binomial counts available at https://CRAN.R-project.org/package=bayesDP. The posterior distribution of the mean outcome with historical data incorporated under vague (flat) priors is $\tilde{\mu} \sim \mathcal{N}ormal\left( \frac{\sigma^2_0N\bar{y} + \tilde{\sigma}^2N_0\bar{y}_0\hat{\alpha}}{N\sigma^2_0 + \tilde{\sigma}^2N_0\hat{\alpha}},\,\frac{\tilde{\sigma}^2\sigma^2_0}{N\sigma^2_0 + \tilde{\sigma}^2N_0\hat{\alpha}} \right).$

Even though there is a closed-form solution for the difference in normally distributed random variables, we use Monte Carlo simulations to estimate the treatment difference.

The estimation of the treatment difference is $$\tilde{\mu_T} - \tilde{\mu_C}$$, where $$\mu_T$$ is the posterior mean outcome in the treatment group and $$\mu_C$$ is the posterior mean outcome in the control group.

# Wrapper Function for Design and Analysis

Unlike traditional R functions, the bayesCT package depends on pipe inputs with different wrapper functions. All the wrapper functions are illustrated below along with details on the arguments for the simulation and analysis.

• normal_outcome - wrapper function for specifying observations. Inputs:
• mu_treatment - mean of the treatment group. Required input.
• sd_treatment - standard deviation of the treatment group. Required input.
• mu_control - mean of the control group. Default value is NULL.
• sd_control - standard deviation of the control group. Default value is NULL.
• study_details - wrapper function for specifying sample size, study length, interim looks, and loss to follow-up. Inputs:
• total_sample_size - sample size required to run the simulation
• study_period - length of the study, i.e., the follow-up period
• interim_look - enrollment values where interim looks are carried out. Specify where the trials are stopped to check for early success or futility. For more details on the early stopping for futility or success, please see https://thevaachandereng.github.io/bayesCT/articles/bayesCT.html#early-stopping-for-futility-or-success. Each interim_look value must be smaller than total_sample_size.
• prop_loss_to_followup - proportion of subjects loss to follow-up. Must be a value be 0 and 1. Default value is 0.10.
• enrollment_rate - wrapper function for specifying enrollment. Inputs:
• hypothesis - wrapper function for both simulation and analysis. Inputs:
• randomize - wrapper function for specifying the randomization scheme. Inputs:
• impute - wrapper function specifying imputation of outcomes at the interim looks. Inputs:
• no_of_impute - number of imputations. Used to impute outcomes for subjects potentially loss to follow-up when estimating the probability of stopping for early success and/or futility. In the case of the futility calculation, subjects not yet enrolled or imputed as well. Default value is 10 so that the vignette can be generated quickly. However, much larger values should be used in practice, e.g., 10000.
• number_mcmc - number of Monte Carlo iterations for sampling from posterior distributions during imputation. Default value is 10000.
• historical_normal - wrapper function for specifying historical data, if available. This function should only be used when historical data is incorporated in the analysis or simulation. For more details on the method and computation, please see https://CRAN.R-project.org/package=bayesDP. Inputs:
• mu0_treatment - mean of the historical treatment group.
• sd0_treatment - standard deviation of the historical treatment group.
• N0_treatment - sample size of the historical treatment group.
• mu0_control - mean of the historical control group.
• sd0_control - standard deviation of the historical control group.
• N0_control - sample size of the historical control group.
• discount_function - discount function to use for controlling the weight given to the historical data. Currently supports weibull, scaledweibull, and identity. The discount function scaledweibull scales the output of the Weibull CDF to have a max value of 1. The identity discount function uses the posterior probability directly as the discount weight. Default value is "identity".
• alpha_max - maximum weight the discount function can apply. Default is 1. For a two-arm trial, users may specify a vector of two values where the first value is used to weight the historical treatment group and the second value is used to weight the historical control group.
• fix_alpha - Should alpha be set to alpha_max? Default value is FALSE.
• weibull_scale - scale parameter of the Weibull discount function used to compute alpha, the weight parameter of the historical data. Default value is 0.135. For a two-arm trial, users may specify a vector of two values where the first value is used to estimate the weight of the historical treatment group and the second value is used to estimate the weight of the historical control group. Not used when discount_function = "identity".
• weibull_shape - shape parameter of the Weibull discount function used to compute alpha, the weight parameter of the historical data. Default value is 3. For a two-arm trial, users may specify a vector of two values where the first value is used to estimate the weight of the historical treatment group and the second value is used to estimate the weight of the historical control group. Not used when discount_function = "identity".
• method - analysis method with respect to estimation of the weight paramter alpha. Default method mc estimates alpha for each Monte Carlo iteration. Alternate value fixed estimates alpha once and holds it fixed throughout the analysis. See the the bdpsurvival vignette vignette("bdpnormal-vignette", package="bayesDP" for more details.
• data_normal - wrapper function for inputting data with a normal outcome. Inputs:
• treatment - assignment for patients in the trial, 1 for treatment group and 0 for control group. Input should be a vector.
• outcome - continuous valued outcome of the trial. Input should be a vector with same length as treatment.
• complete - indicator for whether outcome is observed; 1 for complete outcome, 0 for loss to follow up or not yet observed. If a value is 0, the outcome of the patient is ignored. If not provided, the dataset is assumed to be complete.
• analysis - wrapper function for analyzing a trial. Inputs:
• input - list of input to analyze the adaptive trial.
• type - outcome type of the trial. Use type="normal" for a trial with normal outcomes.
• simulate - wrapper function for simulating trials. Should be used as the terminal end of the pipe. Inputs:
• input - list of inputs specifying the adaptive trial simulation set-up.
• no_of_sim - number of trials to simulate. Default value is 10000.

In the following section, we will discuss the design of adaptive trials using bayesCT for normal outcomes. We illustrate an example for one-arm trial and two-arm trials using the wrapper functions described above.

## One-arm Trial

In the example below, we will illustrate how to compute power, type 1 error, and other clinical trial characteristics for an objective performance criterion (OPC) trial with mean outcome and hypothesis described as follows, $H_0: \mu_{treatment} \geq 120 \qquad H_A:\mu_{treatment} < 120.$

The most important wrapper functions are study_details and normal_outcome (especially since there are no default values).

The normal mean outcomes are simulated using a mean value of 120 and standard deviation of 5.5. The total sample size is 400 with a study length of 60 days. A 10% loss to follow up is assumed. Based on this information, the adaptive trials are simulated 10 times to obtain the following output (NOTE: for the purpose of reproducing the vignette quickly, we reduce the number of simulations to 10, you should use a much larger value, e.g., 10000). The aforementioned inputs were chosen for illustration purposes only.

value <- normal_outcome(mu_treatment = 120,
sd_treatment = 5.5) %>%
study_details(total_sample_size     = 400,
study_period          = 60,
interim_look          = NULL,
prop_loss_to_followup = 0.10)

# Simulate 10 trials
output <- value %>%
simulate(no_of_sim = 10)

# Structure of the simulation output
str(output)
#> List of 8
#>  $input :List of 5 #> ..$ mu_treatment         : num 120
#>   ..$sd_treatment : num 5.5 #> ..$ N_total              : num 400
#>   ..$EndofStudy : num 60 #> ..$ prop_loss_to_followup: num 0.1
#>  $power :'data.frame': 1 obs. of 2 variables: #> ..$ interim_looks: num 400
#>   ..$power : num 0 #>$ type1_error                 : num 0
#>  $est_final : num [1:10] 120 120 120 120 120 ... #>$ post_prob_accept_alternative: num [1:10] 0.513 0.48 0.471 0.483 0.51 ...
#>  $N_enrolled : num [1:10] 400 400 400 400 400 400 400 400 400 400 #>$ stop_expect_success         : num [1:10] 0 0 0 0 0 0 0 0 0 0
#>  $stop_futility : num [1:10] 0 0 0 0 0 0 0 0 0 0 To allow for early stopping for success or futility, we add interim looks to the design. We’ll check for success or futility at the enrollment of the 350th and 380th subject. Upon adding this interim look requirement, the trial is simulated 10 times to obtain the output. # adding interim look value <- value %>% study_details(total_sample_size = 400, study_period = 60, interim_look = c(350, 380), prop_loss_to_followup = 0.10) # Simulate 10 trials output <- value %>% simulate(no_of_sim = 10) # Structure of the simulation output str(output) #> List of 8 #>$ input                       :List of 6
#>   ..$mu_treatment : num 120 #> ..$ sd_treatment         : num 5.5
#>   ..$N_total : num 400 #> ..$ EndofStudy           : num 60
#>   ..$prop_loss_to_followup: num 0.1 #> ..$ interim_look         : num [1:2] 350 380
#>  $power :'data.frame': 3 obs. of 2 variables: #> ..$ interim_looks: num [1:3] 350 380 400
#>   ..$power : num [1:3] 0 0 0 #>$ type1_error                 : num 0
#>  $est_final : num [1:10] 120 120 120 120 120 ... #>$ post_prob_accept_alternative: num [1:10] 0.497 0.514 0.479 0.53 0.484 ...
#>  $N_enrolled : num [1:10] 350 350 350 350 350 350 350 350 350 350 #>$ stop_expect_success         : num [1:10] 0 0 0 0 0 0 0 0 0 0
#>  $stop_futility : num [1:10] 1 1 1 1 1 1 1 1 1 1 Patient enrollment is assumed to follow a Poisson process. The default enrollment rate is 0.3 patients per day. In this simulation we’ll introduce a step-wise Poisson process with rate $$\lambda$$ as follows: $\lambda = \left\{ \begin{array}{ll} 0.4 & \text(time) \in [0, 40) \\ 0.7 & \text(time) \in [40, \infty) \\ \end{array} \right.$ This enrollment scheme is illustrated below. value <- value %>% enrollment_rate(lambda = c(0.4, 0.7), time = 40) output <- value %>% simulate(no_of_sim = 10) str(output) #> List of 8 #>$ input                       :List of 8
#>   ..$mu_treatment : num 120 #> ..$ sd_treatment         : num 5.5
#>   ..$N_total : num 400 #> ..$ EndofStudy           : num 60
#>   ..$prop_loss_to_followup: num 0.1 #> ..$ interim_look         : num [1:2] 350 380
#>   ..$lambda : num [1:2] 0.4 0.7 #> ..$ lambda_time          : num 40
#>  $power :'data.frame': 3 obs. of 2 variables: #> ..$ interim_looks: num [1:3] 350 380 400
#>   ..$power : num [1:3] 0 0 0 #>$ type1_error                 : num 0
#>  $est_final : num [1:10] 120 120 120 119 121 ... #>$ post_prob_accept_alternative: num [1:10] 0.489 0.53 0.502 0.456 0.555 ...
#>  $N_enrolled : num [1:10] 350 350 350 350 350 350 350 350 350 350 #>$ stop_expect_success         : num [1:10] 0 0 0 0 0 0 0 0 0 0
#>  $stop_futility : num [1:10] 1 1 1 1 1 1 1 1 1 1 The hypothesis is an important wrapper function which controls the probability of futility, probability of accepting the alternative hypothesis, probability of early success, the alternative hypothesis, and the treatment difference margin. Since, in an OPC trial, the outcomes in the treatment group are simulated using the input provided, delta controls the maximum threshold allowed for trial to succeed/fail. The default value of delta is 0. Here, we’ll use delta = -10 (i.e $$120 - \hat{\mu}_{treatment} < -10$$). We’ll further set the futility probability to 0.10, the expected success probability for early stopping to 0.85, and the final probability of accepting the alternative to 0.95. The alternative is "less" due to the hypothesis function specified above. value <- value %>% hypothesis(delta = -10, futility_prob = 0.10, prob_accept_ha = 0.95, expected_success_prob = 0.85, alternative = "less") output <- value %>% simulate(no_of_sim = 10) str(output) #> List of 8 #>$ input                       :List of 13
#>   ..$mu_treatment : num 120 #> ..$ sd_treatment         : num 5.5
#>   ..$N_total : num 400 #> ..$ EndofStudy           : num 60
#>   ..$prop_loss_to_followup: num 0.1 #> ..$ interim_look         : num [1:2] 350 380
#>   ..$lambda : num [1:2] 0.4 0.7 #> ..$ lambda_time          : num 40
#>   ..$h0 : num -10 #> ..$ futility_prob        : num 0.1
#>   ..$prob_ha : num 0.95 #> ..$ expected_success_prob: num 0.85
#>   ..$alternative : chr "less" #>$ power                       :'data.frame':    3 obs. of  2 variables:
#>   ..$interim_looks: num [1:3] 350 380 400 #> ..$ power        : num [1:3] 1 1 1
#>  $type1_error : num 0 #>$ est_final                   : num [1:10] 120 120 121 120 120 ...
#>  $post_prob_accept_alternative: num [1:10] 0.964 0.963 0.956 0.967 0.954 ... #>$ N_enrolled                  : num [1:10] 350 350 350 350 350 350 350 350 350 350
#>  $stop_expect_success : num [1:10] 1 1 1 1 1 1 1 1 1 1 #>$ stop_futility               : num [1:10] 0 0 0 0 0 0 0 0 0 0

Next, we’ll illustrate imputations for imputing outcomes for subjects loss to follow up. We’ll carry out 20 imputations and draw 2000 values from the posterior of each imputation.

value <- value %>%
impute(no_of_impute = 20,
number_mcmc  = 2000)

output <- value %>%
simulate(no_of_sim = 10)

str(output)
#> List of 8
#>  $input :List of 15 #> ..$ mu_treatment         : num 120
#>   ..$sd_treatment : num 5.5 #> ..$ N_total              : num 400
#>   ..$EndofStudy : num 60 #> ..$ prop_loss_to_followup: num 0.1
#>   ..$interim_look : num [1:2] 350 380 #> ..$ lambda               : num [1:2] 0.4 0.7
#>   ..$lambda_time : num 40 #> ..$ h0                   : num -10
#>   ..$futility_prob : num 0.1 #> ..$ prob_ha              : num 0.95
#>   ..$expected_success_prob: num 0.85 #> ..$ alternative          : chr "less"
#>   ..$N_impute : num 20 #> ..$ number_mcmc          : num 2000
#>  $power :'data.frame': 3 obs. of 2 variables: #> ..$ interim_looks: num [1:3] 350 380 400
#>   ..$power : num [1:3] 1 1 1 #>$ type1_error                 : num 0
#>  $est_final : num [1:10] 119 120 120 120 120 ... #>$ post_prob_accept_alternative: num [1:10] 0.976 0.965 0.969 0.962 0.963 ...
#>  $N_enrolled : num [1:10] 350 350 350 350 350 350 350 350 350 350 #>$ stop_expect_success         : num [1:10] 1 1 1 1 1 1 1 1 1 1
#>  $stop_futility : num [1:10] 0 0 0 0 0 0 0 0 0 0 The above flow was for illustrative purposes. Instead of inputting parameters step by step, the trial parameters can be filled in all at once as illustrated below. The pipe function connects all inputs together and the trial is simulated 20 times to obtain results. value <- normal_outcome(mu_treatment = 120, sd_treatment = 5.5) %>% study_details(total_sample_size = 400, study_period = 60, interim_look = c(350, 380), prop_loss_to_followup = 0.10) %>% hypothesis(delta = -10, futility_prob = 0.10, prob_accept_ha = 0.95, expected_success_prob = 0.85, alternative = "less") %>% enrollment_rate(lambda = c(0.4, 0.7), time = 4) %>% randomize(block_size = c(10, 20), randomization_ratio = c(1, 1)) %>% impute(no_of_impute = 20, number_mcmc = 2000) str(value) #> List of 17 #>$ mu_treatment         : num 120
#>  $sd_treatment : num 5.5 #>$ N_total              : num 400
#>  $EndofStudy : num 60 #>$ interim_look         : num [1:2] 350 380
#>  $prop_loss_to_followup: num 0.1 #>$ h0                   : num -10
#>  $futility_prob : num 0.1 #>$ prob_ha              : num 0.95
#>  $expected_success_prob: num 0.85 #>$ alternative          : chr "less"
#>  $lambda : num [1:2] 0.4 0.7 #>$ lambda_time          : num 4
#>  $block : num [1:2] 10 20 #>$ rand_ratio           : num [1:2] 1 1
#>  $N_impute : num 20 #>$ number_mcmc          : num 2000

## Two-arm Trial

In this section, we will illustrate how to perform the design of a two-arm trial with the incorporation of historical data. The example will compute the type 1 error, power, and other outputs for a superiority trial. The study hypothesis is $H_0: \mu_{treatment} - \mu_{control} \leq 0 \qquad H_A: \mu_{treatment} - \mu_{control} > 0.$

Unlike the OPC trial above, we will not include interim looks. The normal mean outcomes are simulated using a mean value of 13 and standard deviation of 1.4 for the treatment group and a mean value of 15 and standard deviation of 1.9 for the control group. The total sample size is 400 with a study length of 30 days. A 15% loss to follow up is assumed. The following code simulates a trial 10 times using the piping procedure.

value <- normal_outcome(mu_treatment = 13,
mu_control = 16,
sd_treatment = 1.4,
sd_control = 1.9) %>%
study_details(total_sample_size     = 300,
study_period          = 50,
interim_look          = NULL,
prop_loss_to_followup = 0.10) %>%
hypothesis(delta                 = 0,
futility_prob         = 0,
prob_accept_ha        = 0.95,
expected_success_prob = 1,
alternative           = "less") %>%
impute(no_of_impute = 25,
number_mcmc  = 5000) %>%
enrollment_rate(lambda = c(0.8),
time = NULL) %>%
randomize(block_size          = c(4, 6),
randomization_ratio = c(1, 1)) %>%
historical_normal(mu0_treatment     = 13,
sd0_treatment     = 5,
N0_treatment      = 100,
mu0_control       = 12,
sd0_control       = 3,
N0_control        = 120,
discount_function = "scaledweibull",
alpha_max         = FALSE,
fix_alpha         = 1,
weibull_scale     = 0.135,
weibull_shape     = 3,
method            = "fixed") %>%
simulate(no_of_sim = 10)

str(value)
#> List of 8
#>  $input :List of 29 #> ..$ mu_control           : num 16
#>   ..$sd_control : num 1.9 #> ..$ mu_treatment         : num 13
#>   ..$sd_treatment : num 1.4 #> ..$ N_total              : num 300
#>   ..$EndofStudy : num 50 #> ..$ prop_loss_to_followup: num 0.1
#>   ..$h0 : num 0 #> ..$ futility_prob        : num 0
#>   ..$prob_ha : num 0.95 #> ..$ expected_success_prob: num 1
#>   ..$alternative : chr "less" #> ..$ N_impute             : num 25
#>   ..$number_mcmc : num 5000 #> ..$ lambda               : num 0.8
#>   ..$block : num [1:2] 4 6 #> ..$ rand_ratio           : num [1:2] 1 1
#>   ..$mu0_treatment : num 13 #> ..$ sd0_treatment        : num 5
#>   ..$N0_treatment : num 100 #> ..$ mu0_control          : num 12
#>   ..$sd0_control : num 3 #> ..$ N0_control           : num 120
#>   ..$discount_function : chr "scaledweibull" #> ..$ alpha_max            : logi FALSE
#>   ..$fix_alpha : num 1 #> ..$ weibull_scale        : num 0.135
#>   ..$weibull_shape : num 3 #> ..$ method               : chr "fixed"
#>  $power :'data.frame': 1 obs. of 2 variables: #> ..$ interim_looks: num 300
#>   ..$power : num 1 #>$ type1_error                 : num 0.3
#>  $est_final : num [1:10] -3.03 -2.77 -2.83 -3.44 -2.72 ... #>$ post_prob_accept_alternative: num [1:10] 1 1 1 1 1 1 1 1 1 1
#>  $N_enrolled : num [1:10] 300 300 300 300 300 300 300 300 300 300 #>$ stop_expect_success         : num [1:10] 0 0 0 0 0 0 0 0 0 0
#>  $stop_futility : num [1:10] 0 0 0 0 0 0 0 0 0 0 # Analysis In this section, we will demonstrate how to run an adaptive Bayesian trial using bayesCT. A sample dataset is provided in the package. The dataset normaldata contains the results of 300 subjects from a two-arm trial with normal outcome. The complete column indicates whether the outcome was observed, i.e., loss to follow-up. data(normaldata) head(normaldata) #> id treatment outcome complete #> 1 Patient1 1 13.702191 1 #> 2 Patient2 1 11.925518 1 #> 3 Patient3 1 10.586943 1 #> 4 Patient4 1 10.450031 1 #> 5 Patient5 1 7.154615 1 #> 6 Patient6 1 8.874396 0 The minimum input needed to run an adaptive Bayesian trial is the data itself. The data_normal input allows the input of the data. The treatment group (0 for control, 1 for treatment) and outcome input are essential for the analysis. However, if the complete input is not provided, the function assumes the outcome data is complete. A default analysis is carried out below. input <- data_normal(treatment = normaldata$treatment,
outcome   = normaldata$outcome, complete = normaldata$complete)

out <- input %>%
analysis(type = "normal")

str(out)
#> List of 11
#>  $prob_of_accepting_alternative: num 0.95 #>$ margin                       : num 0
#>  $alternative : chr "greater" #>$ N_treatment                  : num 139
#>  $N_control : int 127 #>$ N_complete                   : num 266
#>  $N_enrolled : int 300 #>$ post_prob_accept_alternative : num 0.138
#>  $est_final : num -3.05 #>$ stop_futility                : num 1
#>  $stop_expected_success : num 0 We’ll now illustrate using piping to carry out the complete analysis. First, we’ll assume the following hypothesis: $H_0:\mu_{treatment} - \mu_{control} > 0 \quad H_A: \mu_{treatment} - \mu_{control} < 0$ The delta and alternative used to analyze the trial is 0 and “less” respectively. The probability of accepting the alternative is 0.95, the probability of stopping for futility is 0.05, and the probability of stopping for success is 0.90. We will carry out imputations on subjects loss to follow up. Additionally, we will incorporate historical data on the treatment arm. out <- data_normal(treatment = normaldata$treatment,
outcome   = normaldata$outcome, complete = normaldata$complete) %>%
hypothesis(delta = 0,
futility_prob         = 0.05,
prob_accept_ha        = 0.95,
expected_success_prob = 0.90,
alternative           = "less") %>%
impute(no_of_impute = 40,
number_mcmc  = 8000) %>%
analysis(type = "normal")

str(out)
#> List of 11
#>  $prob_of_accepting_alternative: num 0.95 #>$ margin                       : num 0
#>  $alternative : chr "less" #>$ N_treatment                  : num 139
#>  $N_control : int 127 #>$ N_complete                   : num 266
#>  $N_enrolled : int 300 #>$ post_prob_accept_alternative : num 0.866
#>  $est_final : num -3.09 #>$ stop_futility                : num 1
#>  \$ stop_expected_success        : num 0