Binomial Outcome

Thevaa Chandereng, Donald Musgrove, Tarek Haddad, Graeme Hickey, Timothy Hanson, Theodore Lystig

Introduction

The purpose of this vignette is to introduce how to simulate and analyze an adaptive Bayesian clinical trial for binomial 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 beta-binomial conjugate prior for the estimation of posterior probabilities. Available historical data can be used as an informative prior; we use the bayesDP package as the engine for incorporating the historical data. By default, the model uses a non-informative prior of \(Beta(a = 1, b = 1)\) with or without the incorporation of historical data. Instead of using traditional R function, we use pipes to input our parameters.

Estimation of Treatment Difference

Let \(y\) and \(N\) denote the number of events and the sample size, respectively. Let \(a_0\) and \(b_0\) denote the rate parameters of a Beta distribution. Then, the posterior distribution of the event rate under vague (flat) priors is

\[ \tilde{\theta}\mid y,N \; \sim \; \mathcal{B}eta\left(y+a_0,\,N-y+b_0 \right).\] When historical data is present, \(y_0\) and \(N_0\) denote the number of events and sample size of the historical data. The weight of the historical data included in the study design and analysis is denoted by \(\hat\alpha\). For more details on the computation of \(\hat{\alpha}\), please refer to the vignette of binomial counts avaialable at https://CRAN.R-project.org/package=bayesDP. The posterior distribution of the event rate with historical data incorporated under vague (flat) priors is

\[\tilde{\theta} \mid y,N, y_0, N_0 \; \sim \; \mathcal{B}eta\left(y+y_0\hat{\alpha}+a_0,\, N-y+\hat{\alpha}(N_0-y_0)+b_0 \right),\] Since there is no closed-form solution for the difference in beta distributed random variables, we use Monte Carlo simulations to estimate the posterior of the treatment difference.

The estimation of the difference in proportions is \(\tilde{\theta_T} - \tilde{\theta_C}\), where \(\theta_T\) is the posterior of the event rates in the treatment group and \(\theta_C\) is the posterior of the event rates in the control group.

The following section lays out each of the functions and inputs for carrying out simulations and analyses of Bayesian adaptive trials.

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.

Design of Adaptive Trials

In the following section, we will discuss the design of adaptive trials using bayesCT for binomial 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 characteristics for an objective performace criterion (OPC) trial with an observed proportion of events and hypothesis described as follows, \[H_0: \theta_{treatment} \geq 0.08 \qquad H_A:\theta_{treatment} < 0.08.\]

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

Binomial events are simulated using an event rate of 0.08. The total sample size is 900 with a study length of 50 days. A 10% loss to follow-up ais 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 <- binomial_outcome(p_treatment = 0.08) %>%
   study_details(total_sample_size     = 900, 
                 study_period          = 50,
                 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 4
#>   ..$ p_treatment          : num 0.08
#>   ..$ N_total              : num 900
#>   ..$ EndofStudy           : num 50
#>   ..$ prop_loss_to_followup: num 0.1
#>  $ power                       :'data.frame':    1 obs. of  2 variables:
#>   ..$ interim_looks: num 900
#>   ..$ power        : num 0
#>  $ type1_error                 : num 0.1
#>  $ est_final                   : num [1:10] 0.0813 0.0863 0.0801 0.0875 0.0838 ...
#>  $ post_prob_accept_alternative: num [1:10] 0.542 0.727 0.491 0.768 0.64 ...
#>  $ N_enrolled                  : num [1:10] 900 900 900 900 900 900 900 900 900 900
#>  $ 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 can add interim looks to the design. We’ll check for success or futility at the enrollment of the 600th, 700th and 800th subject. Upon adding this interim look requirement, the trial is simulated 10 times to obtain the output.

# Adding interim looks
value <- value %>%
  study_details(total_sample_size     = 900, 
                study_period          = 50,
                interim_look          = c(600, 700, 800),
                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
#>   ..$ p_treatment          : num 0.08
#>   ..$ N_total              : num 900
#>   ..$ EndofStudy           : num 50
#>   ..$ prop_loss_to_followup: num 0.1
#>   ..$ interim_look         : num [1:3] 600 700 800
#>  $ power                       :'data.frame':    4 obs. of  2 variables:
#>   ..$ interim_looks: num [1:4] 600 700 800 900
#>   ..$ power        : num [1:4] 0 0 0 0
#>  $ type1_error                 : num 0.1
#>  $ est_final                   : num [1:10] 0.0817 0.0767 0.0799 0.0879 0.0936 ...
#>  $ post_prob_accept_alternative: num [1:10] 0.541 0.367 0.477 0.752 0.911 ...
#>  $ N_enrolled                  : num [1:10] 600 600 600 700 900 600 800 800 600 600
#>  $ 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 0 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.3 & \text(time) \in [0, 25) \\ 1 & \text(time) \in [25, \infty) \\ \end{array} \right. \]

This enrollment scheme is illustrated below.

value <- value %>%
  enrollment_rate(lambda = c(0.3, 1), 
                  time   = 25)

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

str(output)
#> List of 8
#>  $ input                       :List of 7
#>   ..$ p_treatment          : num 0.08
#>   ..$ N_total              : num 900
#>   ..$ EndofStudy           : num 50
#>   ..$ prop_loss_to_followup: num 0.1
#>   ..$ interim_look         : num [1:3] 600 700 800
#>   ..$ lambda               : num [1:2] 0.3 1
#>   ..$ lambda_time          : num 25
#>  $ power                       :'data.frame':    4 obs. of  2 variables:
#>   ..$ interim_looks: num [1:4] 600 700 800 900
#>   ..$ power        : num [1:4] 0 0 0 0
#>  $ type1_error                 : num 0
#>  $ est_final                   : num [1:10] 0.0514 0.0748 0.0886 0.0904 0.094 ...
#>  $ post_prob_accept_alternative: num [1:10] 0.0034 0.3021 0.7868 0.7969 0.8724 ...
#>  $ N_enrolled                  : num [1:10] 600 700 800 600 600 600 600 600 600 600
#>  $ 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, probablity of early success, the alternative hypothesis, and the treatment difference margin.

Since, in an OPC trial, the proportion of events in the treatment group are simulated using the input provided, delta controls the maximum threshold allowed for the trial to succeed/fail. The default value of delta is 0. Here, we’ll use delta = -0.03 (i.e \(\hat\theta_{treatment} - 0.08 > -0.03\)).

We’ll further set the futility probability to 0.05, the expected success probability for early stopping to 0.90, 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                 = -0.03, 
             futility_prob         = 0.05, 
             prob_accept_ha        = 0.95,
             expected_success_prob = 0.90, 
             alternative           = "less")

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

str(output)
#> List of 8
#>  $ input                       :List of 12
#>   ..$ p_treatment          : num 0.08
#>   ..$ N_total              : num 900
#>   ..$ EndofStudy           : num 50
#>   ..$ prop_loss_to_followup: num 0.1
#>   ..$ interim_look         : num [1:3] 600 700 800
#>   ..$ lambda               : num [1:2] 0.3 1
#>   ..$ lambda_time          : num 25
#>   ..$ h0                   : num -0.03
#>   ..$ futility_prob        : num 0.05
#>   ..$ prob_ha              : num 0.95
#>   ..$ expected_success_prob: num 0.9
#>   ..$ alternative          : chr "less"
#>  $ power                       :'data.frame':    4 obs. of  2 variables:
#>   ..$ interim_looks: num [1:4] 600 700 800 900
#>   ..$ power        : num [1:4] 0.6 0.7 0.9 0.9
#>  $ type1_error                 : num 0
#>  $ est_final                   : num [1:10] 0.0863 0.0757 0.0801 0.0706 0.0881 ...
#>  $ post_prob_accept_alternative: num [1:10] 0.983 0.997 0.992 0.999 0.967 ...
#>  $ N_enrolled                  : num [1:10] 800 600 600 600 700 800 600 900 600 600
#>  $ stop_expect_success         : num [1:10] 1 1 1 1 1 1 1 0 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 25 imputations and draw 1000 values from the posterior of each imputation.

value <- value %>%
  impute(no_of_impute = 25, 
         number_mcmc  = 1000)

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

str(output)
#> List of 8
#>  $ input                       :List of 14
#>   ..$ p_treatment          : num 0.08
#>   ..$ N_total              : num 900
#>   ..$ EndofStudy           : num 50
#>   ..$ prop_loss_to_followup: num 0.1
#>   ..$ interim_look         : num [1:3] 600 700 800
#>   ..$ lambda               : num [1:2] 0.3 1
#>   ..$ lambda_time          : num 25
#>   ..$ h0                   : num -0.03
#>   ..$ futility_prob        : num 0.05
#>   ..$ prob_ha              : num 0.95
#>   ..$ expected_success_prob: num 0.9
#>   ..$ alternative          : chr "less"
#>   ..$ N_impute             : num 25
#>   ..$ number_mcmc          : num 1000
#>  $ power                       :'data.frame':    4 obs. of  2 variables:
#>   ..$ interim_looks: num [1:4] 600 700 800 900
#>   ..$ power        : num [1:4] 0.3 0.4 0.7 0.9
#>  $ type1_error                 : num 0
#>  $ est_final                   : num [1:10] 0.091 0.0938 0.0886 0.0809 0.0776 ...
#>  $ post_prob_accept_alternative: num [1:10] 0.965 0.94 0.971 0.997 0.997 0.982 0.993 0.996 1 0.985
#>  $ N_enrolled                  : num [1:10] 900 900 900 800 800 600 600 800 600 700
#>  $ stop_expect_success         : num [1:10] 0 0 0 1 1 1 1 1 1 1
#>  $ stop_futility               : num [1:10] 0 0 0 0 0 0 0 0 0 0

The default non-informative beta prior used in the simulation is \(\mathcal{B}eta(1, 1)\). In our OPC trial simulation, we’ll change the default to \(\mathcal{B}eta(5,5)\). This will increase the weight of the non-informative prior in the simulation. This non-informative beta prior is implemented using beta_prior wrapper function.

value <- value %>%
  beta_prior(a0 = 5, 
             b0 = 5)

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

str(output)
#> List of 8
#>  $ input                       :List of 15
#>   ..$ p_treatment          : num 0.08
#>   ..$ N_total              : num 900
#>   ..$ EndofStudy           : num 50
#>   ..$ prop_loss_to_followup: num 0.1
#>   ..$ interim_look         : num [1:3] 600 700 800
#>   ..$ lambda               : num [1:2] 0.3 1
#>   ..$ lambda_time          : num 25
#>   ..$ h0                   : num -0.03
#>   ..$ futility_prob        : num 0.05
#>   ..$ prob_ha              : num 0.95
#>   ..$ expected_success_prob: num 0.9
#>   ..$ alternative          : chr "less"
#>   ..$ N_impute             : num 25
#>   ..$ number_mcmc          : num 1000
#>   ..$ prior                : num [1:2] 5 5
#>  $ power                       :'data.frame':    4 obs. of  2 variables:
#>   ..$ interim_looks: num [1:4] 600 700 800 900
#>   ..$ power        : num [1:4] 0.4 0.7 0.8 0.9
#>  $ type1_error                 : num 0
#>  $ est_final                   : num [1:10] 0.0719 0.0976 0.0794 0.0922 0.0829 ...
#>  $ post_prob_accept_alternative: num [1:10] 1 0.881 0.999 0.958 0.985 1 0.992 0.968 0.987 0.993
#>  $ N_enrolled                  : num [1:10] 600 900 700 900 700 600 800 600 600 700
#>  $ stop_expect_success         : num [1:10] 1 0 1 0 1 1 1 1 1 1
#>  $ stop_futility               : num [1:10] 0 0 0 0 0 0 0 0 0 0

Historical data is not required to compute the simulation. However, if historical data is avaialble, it can be incorporated into the analysis using the discount prior approach as implemented in the bayesDP R packageg.

In our OPC trial, we’ll illustrate historical data incorporation. We’ll assume that the historical data had 5 events in 55 subjects. We’ll incorporate this historical data using the identity discount function.

For more details on the historical data incorporation method and computation, please see https://CRAN.R-project.org/package=bayesDP.

value <- value %>%
  historical_binomial(y0_treatment      = 5, 
                      N0_treatment      = 55,
                      discount_function = "identity",
                      y0_control        = NULL, 
                      N0_control        = NULL,
                      alpha_max         = 1, 
                      fix_alpha         = FALSE,
                      weibull_scale     = 0.135, 
                      weibull_shape     = 3)

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

str(output)
#> List of 8
#>  $ input                       :List of 22
#>   ..$ p_treatment          : num 0.08
#>   ..$ N_total              : num 900
#>   ..$ EndofStudy           : num 50
#>   ..$ prop_loss_to_followup: num 0.1
#>   ..$ interim_look         : num [1:3] 600 700 800
#>   ..$ lambda               : num [1:2] 0.3 1
#>   ..$ lambda_time          : num 25
#>   ..$ h0                   : num -0.03
#>   ..$ futility_prob        : num 0.05
#>   ..$ prob_ha              : num 0.95
#>   ..$ expected_success_prob: num 0.9
#>   ..$ alternative          : chr "less"
#>   ..$ N_impute             : num 25
#>   ..$ number_mcmc          : num 1000
#>   ..$ prior                : num [1:2] 5 5
#>   ..$ y0_treatment         : num 5
#>   ..$ N0_treatment         : num 55
#>   ..$ discount_function    : chr "identity"
#>   ..$ alpha_max            : num 1
#>   ..$ fix_alpha            : logi FALSE
#>   ..$ weibull_scale        : num 0.135
#>   ..$ weibull_shape        : num 3
#>  $ power                       :'data.frame':    4 obs. of  2 variables:
#>   ..$ interim_looks: num [1:4] 600 700 800 900
#>   ..$ power        : num [1:4] 0.5 0.6 0.6 0.7
#>  $ type1_error                 : num 0
#>  $ est_final                   : num [1:10] 0.087 0.0852 0.0773 0.0917 0.0747 ...
#>  $ post_prob_accept_alternative: num [1:10] 0.966 0.979 0.999 0.966 0.996 0.984 0.897 0.997 0.873 0.68
#>  $ N_enrolled                  : num [1:10] 700 600 600 900 600 600 600 600 900 900
#>  $ stop_expect_success         : num [1:10] 1 1 1 0 1 1 1 1 0 0
#>  $ 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 10 times to obtain results.

value <- binomial_outcome(p_treatment = 0.08) %>%
  enrollment_rate(lambda = c(0.3, 1), 
                  time   = 25) %>%
  study_details(total_sample_size     = 900, 
                study_period          = 50,
                interim_look          = c(600, 700, 800),
                prop_loss_to_followup = 0.10) %>%
  hypothesis(delta                 = -0.03, 
             futility_prob         = 0.05, 
             prob_accept_ha        = 0.95,
             expected_success_prob = 0.90, 
             alternative           = "less") %>%
  impute(no_of_impute = 25, 
         number_mcmc  = 1000) %>%
  beta_prior(a0 = 5, 
             b0 = 5) %>%
  historical_binomial(y0_treatment      = 5, 
                      N0_treatment      = 55,
                      discount_function = "identity",
                      y0_control        = NULL, 
                      N0_control        = NULL,
                      alpha_max         = 1, 
                      fix_alpha         = FALSE,
                      weibull_scale     = 0.135, 
                      weibull_shape     = 3) %>%
  simulate(no_of_sim = 10)

str(value)
#> List of 8
#>  $ input                       :List of 22
#>   ..$ p_treatment          : num 0.08
#>   ..$ lambda               : num [1:2] 0.3 1
#>   ..$ lambda_time          : num 25
#>   ..$ N_total              : num 900
#>   ..$ EndofStudy           : num 50
#>   ..$ interim_look         : num [1:3] 600 700 800
#>   ..$ prop_loss_to_followup: num 0.1
#>   ..$ h0                   : num -0.03
#>   ..$ futility_prob        : num 0.05
#>   ..$ prob_ha              : num 0.95
#>   ..$ expected_success_prob: num 0.9
#>   ..$ alternative          : chr "less"
#>   ..$ N_impute             : num 25
#>   ..$ number_mcmc          : num 1000
#>   ..$ prior                : num [1:2] 5 5
#>   ..$ y0_treatment         : num 5
#>   ..$ N0_treatment         : num 55
#>   ..$ discount_function    : chr "identity"
#>   ..$ alpha_max            : num 1
#>   ..$ fix_alpha            : logi FALSE
#>   ..$ weibull_scale        : num 0.135
#>   ..$ weibull_shape        : num 3
#>  $ power                       :'data.frame':    4 obs. of  2 variables:
#>   ..$ interim_looks: num [1:4] 600 700 800 900
#>   ..$ power        : num [1:4] 0.5 0.7 0.8 0.9
#>  $ type1_error                 : num 0
#>  $ est_final                   : num [1:10] 0.0999 0.0869 0.0873 0.0781 0.078 ...
#>  $ post_prob_accept_alternative: num [1:10] 0.845 0.991 0.968 0.997 0.996 0.994 0.991 1 0.983 0.994
#>  $ N_enrolled                  : num [1:10] 800 800 700 700 600 600 600 600 900 600
#>  $ stop_expect_success         : num [1:10] 0 1 1 1 1 1 1 1 0 1
#>  $ stop_futility               : num [1:10] 1 0 0 0 0 0 0 0 0 0

Two-arm trial

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

The binomial events are simulated using an event rate of 0.15 for the treatment group and 0.12 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. Further, we will illustrate block randomization. The following code simulates a trial 10 times using the piping procedure.


value <- binomial_outcome(p_treatment = 0.15, 
                          p_control   = 0.12) %>%
  study_details(total_sample_size     = 400, 
                study_period          = 30,
                interim_look          = 350,
                prop_loss_to_followup = 0.15) %>%
  hypothesis(delta                 = 0, 
             futility_prob         = 0.10, 
             prob_accept_ha        = 0.975,
             expected_success_prob = 1, 
             alternative           = "greater") %>%
  randomize(block_size          = 9, 
            randomization_ratio = c(4, 5)) %>%
  impute(no_of_impute = 20, 
         number_mcmc  = 5000) %>%
  beta_prior(a0 = 0, 
             b0 = 0) %>%
  simulate(no_of_sim = 10)

str(value)
#> List of 8
#>  $ input                       :List of 16
#>   ..$ p_treatment          : num 0.15
#>   ..$ p_control            : num 0.12
#>   ..$ N_total              : num 400
#>   ..$ EndofStudy           : num 30
#>   ..$ interim_look         : num 350
#>   ..$ prop_loss_to_followup: num 0.15
#>   ..$ h0                   : num 0
#>   ..$ futility_prob        : num 0.1
#>   ..$ prob_ha              : num 0.975
#>   ..$ expected_success_prob: num 1
#>   ..$ alternative          : chr "greater"
#>   ..$ block                : num 9
#>   ..$ rand_ratio           : num [1:2] 4 5
#>   ..$ N_impute             : num 20
#>   ..$ number_mcmc          : num 5000
#>   ..$ prior                : num [1:2] 0 0
#>  $ power                       :'data.frame':    2 obs. of  2 variables:
#>   ..$ interim_looks: num [1:2] 350 400
#>   ..$ power        : num [1:2] 0 0
#>  $ type1_error                 : num 0
#>  $ est_final                   : num [1:10] -0.019967 0.000161 0.031991 0.029932 0.056089 ...
#>  $ post_prob_accept_alternative: num [1:10] 0.316 0.504 0.834 0.779 0.942 ...
#>  $ N_enrolled                  : num [1:10] 350 350 350 350 400 400 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 0 0 1 1 1 1

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 binomialdata contains the results of 300 subjects from a two-arm trial with binomial outcome. The complete column indicates whether the outcome was observed, i.e., loss to follow-up.

data(binomialdata)

head(binomialdata)
#>         id treatment outcome complete
#> 1 Patient1         1       0        0
#> 2 Patient2         0       0        1
#> 3 Patient3         1       0        1
#> 4 Patient4         0       0        1
#> 5 Patient5         0       1        1
#> 6 Patient6         0       0        1

The minimum input needed to run an adaptive Bayesian trial is the data itself. The data_binomial 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_binomial(treatment = binomialdata$treatment,
                       outcome   = binomialdata$outcome,
                       complete  = binomialdata$complete)

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

str(out)
#> List of 11
#>  $ prob_of_accepting_alternative: num 0.95
#>  $ margin                       : num 0
#>  $ alternative                  : chr "greater"
#>  $ N_treatment                  : num 125
#>  $ N_control                    : int 137
#>  $ N_complete                   : num 262
#>  $ N_enrolled                   : int 300
#>  $ post_prob_accept_alternative : num 0.994
#>  $ est_final                    : num 0.124
#>  $ stop_futility                : num 0
#>  $ stop_expected_success        : num 1

We’ll now illustrate using piping to carry out the complete analysis. First, we’ll assume the following hypothesis: \[H_0:\theta_{treatment} - \theta_{control} <= 0.02 \quad H_A: \theta_{treatment} - \theta_{control} > 0.02\] The delta and alternative used to analyze the trial is 0.02 and “greater” 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_binomial(treatment = binomialdata$treatment,
                     outcome   = binomialdata$outcome,
                     complete  = binomialdata$complete) %>%
  hypothesis(delta                 = 0.02, 
             futility_prob         = 0.05, 
             prob_accept_ha        = 0.95,
             expected_success_prob = 0.90, 
             alternative           = "greater") %>%
  impute(no_of_impute = 50, 
         number_mcmc  = 10000) %>%
  beta_prior(a0 = 3, 
             b0 = 3) %>%
  historical_binomial(y0_treatment      = 12, 
                      N0_treatment      = 100,
                      y0_control        = NULL, 
                      N0_control        = NULL,
                      discount_function = "weibull",
                      alpha_max         = 1, 
                      fix_alpha         = FALSE,
                      weibull_scale     = 0.135, 
                      weibull_shape     = 3) %>%
  analysis(type = "binomial")

str(out)
#> List of 11
#>  $ prob_of_accepting_alternative: num 0.95
#>  $ margin                       : num 0.02
#>  $ alternative                  : chr "greater"
#>  $ N_treatment                  : num 125
#>  $ N_control                    : int 137
#>  $ N_complete                   : num 262
#>  $ N_enrolled                   : int 300
#>  $ post_prob_accept_alternative : num 0.956
#>  $ est_final                    : num 0.108
#>  $ stop_futility                : num 0
#>  $ stop_expected_success        : num 0