Nothing
R/binomial.R
In bayesCT: Simulation and Analysis of Adaptive Bayesian Clinical Trials
Defines functions
data_binomial
binomial_analysis
beta_prior
historical_binomial
binomial_outcome
binomialBACT
Documented in
beta_prior
binomial_analysis
binomialBACT
binomial_outcome
data_binomial
historical_binomial
#' @title Binomial counts for Bayesian Adaptive Trials
#'
#' @description Simulation for binomial counts for Bayesian Adaptive trial with
#' different inputs to control for power, sample size, type 1 error rate, etc.
#' @param p_treatment scalar. Proportion of events under the treatment arm.
#' @param p_control scalar. Proportion of events under the control arm.
#' @param y0_treatment scalar. Number of events for the historical treatment
#' arm.
#' @param N0_treatment scalar. Sample size of the historical treatment arm.
#' @param y0_control scalar. Number of events for the historical control arm.
#' @param N0_control scalar. Sample size of the historical control arm.
#' @param prior vector. Prior values of beta rate, Beta(a0, b0). The default is
#' set to Beta(1, 1).
#' @inheritParams normalBACT
#'
#' @return a list of output for a single trial simulation.
#' \describe{
#' \item{\code{p_treatment}}{
#' scalar. The input parameter of proportion of events in the
#' treatment group.}
#' \item{\code{p_control}}{
#' scalar. The input parameter of proportion of events in the
#' control group.}
#' \item{\code{prob_of_accepting_alternative}}{
#' scalar. The input parameter of probability threshold of accepting the
#' alternative.}
#' \item{\code{margin}}{
#' scalar. The margin input value of difference between mean estimate of treatment
#' and mean estimate of the control.}
#' \item{\code{alternative}}{
#' character. The input parameter of alternative hypothesis. }
#' \item{\code{interim_look}}{
#' vector. The sample size for each interim look.}
#' \item{\code{N_treatment}}{
#' scalar. The number of patients enrolled in the experimental group for
#' each simulation.}
#' \item{\code{N_control}}{
#' scalar. The number of patients enrolled in the control group for
#' each simulation.}
#' \item{\code{N_enrolled}}{
#' vector. The number of patients enrolled in the trial (sum of control
#' and experimental group for each simulation. )}
#' \item{\code{N_complete}}{
#' scalar. The number of patients who completed the trial and had no
#' loss to follow-up.}
#' \item{\code{post_prob_accept_alternative}}{
#' vector. The final probability of accepting the alternative
#' hypothesis after the analysis is done.}
#' \item{\code{est_final}}{
#' scalar. The final estimate of the difference in posterior estimate of
#' treatment and posterior estimate of the control group.}
#' \item{\code{stop_futility}}{
#' scalar. Did the trial stop for futility during imputation of patient
#' who had loss to follow up? 1 for yes and 0 for no.}
#' \item{\code{stop_expected_success}}{
#' scalar. Did the trial stop for early success during imputation of patient
#' who had loss to follow up? 1 for yes and 0 for no.}
#' \item{\code{est_interim}}{
#' scalar. The interim estimate of the difference in posterior estimate of
#' treatment and posterior estimate of the control group.}
#' }
#'
#'
#' @importFrom stats rbinom glm
#' @importFrom dplyr mutate filter group_by bind_rows select n
#' @importFrom bayesDP bdpbinomial
#' @export binomialBACT
binomialBACT <- function(
p_treatment,
p_control = NULL,
y0_treatment = NULL,
N0_treatment = NULL,
y0_control = NULL,
N0_control = NULL,
N_total,
lambda = 0.3,
lambda_time = NULL,
interim_look = NULL,
EndofStudy,
prior = c(1, 1),
block = 2, # block size for randomization
rand_ratio = c(1, 1), # randomization ratio in control to treatament (default 1:1)
prop_loss_to_followup = 0.10, # Proportion of loss in data
alternative = "greater", # the alternative hypothesis (either two-sided, greater, less)
h0 = 0, # Null hypothesis value
futility_prob = 0.05, # Futility probability
expected_success_prob = 0.9, # Expected success probability
prob_ha = 0.95, # Posterior probability of accepting alternative hypothesis
N_impute = 10, # Number of imputation simulations for predictive distribution
number_mcmc = 10000, # Number of posterior sampling
discount_function = "identity",
alpha_max = 1, # max weight on incorporating historical data
fix_alpha = FALSE, # fix alpha set weight of historical data to alpha_max
weibull_scale = 0.135, # weibull parameter
weibull_shape = 3, # weibull parameter
method = "fixed"
){
# checking p_control
if(!is.null(p_control)){stopifnot(p_control > 0 & p_control < 1)}
# checking historical data for treatment group
if(!is.null(y0_treatment) | !is.null(N0_treatment)){
stopifnot((y0_treatment > 0 & N0_treatment > 0), y0_treatment <= N0_treatment,
(discount_function == "identity" | discount_function == "weibull" |
discount_function == "scaledweibull"))
}
# checking historical data for control group
if(!is.null(y0_control) | !is.null(N0_control)){
stopifnot((y0_control > 0 & N0_control > 0), y0_control <= N0_control)
}
# checking interim_look
if(!is.null(interim_look)){
stopifnot(all(N_total > interim_look))
}
# checking other inputs
stopifnot((p_treatment < 1 & p_treatment > 0),
length(lambda) == (length(lambda_time) + 1),
EndofStudy > 0, block %% sum(rand_ratio) == 0,
(prop_loss_to_followup >= 0 & prop_loss_to_followup < 0.75),
(h0 >= -1 & h0 < 1), (futility_prob < 0.20 & futility_prob >= 0),
(expected_success_prob > 0.70 & expected_success_prob <= 1),
(prob_ha > 0.70 & prob_ha < 1), N_impute > 0)
#checking if alternative is right
if(alternative != "two-sided" & alternative != "greater" & alternative != "less"){
stop("The input for alternative is wrong!")
}
# assigining interim look and final look
analysis_at_enrollnumber <- c(interim_look, N_total)
# assignment of enrollment based on the enrollment function
enrollment <- enrollment(param = lambda, N_total = N_total, time = lambda_time)
# simulating group and treatment group assignment
if(!is.null(p_control)){
group <- randomization(N_total = N_total, block = block, allocation = rand_ratio)
}
else{
group <- rep(1, N_total)
}
# simulate binomial outcome
if(!is.null(p_control)){
sim <- rbinom(N_total, 1, prob = group * p_treatment + (1 - group) * p_control)
# dividing treatment and control
control <- sim[group == 0]
}
else{
sim <- rbinom(N_total, 1, p_treatment)
}
treatment <- sim[group == 1]
# simulate loss to follow-up
n_loss_to_fu <- ceiling(prop_loss_to_followup * N_total)
loss_to_fu <- rep(FALSE, N_total)
loss_to_fu[sample(1:N_total, n_loss_to_fu)] <- TRUE
# creating a new data.frame for all the variables
data_total <- data.frame(
Y = sim,
treatment = group,
enrollment = enrollment,
id = 1:N_total,
loss_to_fu = loss_to_fu)
# assigning stop_futility and expected success
stop_futility <- 0
stop_expected_success <- 0
if(length(analysis_at_enrollnumber) > 1){
for(i in 1:(length(analysis_at_enrollnumber) - 1)){
# Analysis at the `analysis_at_enrollnumber` look
# Indicators for subject type:
# - subject_enrolled: subject has data present in the current look
# - subject_impute_success: subject has data present in the current look but has not
# reached end of study or subject is loss to follow-up: needs BP change imputed
# - subject_impute_futility: subject has no data present in the current look;
# needs baseline BP and BP change imputed
data_interim <- data_total %>%
mutate(subject_enrolled = id <= analysis_at_enrollnumber[i],
subject_impute_futility = !subject_enrolled) %>%
group_by(subject_enrolled) %>%
mutate(subject_impute_success = (enrollment[analysis_at_enrollnumber[i]] - enrollment <= EndofStudy & subject_enrolled) |
(subject_enrolled & loss_to_fu))
# Carry out interim analysis on patients with complete data only
# - Set-up `new data` data frame
data <- data_interim %>%
filter(subject_enrolled,
!subject_impute_success)
# MLE of data at interim analysis
# MLE_int <- glm(Y ~ treatment, data = data, family = "binomial")
# assigning input for control arm given it is a single or double arm
if(!is.null(p_control)){
y_c <- sum(data$Y[data$treatment == 0])
N_c <- length(data$Y[data$treatment == 0])
}
else{
y_c <- NULL
N_c <- NULL
}
# Analyze data using discount funtion via binomial
post <- bdpbinomial(y_t = sum(data$Y[data$treatment == 1]),
N_t = length(data$Y[data$treatment == 1]),
y_c = y_c,
N_c = N_c,
y0_t = y0_treatment,
N0_t = N0_treatment,
y0_c = y0_control,
N0_c = N0_control,
discount_function = discount_function,
number_mcmc = number_mcmc,
a0 = prior[1],
b0 = prior[2],
alpha_max = alpha_max,
fix_alpha = fix_alpha,
weibull_scale = weibull_scale,
weibull_shape = weibull_shape,
method = method)
# Imputation phase futility and expected success - initialize counters
# for the current imputation phase
futility_test <- 0
expected_success_test <- 0
for(j in 1:N_impute){
#imputing the success for control group
data_control_success_impute <- data_interim %>%
filter(treatment == 0) %>%
mutate(Y_impute = ifelse(subject_impute_success & subject_enrolled,
rbinom(n(), 1, p_control),
Y))
# imputing success for treatment group
data_treatment_success_impute <- data_interim %>%
filter(treatment == 1) %>%
mutate(Y_impute = ifelse(subject_impute_success & subject_enrolled,
rbinom(n(), 1, p_treatment),
Y))
# combine the treatment and control imputed datasets
data_success_impute <- bind_rows(data_control_success_impute,
data_treatment_success_impute) %>%
mutate(Y = Y_impute) %>%
select(-Y_impute)
# create enrolled subject data frame for discount function analysis
data <- data_success_impute %>%
filter(subject_enrolled)
# assigning input for control arm given it is a single or double arm
if(!is.null(p_control)){
y_c <- sum(data$Y[data$treatment == 0])
N_c <- length(data$Y[data$treatment == 0])
}
else{
y_c <- NULL
N_c <- NULL
}
# analyze complete+imputed data using discount funtion via binomial
post_imp <- bdpbinomial(y_t = sum(data$Y[data$treatment == 1]),
N_t = length(data$Y[data$treatment == 1]),
y_c = y_c,
N_c = N_c,
y0_t = y0_treatment,
N0_t = N0_treatment,
y0_c = y0_control,
N0_c = N0_control,
discount_function = discount_function,
number_mcmc = number_mcmc,
a0 = prior[1],
b0 = prior[2],
alpha_max = alpha_max,
fix_alpha = fix_alpha,
weibull_scale = weibull_scale,
weibull_shape = weibull_shape,
method = method)
# estimation of the posterior effect for difference between test and
# control - If expected success, add 1 to the counter
if(!is.null(p_control)){
if(alternative == "two-sided"){
effect_imp <- post_imp$posterior_treatment$posterior - post_imp$posterior_control$posterior
success <- max(c(mean(effect_imp > h0), mean(-effect_imp > h0)))
}
else if(alternative == "greater"){
effect_imp <- post_imp$posterior_treatment$posterior - post_imp$posterior_control$posterior
success <- mean(effect_imp > h0)
}
else{
effect_imp <- post_imp$posterior_treatment$posterior - post_imp$posterior_control$posterior
success <- mean(-effect_imp > h0)
}
}
else{
effect_imp <- post_imp$final$posterior
if(alternative == "two-sided"){
success <- max(c(mean(effect_imp - p_treatment > h0), mean(p_treatment - effect_imp > h0)))
}
else if(alternative == "greater"){
success <- mean(effect_imp - p_treatment > h0)
}
else{
success <- mean(p_treatment - effect_imp > h0)
}
}
if(success > prob_ha){
expected_success_test <- expected_success_test + 1
}
##########################################################################
### Futility computations
##########################################################################
# For patients not enrolled, impute the outcome
data_control_futility_impute <- data_success_impute %>%
filter(treatment == 0) %>%
mutate(Y_impute = ifelse(subject_impute_futility,
rbinom(n(), 1, p_control),
Y))
data_treatment_futility_impute <- data_success_impute %>%
filter(treatment == 1) %>%
mutate(Y_impute = ifelse(subject_impute_futility,
rbinom(n(), 1, p_treatment),
Y))
# Combine the treatment and control imputed datasets
data_futility_impute <- bind_rows(data_control_futility_impute,
data_treatment_futility_impute) %>%
mutate(Y = Y_impute) %>%
select(-Y_impute)
# Create enrolled subject data frame for discount function analysis
data <- data_futility_impute
if(!is.null(p_control)){
y_c <- sum(data$Y[data$treatment == 0])
N_c <- length(data$Y[data$treatment == 0])
}
else{
y_c <- NULL
N_c <- NULL
}
# Analyze complete+imputed data using discount funtion via binomial
post_imp <- bdpbinomial(y_t = sum(data$Y[data$treatment == 1]),
N_t = length(data$Y[data$treatment == 1]),
y_c = y_c,
N_c = N_c,
y0_t = y0_treatment,
N0_t = N0_treatment,
y0_c = y0_control,
N0_c = N0_control,
discount_function = discount_function,
number_mcmc = number_mcmc,
a0 = prior[1],
b0 = prior[2],
alpha_max = alpha_max,
fix_alpha = fix_alpha,
weibull_scale = weibull_scale,
weibull_shape = weibull_shape,
method = method)
# Estimation of the posterior effect for difference between test and
# control
if(!is.null(p_control)){
if(alternative == "two-sided"){
effect_imp <- post_imp$posterior_treatment$posterior - post_imp$posterior_control$posterior
success <- max(c(mean(effect_imp > h0), mean(-effect_imp > h0)))
}
else if(alternative == "greater"){
effect_imp <- post_imp$posterior_treatment$posterior - post_imp$posterior_control$posterior
success <- mean(effect_imp > h0)
}
else{
effect_imp <- post_imp$posterior_treatment$posterior - post_imp$posterior_control$posterior
success <- mean(-effect_imp > h0)
}
}
else{
effect_imp <- post_imp$final$posterior
if(alternative == "two-sided"){
success <- max(c(mean(effect_imp - p_treatment > h0), mean(p_treatment - effect_imp > h0)))
}
else if(alternative == "greater"){
success <- mean(effect_imp - p_treatment > h0)
}
else{
success <- mean(p_treatment - effect_imp > h0)
}
}
# Increase futility counter by 1 if P(effect_imp < h0) > ha
if(success > prob_ha){
futility_test <- futility_test + 1
}
}
# print(analysis_at_enrollnumber[i])
# Test if expected success criteria met
if(expected_success_test / N_impute > expected_success_prob ){
stop_expected_success <- 1
stage_trial_stopped <- analysis_at_enrollnumber[i]
break
}
# Test if futility success criteria is met
if(futility_test / N_impute < futility_prob){
stop_futility <- 1
stage_trial_stopped <- analysis_at_enrollnumber[i]
break
}
# Stop study if at last interim look
if(analysis_at_enrollnumber[i + 1] == N_total){
stage_trial_stopped <- analysis_at_enrollnumber[i + 1]
break
}
}
##############################################################################
### Final analysis
##############################################################################
# Estimation of the posterior of the difference
if(!is.null(p_control)){
if(alternative == "two-sided"){
effect_int <- post$posterior_treatment$posterior - post$posterior_control$posterior
}
else if(alternative == "greater"){
effect_int <- post$posterior_treatment$posterior - post$posterior_control$posterior
}
else{
effect_int <- post$posterior_treatment$posterior - post$posterior_control$posterior
}
}
else{
effect_int <- post$final$posterior
}
# Number of patients enrolled at trial stop
N_enrolled <- nrow(data_interim[data_interim$id <= stage_trial_stopped, ])
}
# assigning stage trial stopped given no interim look
else{
N_enrolled <- N_total
stage_trial_stopped <- N_total
stop_futility <- 0
stop_expected_success <- 0
}
# All patients that have made it to the end of study
# - Subset out patients loss to follow-up
data_final <- data_total %>%
filter(id <= stage_trial_stopped,
!loss_to_fu)
# Compute the final MLE for the complete data using GLM function
# MLE <- glm(Y ~ treatment, data = data_final, family = "binomial")
if(!is.null(p_control)){
y_c <- sum(data_final$Y[data_final$treatment == 0])
N_c <- length(data_final$Y[data_final$treatment == 0])
}
else{
y_c <- NULL
N_c <- NULL
}
# Analyze complete data using discount funtion via binomial
post_final <- bdpbinomial(y_t = sum(data_final$Y[data_final$treatment == 1]),
N_t = length(data_final$Y[data_final$treatment == 1]),
y_c = y_c,
N_c = N_c,
y0_t = y0_treatment,
N0_t = N0_treatment,
y0_c = y0_control,
N0_c = N0_control,
number_mcmc = number_mcmc,
discount_function = discount_function,
a0 = prior[1],
b0 = prior[2],
alpha_max = alpha_max,
fix_alpha = fix_alpha,
weibull_scale = weibull_scale,
weibull_shape = weibull_shape,
method = method)
### Format and output results
# Posterior effect size: test vs control or treatment itself
if(!is.null(p_control)){
if(alternative == "two-sided"){
effect <- post_final$posterior_treatment$posterior - post_final$posterior_control$posterior
post_paa <- max(c(mean(effect > h0), mean(-effect > h0)))
}
else if(alternative == "greater"){
effect <- post_final$posterior_treatment$posterior - post_final$posterior_control$posterior
post_paa <- mean(effect > h0)
}
else{
effect <- post_final$posterior_treatment$posterior - post_final$posterior_control$posterior
post_paa <- mean(-effect > h0)
}
}
else{
effect <- post_final$final$posterior
if(alternative == "two-sided"){
post_paa <- max(c(mean(effect - p_treatment > h0), mean(p_treatment - effect > h0)))
}
else if(alternative == "greater"){
post_paa <- mean(effect - p_treatment > h0)
}
else{
post_paa <- mean(p_treatment - effect > h0)
}
}
N_treatment <- sum(data_final$treatment) # Total sample size analyzed - test group
N_control <- sum(!data_final$treatment) # Total sample size analyzed - control group
## output
results_list <- list(
p_treatment = p_treatment, # probability of treatment in binomial
p_control = p_control, # probability of control in binomial
prob_of_accepting_alternative = prob_ha,
margin = h0, # margin for error
alternative = alternative, # alternative hypothesis
interim_look = interim_look, # print interim looks
N_treatment = N_treatment,
N_control = N_control,
N_complete = N_treatment + N_control,
N_enrolled = N_enrolled, # Total sample size enrolled when trial stopped
N_max = N_total, # Total potential sample size
post_prob_accept_alternative = post_paa, # Posterior probability that alternative hypothesis is true
est_final = mean(effect), # Posterior Mean of treatment effect
stop_futility = stop_futility, # Did the trial stop for futility
stop_expected_success = stop_expected_success # Did the trial stop for expected success
#MLE_est = MLE$coe[2], # Treatment effect useing MLE
#MLE_est_interim = MLE_int$coe[2] # Treatment effect useing MLE at interim analysis
)
if(length(analysis_at_enrollnumber) > 1){
results_list <- c(results_list,
est_interim = mean(effect_int)) # Posterior Mean of treatment effect at interim analysis
}
#return results
results_list
}
## quiets concerns of R CMD check re: the .'s that appear in pipelines
if(getRversion() >= "2.15.1") utils::globalVariables(c("Y", "Y_impute", "id", "subject_enrolled",
"subject_impute_success", "subject_impute_futility"))
#' @title Proportion of an event in control and treatment
#'
#' @description Wrapper function for proportion of an event in control and treatment group with binomial outcome.
#'
#' @param p_treatment numeric. The proportion of an event in the treatment group, 0 < $p_treatment$ < 1.
#' @param p_control numeric. The proportion of an event in the control group, 0 < $p_control$ < 1.
#' @param .data NULL. stores the proportion of control and treatment, please do not fill it in.
#'
#' @return a list with proportion of control and treatment group.
#'
#' @examples binomial_outcome(p_control = 0.12, p_treatment = 0.08)
#' @export binomial_outcome
binomial_outcome <- function(p_treatment = NULL, p_control = NULL, .data = NULL){
.data$p_treatment <- p_treatment
.data$p_control <- p_control
.data
}
#' @title Historical data for binomial distribution
#'
#' @description Wrapper function for historical data from binomial outcome.
#'
#' @inheritParams normalBACT
#' @inheritParams binomialBACT
#' @param .data NULL. stores the proportion of control and treatment, please do not fill it in.
#'
#' @return a list with historical data for control and treatment group with the discount function.
#'
#' @examples
#' historical_binomial(y0_treatment = 5, N0_treatment = 10, y0_control = 15, N0_control = 23)
#' historical_binomial(y0_treatment = 5, N0_treatment = 10, y0_control = 15, N0_control = 23,
#' discount_function = "weibull", alpha_max = 1, fix_alpha = FALSE,
#' weibull_scale = 0.135, weibull_shape = 3)
#' @export historical_binomial
historical_binomial <- function(y0_treatment = NULL,
N0_treatment = NULL,
discount_function = "identity",
y0_control = NULL,
N0_control = NULL,
alpha_max = 1, # max weight on incorporating historical data
fix_alpha = FALSE, # fix alpha set weight of historical data to alpha_max
weibull_scale = 0.135, # weibull parameter
weibull_shape = 3, # weibull parameter
method = "fixed",
.data = NULL
){
.data$y0_treatment <- y0_treatment
.data$N0_treatment <- N0_treatment
.data$y0_control <- y0_control
.data$N0_control <- N0_control
.data$discount_function <- discount_function
.data$alpha_max <- alpha_max
.data$fix_alpha <- fix_alpha
.data$weibull_scale <- weibull_scale
.data$weibull_shape <- weibull_shape
.data$method <- method
.data
}
#' @title Beta prior for for control and treatment group
#'
#' @description Wrapper function for beta prior \code{beta(a0, b0)}.
#'
#' @param a0 numeric. The first shape parameter in the beta distribution
#' (\code{beta(a0, b0)}).
#' @param b0 numeric. The second shape parameter in the beta distribution
#' (\code{beta(a0, b0)}).
#' @param .data NULL. stores the beta prior rate, please do not fill it in.
#'
#' @return a list with vector of beta rate for the beta prior for treatment and control group.
#'
#' @examples beta_prior(a0 = 1, b0 = 1)
#' @export beta_prior
beta_prior <- function(a0 = 1, b0 = 1, .data = NULL){
.data$prior <- c(a0, b0)
.data
}
#' @title Analyzing Bayesian trial for binomial counts
#'
#' @description Function to analyze Bayesian trial for binomial count data
#' which allows early stopping and incorporation of historical data using
#' the discount function approach
#'
#' @param treatment vector. treatment assignment for patients, 1 for treatment group and
#' 0 for control group
#' @param outcome vector. binomial outcome of the trial, 1 for response (success or
#' failure), 0 for no response.
#' @param complete vector. similar length as treatment and outcome variable,
#' 1 for complete outcome, 0 for loss to follow up. If complete is not provided,
#' the dataset is assumed to be complete.
#' @inheritParams normalBACT
#' @inheritParams binomialBACT
#'
#' @importFrom stats rbinom glm
#' @importFrom dplyr mutate filter group_by bind_rows select n summarize
#' @importFrom bayesDP bdpbinomial
#'
#' @return a list of output for the Bayesian trial for binomial count.
#'
#' \describe{
#' \item{\code{prob_of_accepting_alternative}}{
#' scalar. The input parameter of probability of accepting the alternative.}
#' \item{\code{margin}}{
#' scalar. The margin input value of difference between mean estimate of treatment
#' and mean estimate of the control.}
#' \item{\code{alternative}}{
#' character. The input parameter of alternative hypothesis. }
#' \item{\code{N_treatment}}{
#' scalar. The number of patients enrolled in the experimental group for
#' each simulation.}
#' \item{\code{N_control}}{
#' scalar. The number of patients enrolled in the control group for
#' each simulation.}
#' \item{\code{N_enrolled}}{
#' vector. The number of patients enrolled in the trial (sum of control
#' and experimental group for each simulation. )}
#' \item{\code{N_complete}}{
#' scalar. The number of patients who completed the trial and had no
#' loss to follow-up.}
#' \item{\code{post_prob_accept_alternative}}{
#' vector. The final probability of accepting the alternative
#' hypothesis after the analysis is done.}
#' \item{\code{est_final}}{
#' scalar. The final estimate of the difference in posterior estimate of
#' treatment and posterior estimate of the control group.}
#' \item{\code{stop_futility}}{
#' scalar. Did the trial stop for futility during imputation of patient
#' who had loss to follow up? 1 for yes and 0 for no.}
#' \item{\code{stop_expected_success}}{
#' scalar. Did the trial stop for early success during imputation of patient
#' who had loss to follow up? 1 for yes and 0 for no.}
#'
#' }
#'
#' @export binomial_analysis
binomial_analysis <- function(
treatment,
outcome,
complete = NULL,
y0_treatment = NULL,
N0_treatment = NULL,
y0_control = NULL,
N0_control = NULL,
alternative = "greater",
N_impute = 10,
h0 = 0,
number_mcmc = 10000,
prob_ha = 0.95,
futility_prob = 0.10,
expected_success_prob = 0.90,
prior = c(1, 1),
discount_function = "identity",
fix_alpha = FALSE,
alpha_max = 1,
weibull_scale = 0.135,
weibull_shape = 3,
method = "fixed"
){
#if complete is NULL, assume the data is complete
if(is.null(complete)){
complete <- rep(1, length(outcome))
}
#reading the data
data_total <- data.frame(cbind(treatment, outcome, complete))
data_interim <- data_total %>%
mutate(futility = (complete == 0))
data <- data_interim %>%
filter(!futility)
prop <- data %>%
group_by(treatment) %>%
summarize(p_outcome = mean(outcome))
if(sum(data$treatment == 0) != 0){
y_c <- sum(data$outcome[data$treatment == 0])
N_c <- length(data$outcome[data$treatment == 0])
}
else{
y_c <- NULL
N_c <- NULL
}
# analyze the data using bayesDp
post <- bdpbinomial(y_t = sum(data$outcome[data$treatment == 1]),
N_t = length(data$outcome[data$treatment == 1]),
y_c = y_c,
N_c = N_c,
y0_t = y0_treatment,
N0_t = N0_treatment,
y0_c = y0_control,
N0_c = N0_control,
discount_function = discount_function,
number_mcmc = number_mcmc,
a0 = prior[1],
b0 = prior[2],
alpha_max = alpha_max,
fix_alpha = fix_alpha,
weibull_scale = weibull_scale,
weibull_shape = weibull_shape,
method = method)
# assigning stop_futility and expected success
stop_futility <- 0
stop_expected_success <- 0
expected_success_test <- 0
for(i in 1:N_impute){
data_control_success_impute <- data_interim %>%
filter(treatment == 0) %>%
mutate(outcome_impute = ifelse(futility,
rbinom(n(), 1, prop$p_outcome[1]),
outcome))
# imputing success for treatment group
data_treatment_success_impute <- data_interim %>%
filter(treatment == 1) %>%
mutate(outcome_impute = ifelse(futility,
rbinom(n(), 1, prop$p_outcome[2]),
outcome))
# combine the treatment and control imputed datasets
data_success_impute <- bind_rows(data_control_success_impute,
data_treatment_success_impute) %>%
mutate(outcome = outcome_impute) %>%
select(-outcome_impute)
# Create enrolled subject data frame for discount function analysis
data <- data_success_impute
# assigning input for control arm given it is a single or double arm
if(sum(data$treatment == 0) != 0){
y_c <- sum(data$outcome[data$treatment == 0])
N_c <- length(data$outcome[data$treatment == 0])
}
else{
y_c <- NULL
N_c <- NULL
}
# analyze complete+imputed data using discount funtion via binomial
post_imp <- bdpbinomial(y_t = sum(data$outcome[data$treatment == 1]),
N_t = length(data$outcome[data$treatment == 1]),
y_c = y_c,
N_c = N_c,
y0_t = y0_treatment,
N0_t = N0_treatment,
y0_c = y0_control,
N0_c = N0_control,
discount_function = discount_function,
number_mcmc = number_mcmc,
a0 = prior[1],
b0 = prior[2],
alpha_max = alpha_max,
fix_alpha = fix_alpha,
weibull_scale = weibull_scale,
weibull_shape = weibull_shape,
method = method)
if(sum(data$treatment == 0) != 0){
if(alternative == "two-sided"){
effect_imp <- post_imp$posterior_treatment$posterior - post_imp$posterior_control$posterior
success <- max(c(mean(effect_imp > h0), mean(-effect_imp > h0)))
}
else if(alternative == "greater"){
effect_imp <- post_imp$posterior_treatment$posterior - post_imp$posterior_control$posterior
success <- mean(effect_imp > h0)
}
else{
effect_imp <- post_imp$posterior_treatment$posterior - post_imp$posterior_control$posterior
success <- mean(-effect_imp > h0)
}
}
else{
effect_imp <- post_imp$final$posterior
if(alternative == "two-sided"){
success <- max(c(mean(effect_imp > h0), mean(effect_imp < h0)))
}
else if(alternative == "greater"){
success <- mean(effect_imp > h0)
}
else{
success <- mean(effect_imp < h0)
}
}
if(success > prob_ha){
expected_success_test <- expected_success_test + 1
}
}
if(expected_success_test / N_impute < futility_prob){
stop_futility <- 1
}
# Test if expected success criteria met
if(expected_success_test / N_impute > expected_success_prob ){
stop_expected_success <- 1
}
data_final <- data_interim %>%
filter(!futility)
if(sum(data_final$treatment == 0) != 0){
y_c <- sum(data_final$outcome[data_final$treatment == 0])
N_c <- length(data_final$outcome[data_final$treatment == 0])
}
else{
y_c <- NULL
N_c <- NULL
}
# Analyze complete data using discount funtion via binomial
post_final <- bdpbinomial(y_t = sum(data_final$outcome[data_final$treatment == 1]),
N_t = length(data_final$outcome[data_final$treatment == 1]),
y_c = y_c,
N_c = N_c,
y0_t = y0_treatment,
N0_t = N0_treatment,
y0_c = y0_control,
N0_c = N0_control,
number_mcmc = number_mcmc,
discount_function = discount_function,
a0 = prior[1],
b0 = prior[2],
alpha_max = alpha_max,
fix_alpha = fix_alpha,
weibull_scale = weibull_scale,
weibull_shape = weibull_shape,
method = method)
### Format and output results
# Posterior effect size: test vs control or treatment itself
if(sum(data_final$treatment == 0) != 0){
if(alternative == "two-sided"){
effect <- post_final$posterior_treatment$posterior - post_final$posterior_control$posterior
post_paa <- max(c(mean(effect > h0), mean(-effect > h0)))
}
else if(alternative == "greater"){
effect <- post_final$posterior_treatment$posterior - post_final$posterior_control$posterior
post_paa <- mean(effect > h0)
}
else{
effect <- post_final$posterior_treatment$posterior - post_final$posterior_control$posterior
post_paa <- mean(-effect > h0)
}
}
else{
effect <- post_final$final$posterior
if(alternative == "two-sided"){
post_paa <- max(c(mean(effect > h0), mean(effect < h0)))
}
else if(alternative == "greater"){
post_paa <- mean(effect > h0)
}
else{
post_paa <- mean(effect < h0)
}
}
N_treatment <- sum(data_final$treatment) # Total sample size analyzed - test group
N_control <- sum(!data_final$treatment) # Total sample size analyzed - control group
N_enrolled <- dim(data_total)[1]
#estimating prop
prop <- data %>%
group_by(treatment) %>%
summarize(p_outcome = mean(outcome))
## output
results_list <- list(
prob_of_accepting_alternative = prob_ha,
margin = h0, # margin for error
alternative = alternative, # alternative hypothesis
N_treatment = N_treatment,
N_control = N_control,
N_complete = N_treatment + N_control,
N_enrolled = N_enrolled, # Total sample size enrolled when trial stopped
post_prob_accept_alternative = post_paa, # Posterior probability that alternative hypothesis is true
est_final = mean(effect), # Posterior Mean of treatment effect
stop_futility = stop_futility, # Did the trial stop for futility
stop_expected_success = stop_expected_success # Did the trial stop for expected success
#MLE_est = MLE$coe[2], # Treatment effect useing MLE
#MLE_est_interim = MLE_int$coe[2] # Treatment effect useing MLE at interim analysis
)
return(results_list)
}
## quiets concerns of R CMD check re: the .'s that appear in pipelines
if(getRversion() >= "2.15.1") utils::globalVariables(c("complete", "outcome", "outcome_impute", "id",
"futility", "treatment",
"subject_impute_success",
"subject_impute_futility", "p_outcome"))
#' @title Data file for binomial analysis
#'
#' @description Wrapper function for data file in binomial analysis.
#'
#' @inheritParams binomial_analysis
#' @param .data NULL. stores the binomial data for analysis, please do not fill it in.
#'
#' @return a list with treatment, outcome and loss to follow up vector with binomial
#' outcome.
#'
#' @examples data_binomial(treatment = c(0, 1), outcome = c(1, 1), complete = c(1, 1))
#' @export data_binomial
data_binomial <- function(treatment, outcome, complete, .data = NULL){
.data$treatment <- treatment
.data$outcome <- outcome
.data$complete <- complete
.data
}
Try the bayesCT package in your browser
Any scripts or data that you put into this service are public.
bayesCT documentation built on July 2, 2020, 2:34 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions data_binomial binomial_analysis beta_prior historical_binomial binomial_outcome binomialBACT
Documented in beta_prior binomial_analysis binomialBACT binomial_outcome data_binomial historical_binomial
#' @title Binomial counts for Bayesian Adaptive Trials
#'
#' @description Simulation for binomial counts for Bayesian Adaptive trial with
#' different inputs to control for power, sample size, type 1 error rate, etc.
#' @param p_treatment scalar. Proportion of events under the treatment arm.
#' @param p_control scalar. Proportion of events under the control arm.
#' @param y0_treatment scalar. Number of events for the historical treatment
#' arm.
#' @param N0_treatment scalar. Sample size of the historical treatment arm.
#' @param y0_control scalar. Number of events for the historical control arm.
#' @param N0_control scalar. Sample size of the historical control arm.
#' @param prior vector. Prior values of beta rate, Beta(a0, b0). The default is
#' set to Beta(1, 1).
#' @inheritParams normalBACT
#'
#' @return a list of output for a single trial simulation.
#' \describe{
#' \item{\code{p_treatment}}{
#' scalar. The input parameter of proportion of events in the
#' treatment group.}
#' \item{\code{p_control}}{
#' scalar. The input parameter of proportion of events in the
#' control group.}
#' \item{\code{prob_of_accepting_alternative}}{
#' scalar. The input parameter of probability threshold of accepting the
#' alternative.}
#' \item{\code{margin}}{
#' scalar. The margin input value of difference between mean estimate of treatment
#' and mean estimate of the control.}
#' \item{\code{alternative}}{
#' character. The input parameter of alternative hypothesis. }
#' \item{\code{interim_look}}{
#' vector. The sample size for each interim look.}
#' \item{\code{N_treatment}}{
#' scalar. The number of patients enrolled in the experimental group for
#' each simulation.}
#' \item{\code{N_control}}{
#' scalar. The number of patients enrolled in the control group for
#' each simulation.}
#' \item{\code{N_enrolled}}{
#' vector. The number of patients enrolled in the trial (sum of control
#' and experimental group for each simulation. )}
#' \item{\code{N_complete}}{
#' scalar. The number of patients who completed the trial and had no
#' loss to follow-up.}
#' \item{\code{post_prob_accept_alternative}}{
#' vector. The final probability of accepting the alternative
#' hypothesis after the analysis is done.}
#' \item{\code{est_final}}{
#' scalar. The final estimate of the difference in posterior estimate of
#' treatment and posterior estimate of the control group.}
#' \item{\code{stop_futility}}{
#' scalar. Did the trial stop for futility during imputation of patient
#' who had loss to follow up? 1 for yes and 0 for no.}
#' \item{\code{stop_expected_success}}{
#' scalar. Did the trial stop for early success during imputation of patient
#' who had loss to follow up? 1 for yes and 0 for no.}
#' \item{\code{est_interim}}{
#' scalar. The interim estimate of the difference in posterior estimate of
#' treatment and posterior estimate of the control group.}
#' }
#'
#'
#' @importFrom stats rbinom glm
#' @importFrom dplyr mutate filter group_by bind_rows select n
#' @importFrom bayesDP bdpbinomial
#' @export binomialBACT
binomialBACT <- function(
p_treatment,
p_control = NULL,
y0_treatment = NULL,
N0_treatment = NULL,
y0_control = NULL,
N0_control = NULL,
N_total,
lambda = 0.3,
lambda_time = NULL,
interim_look = NULL,
EndofStudy,
prior = c(1, 1),
block = 2, # block size for randomization
rand_ratio = c(1, 1), # randomization ratio in control to treatament (default 1:1)
prop_loss_to_followup = 0.10, # Proportion of loss in data
alternative = "greater", # the alternative hypothesis (either two-sided, greater, less)
h0 = 0, # Null hypothesis value
futility_prob = 0.05, # Futility probability
expected_success_prob = 0.9, # Expected success probability
prob_ha = 0.95, # Posterior probability of accepting alternative hypothesis
N_impute = 10, # Number of imputation simulations for predictive distribution
number_mcmc = 10000, # Number of posterior sampling
discount_function = "identity",
alpha_max = 1, # max weight on incorporating historical data
fix_alpha = FALSE, # fix alpha set weight of historical data to alpha_max
weibull_scale = 0.135, # weibull parameter
weibull_shape = 3, # weibull parameter
method = "fixed"
){
# checking p_control
if(!is.null(p_control)){stopifnot(p_control > 0 & p_control < 1)}
# checking historical data for treatment group
if(!is.null(y0_treatment) | !is.null(N0_treatment)){
stopifnot((y0_treatment > 0 & N0_treatment > 0), y0_treatment <= N0_treatment,
(discount_function == "identity" | discount_function == "weibull" |
discount_function == "scaledweibull"))
}
# checking historical data for control group
if(!is.null(y0_control) | !is.null(N0_control)){
stopifnot((y0_control > 0 & N0_control > 0), y0_control <= N0_control)
}
# checking interim_look
if(!is.null(interim_look)){
stopifnot(all(N_total > interim_look))
}
# checking other inputs
stopifnot((p_treatment < 1 & p_treatment > 0),
length(lambda) == (length(lambda_time) + 1),
EndofStudy > 0, block %% sum(rand_ratio) == 0,
(prop_loss_to_followup >= 0 & prop_loss_to_followup < 0.75),
(h0 >= -1 & h0 < 1), (futility_prob < 0.20 & futility_prob >= 0),
(expected_success_prob > 0.70 & expected_success_prob <= 1),
(prob_ha > 0.70 & prob_ha < 1), N_impute > 0)
#checking if alternative is right
if(alternative != "two-sided" & alternative != "greater" & alternative != "less"){
stop("The input for alternative is wrong!")
}
# assigining interim look and final look
analysis_at_enrollnumber <- c(interim_look, N_total)
# assignment of enrollment based on the enrollment function
enrollment <- enrollment(param = lambda, N_total = N_total, time = lambda_time)
# simulating group and treatment group assignment
if(!is.null(p_control)){
group <- randomization(N_total = N_total, block = block, allocation = rand_ratio)
}
else{
group <- rep(1, N_total)
}
# simulate binomial outcome
if(!is.null(p_control)){
sim <- rbinom(N_total, 1, prob = group * p_treatment + (1 - group) * p_control)
# dividing treatment and control
control <- sim[group == 0]
}
else{
sim <- rbinom(N_total, 1, p_treatment)
}
treatment <- sim[group == 1]
# simulate loss to follow-up
n_loss_to_fu <- ceiling(prop_loss_to_followup * N_total)
loss_to_fu <- rep(FALSE, N_total)
loss_to_fu[sample(1:N_total, n_loss_to_fu)] <- TRUE
# creating a new data.frame for all the variables
data_total <- data.frame(
Y = sim,
treatment = group,
enrollment = enrollment,
id = 1:N_total,
loss_to_fu = loss_to_fu)
# assigning stop_futility and expected success
stop_futility <- 0
stop_expected_success <- 0
if(length(analysis_at_enrollnumber) > 1){
for(i in 1:(length(analysis_at_enrollnumber) - 1)){
# Analysis at the `analysis_at_enrollnumber` look
# Indicators for subject type:
# - subject_enrolled: subject has data present in the current look
# - subject_impute_success: subject has data present in the current look but has not
# reached end of study or subject is loss to follow-up: needs BP change imputed
# - subject_impute_futility: subject has no data present in the current look;
# needs baseline BP and BP change imputed
data_interim <- data_total %>%
mutate(subject_enrolled = id <= analysis_at_enrollnumber[i],
subject_impute_futility = !subject_enrolled) %>%
group_by(subject_enrolled) %>%
mutate(subject_impute_success = (enrollment[analysis_at_enrollnumber[i]] - enrollment <= EndofStudy & subject_enrolled) |
(subject_enrolled & loss_to_fu))
# Carry out interim analysis on patients with complete data only
# - Set-up `new data` data frame
data <- data_interim %>%
filter(subject_enrolled,
!subject_impute_success)
# MLE of data at interim analysis
# MLE_int <- glm(Y ~ treatment, data = data, family = "binomial")
# assigning input for control arm given it is a single or double arm
if(!is.null(p_control)){
y_c <- sum(data$Y[data$treatment == 0])
N_c <- length(data$Y[data$treatment == 0])
}
else{
y_c <- NULL
N_c <- NULL
}
# Analyze data using discount funtion via binomial
post <- bdpbinomial(y_t = sum(data$Y[data$treatment == 1]),
N_t = length(data$Y[data$treatment == 1]),
y_c = y_c,
N_c = N_c,
y0_t = y0_treatment,
N0_t = N0_treatment,
y0_c = y0_control,
N0_c = N0_control,
discount_function = discount_function,
number_mcmc = number_mcmc,
a0 = prior[1],
b0 = prior[2],
alpha_max = alpha_max,
fix_alpha = fix_alpha,
weibull_scale = weibull_scale,
weibull_shape = weibull_shape,
method = method)
# Imputation phase futility and expected success - initialize counters
# for the current imputation phase
futility_test <- 0
expected_success_test <- 0
for(j in 1:N_impute){
#imputing the success for control group
data_control_success_impute <- data_interim %>%
filter(treatment == 0) %>%
mutate(Y_impute = ifelse(subject_impute_success & subject_enrolled,
rbinom(n(), 1, p_control),
Y))
# imputing success for treatment group
data_treatment_success_impute <- data_interim %>%
filter(treatment == 1) %>%
mutate(Y_impute = ifelse(subject_impute_success & subject_enrolled,
rbinom(n(), 1, p_treatment),
Y))
# combine the treatment and control imputed datasets
data_success_impute <- bind_rows(data_control_success_impute,
data_treatment_success_impute) %>%
mutate(Y = Y_impute) %>%
select(-Y_impute)
# create enrolled subject data frame for discount function analysis
data <- data_success_impute %>%
filter(subject_enrolled)
# assigning input for control arm given it is a single or double arm
if(!is.null(p_control)){
y_c <- sum(data$Y[data$treatment == 0])
N_c <- length(data$Y[data$treatment == 0])
}
else{
y_c <- NULL
N_c <- NULL
}
# analyze complete+imputed data using discount funtion via binomial
post_imp <- bdpbinomial(y_t = sum(data$Y[data$treatment == 1]),
N_t = length(data$Y[data$treatment == 1]),
y_c = y_c,
N_c = N_c,
y0_t = y0_treatment,
N0_t = N0_treatment,
y0_c = y0_control,
N0_c = N0_control,
discount_function = discount_function,
number_mcmc = number_mcmc,
a0 = prior[1],
b0 = prior[2],
alpha_max = alpha_max,
fix_alpha = fix_alpha,
weibull_scale = weibull_scale,
weibull_shape = weibull_shape,
method = method)
# estimation of the posterior effect for difference between test and
# control - If expected success, add 1 to the counter
if(!is.null(p_control)){
if(alternative == "two-sided"){
effect_imp <- post_imp$posterior_treatment$posterior - post_imp$posterior_control$posterior
success <- max(c(mean(effect_imp > h0), mean(-effect_imp > h0)))
}
else if(alternative == "greater"){
effect_imp <- post_imp$posterior_treatment$posterior - post_imp$posterior_control$posterior
success <- mean(effect_imp > h0)
}
else{
effect_imp <- post_imp$posterior_treatment$posterior - post_imp$posterior_control$posterior
success <- mean(-effect_imp > h0)
}
}
else{
effect_imp <- post_imp$final$posterior
if(alternative == "two-sided"){
success <- max(c(mean(effect_imp - p_treatment > h0), mean(p_treatment - effect_imp > h0)))
}
else if(alternative == "greater"){
success <- mean(effect_imp - p_treatment > h0)
}
else{
success <- mean(p_treatment - effect_imp > h0)
}
}
if(success > prob_ha){
expected_success_test <- expected_success_test + 1
}
##########################################################################
### Futility computations
##########################################################################
# For patients not enrolled, impute the outcome
data_control_futility_impute <- data_success_impute %>%
filter(treatment == 0) %>%
mutate(Y_impute = ifelse(subject_impute_futility,
rbinom(n(), 1, p_control),
Y))
data_treatment_futility_impute <- data_success_impute %>%
filter(treatment == 1) %>%
mutate(Y_impute = ifelse(subject_impute_futility,
rbinom(n(), 1, p_treatment),
Y))
# Combine the treatment and control imputed datasets
data_futility_impute <- bind_rows(data_control_futility_impute,
data_treatment_futility_impute) %>%
mutate(Y = Y_impute) %>%
select(-Y_impute)
# Create enrolled subject data frame for discount function analysis
data <- data_futility_impute
if(!is.null(p_control)){
y_c <- sum(data$Y[data$treatment == 0])
N_c <- length(data$Y[data$treatment == 0])
}
else{
y_c <- NULL
N_c <- NULL
}
# Analyze complete+imputed data using discount funtion via binomial
post_imp <- bdpbinomial(y_t = sum(data$Y[data$treatment == 1]),
N_t = length(data$Y[data$treatment == 1]),
y_c = y_c,
N_c = N_c,
y0_t = y0_treatment,
N0_t = N0_treatment,
y0_c = y0_control,
N0_c = N0_control,
discount_function = discount_function,
number_mcmc = number_mcmc,
a0 = prior[1],
b0 = prior[2],
alpha_max = alpha_max,
fix_alpha = fix_alpha,
weibull_scale = weibull_scale,
weibull_shape = weibull_shape,
method = method)
# Estimation of the posterior effect for difference between test and
# control
if(!is.null(p_control)){
if(alternative == "two-sided"){
effect_imp <- post_imp$posterior_treatment$posterior - post_imp$posterior_control$posterior
success <- max(c(mean(effect_imp > h0), mean(-effect_imp > h0)))
}
else if(alternative == "greater"){
effect_imp <- post_imp$posterior_treatment$posterior - post_imp$posterior_control$posterior
success <- mean(effect_imp > h0)
}
else{
effect_imp <- post_imp$posterior_treatment$posterior - post_imp$posterior_control$posterior
success <- mean(-effect_imp > h0)
}
}
else{
effect_imp <- post_imp$final$posterior
if(alternative == "two-sided"){
success <- max(c(mean(effect_imp - p_treatment > h0), mean(p_treatment - effect_imp > h0)))
}
else if(alternative == "greater"){
success <- mean(effect_imp - p_treatment > h0)
}
else{
success <- mean(p_treatment - effect_imp > h0)
}
}
# Increase futility counter by 1 if P(effect_imp < h0) > ha
if(success > prob_ha){
futility_test <- futility_test + 1
}
}
# print(analysis_at_enrollnumber[i])
# Test if expected success criteria met
if(expected_success_test / N_impute > expected_success_prob ){
stop_expected_success <- 1
stage_trial_stopped <- analysis_at_enrollnumber[i]
break
}
# Test if futility success criteria is met
if(futility_test / N_impute < futility_prob){
stop_futility <- 1
stage_trial_stopped <- analysis_at_enrollnumber[i]
break
}
# Stop study if at last interim look
if(analysis_at_enrollnumber[i + 1] == N_total){
stage_trial_stopped <- analysis_at_enrollnumber[i + 1]
break
}
}
##############################################################################
### Final analysis
##############################################################################
# Estimation of the posterior of the difference
if(!is.null(p_control)){
if(alternative == "two-sided"){
effect_int <- post$posterior_treatment$posterior - post$posterior_control$posterior
}
else if(alternative == "greater"){
effect_int <- post$posterior_treatment$posterior - post$posterior_control$posterior
}
else{
effect_int <- post$posterior_treatment$posterior - post$posterior_control$posterior
}
}
else{
effect_int <- post$final$posterior
}
# Number of patients enrolled at trial stop
N_enrolled <- nrow(data_interim[data_interim$id <= stage_trial_stopped, ])
}
# assigning stage trial stopped given no interim look
else{
N_enrolled <- N_total
stage_trial_stopped <- N_total
stop_futility <- 0
stop_expected_success <- 0
}
# All patients that have made it to the end of study
# - Subset out patients loss to follow-up
data_final <- data_total %>%
filter(id <= stage_trial_stopped,
!loss_to_fu)
# Compute the final MLE for the complete data using GLM function
# MLE <- glm(Y ~ treatment, data = data_final, family = "binomial")
if(!is.null(p_control)){
y_c <- sum(data_final$Y[data_final$treatment == 0])
N_c <- length(data_final$Y[data_final$treatment == 0])
}
else{
y_c <- NULL
N_c <- NULL
}
# Analyze complete data using discount funtion via binomial
post_final <- bdpbinomial(y_t = sum(data_final$Y[data_final$treatment == 1]),
N_t = length(data_final$Y[data_final$treatment == 1]),
y_c = y_c,
N_c = N_c,
y0_t = y0_treatment,
N0_t = N0_treatment,
y0_c = y0_control,
N0_c = N0_control,
number_mcmc = number_mcmc,
discount_function = discount_function,
a0 = prior[1],
b0 = prior[2],
alpha_max = alpha_max,
fix_alpha = fix_alpha,
weibull_scale = weibull_scale,
weibull_shape = weibull_shape,
method = method)
### Format and output results
# Posterior effect size: test vs control or treatment itself
if(!is.null(p_control)){
if(alternative == "two-sided"){
effect <- post_final$posterior_treatment$posterior - post_final$posterior_control$posterior
post_paa <- max(c(mean(effect > h0), mean(-effect > h0)))
}
else if(alternative == "greater"){
effect <- post_final$posterior_treatment$posterior - post_final$posterior_control$posterior
post_paa <- mean(effect > h0)
}
else{
effect <- post_final$posterior_treatment$posterior - post_final$posterior_control$posterior
post_paa <- mean(-effect > h0)
}
}
else{
effect <- post_final$final$posterior
if(alternative == "two-sided"){
post_paa <- max(c(mean(effect - p_treatment > h0), mean(p_treatment - effect > h0)))
}
else if(alternative == "greater"){
post_paa <- mean(effect - p_treatment > h0)
}
else{
post_paa <- mean(p_treatment - effect > h0)
}
}
N_treatment <- sum(data_final$treatment) # Total sample size analyzed - test group
N_control <- sum(!data_final$treatment) # Total sample size analyzed - control group
## output
results_list <- list(
p_treatment = p_treatment, # probability of treatment in binomial
p_control = p_control, # probability of control in binomial
prob_of_accepting_alternative = prob_ha,
margin = h0, # margin for error
alternative = alternative, # alternative hypothesis
interim_look = interim_look, # print interim looks
N_treatment = N_treatment,
N_control = N_control,
N_complete = N_treatment + N_control,
N_enrolled = N_enrolled, # Total sample size enrolled when trial stopped
N_max = N_total, # Total potential sample size
post_prob_accept_alternative = post_paa, # Posterior probability that alternative hypothesis is true
est_final = mean(effect), # Posterior Mean of treatment effect
stop_futility = stop_futility, # Did the trial stop for futility
stop_expected_success = stop_expected_success # Did the trial stop for expected success
#MLE_est = MLE$coe[2], # Treatment effect useing MLE
#MLE_est_interim = MLE_int$coe[2] # Treatment effect useing MLE at interim analysis
)
if(length(analysis_at_enrollnumber) > 1){
results_list <- c(results_list,
est_interim = mean(effect_int)) # Posterior Mean of treatment effect at interim analysis
}
#return results
results_list
}
## quiets concerns of R CMD check re: the .'s that appear in pipelines
if(getRversion() >= "2.15.1") utils::globalVariables(c("Y", "Y_impute", "id", "subject_enrolled",
"subject_impute_success", "subject_impute_futility"))
#' @title Proportion of an event in control and treatment
#'
#' @description Wrapper function for proportion of an event in control and treatment group with binomial outcome.
#'
#' @param p_treatment numeric. The proportion of an event in the treatment group, 0 < $p_treatment$ < 1.
#' @param p_control numeric. The proportion of an event in the control group, 0 < $p_control$ < 1.
#' @param .data NULL. stores the proportion of control and treatment, please do not fill it in.
#'
#' @return a list with proportion of control and treatment group.
#'
#' @examples binomial_outcome(p_control = 0.12, p_treatment = 0.08)
#' @export binomial_outcome
binomial_outcome <- function(p_treatment = NULL, p_control = NULL, .data = NULL){
.data$p_treatment <- p_treatment
.data$p_control <- p_control
.data
}
#' @title Historical data for binomial distribution
#'
#' @description Wrapper function for historical data from binomial outcome.
#'
#' @inheritParams normalBACT
#' @inheritParams binomialBACT
#' @param .data NULL. stores the proportion of control and treatment, please do not fill it in.
#'
#' @return a list with historical data for control and treatment group with the discount function.
#'
#' @examples
#' historical_binomial(y0_treatment = 5, N0_treatment = 10, y0_control = 15, N0_control = 23)
#' historical_binomial(y0_treatment = 5, N0_treatment = 10, y0_control = 15, N0_control = 23,
#' discount_function = "weibull", alpha_max = 1, fix_alpha = FALSE,
#' weibull_scale = 0.135, weibull_shape = 3)
#' @export historical_binomial
historical_binomial <- function(y0_treatment = NULL,
N0_treatment = NULL,
discount_function = "identity",
y0_control = NULL,
N0_control = NULL,
alpha_max = 1, # max weight on incorporating historical data
fix_alpha = FALSE, # fix alpha set weight of historical data to alpha_max
weibull_scale = 0.135, # weibull parameter
weibull_shape = 3, # weibull parameter
method = "fixed",
.data = NULL
){
.data$y0_treatment <- y0_treatment
.data$N0_treatment <- N0_treatment
.data$y0_control <- y0_control
.data$N0_control <- N0_control
.data$discount_function <- discount_function
.data$alpha_max <- alpha_max
.data$fix_alpha <- fix_alpha
.data$weibull_scale <- weibull_scale
.data$weibull_shape <- weibull_shape
.data$method <- method
.data
}
#' @title Beta prior for for control and treatment group
#'
#' @description Wrapper function for beta prior \code{beta(a0, b0)}.
#'
#' @param a0 numeric. The first shape parameter in the beta distribution
#' (\code{beta(a0, b0)}).
#' @param b0 numeric. The second shape parameter in the beta distribution
#' (\code{beta(a0, b0)}).
#' @param .data NULL. stores the beta prior rate, please do not fill it in.
#'
#' @return a list with vector of beta rate for the beta prior for treatment and control group.
#'
#' @examples beta_prior(a0 = 1, b0 = 1)
#' @export beta_prior
beta_prior <- function(a0 = 1, b0 = 1, .data = NULL){
.data$prior <- c(a0, b0)
.data
}
#' @title Analyzing Bayesian trial for binomial counts
#'
#' @description Function to analyze Bayesian trial for binomial count data
#' which allows early stopping and incorporation of historical data using
#' the discount function approach
#'
#' @param treatment vector. treatment assignment for patients, 1 for treatment group and
#' 0 for control group
#' @param outcome vector. binomial outcome of the trial, 1 for response (success or
#' failure), 0 for no response.
#' @param complete vector. similar length as treatment and outcome variable,
#' 1 for complete outcome, 0 for loss to follow up. If complete is not provided,
#' the dataset is assumed to be complete.
#' @inheritParams normalBACT
#' @inheritParams binomialBACT
#'
#' @importFrom stats rbinom glm
#' @importFrom dplyr mutate filter group_by bind_rows select n summarize
#' @importFrom bayesDP bdpbinomial
#'
#' @return a list of output for the Bayesian trial for binomial count.
#'
#' \describe{
#' \item{\code{prob_of_accepting_alternative}}{
#' scalar. The input parameter of probability of accepting the alternative.}
#' \item{\code{margin}}{
#' scalar. The margin input value of difference between mean estimate of treatment
#' and mean estimate of the control.}
#' \item{\code{alternative}}{
#' character. The input parameter of alternative hypothesis. }
#' \item{\code{N_treatment}}{
#' scalar. The number of patients enrolled in the experimental group for
#' each simulation.}
#' \item{\code{N_control}}{
#' scalar. The number of patients enrolled in the control group for
#' each simulation.}
#' \item{\code{N_enrolled}}{
#' vector. The number of patients enrolled in the trial (sum of control
#' and experimental group for each simulation. )}
#' \item{\code{N_complete}}{
#' scalar. The number of patients who completed the trial and had no
#' loss to follow-up.}
#' \item{\code{post_prob_accept_alternative}}{
#' vector. The final probability of accepting the alternative
#' hypothesis after the analysis is done.}
#' \item{\code{est_final}}{
#' scalar. The final estimate of the difference in posterior estimate of
#' treatment and posterior estimate of the control group.}
#' \item{\code{stop_futility}}{
#' scalar. Did the trial stop for futility during imputation of patient
#' who had loss to follow up? 1 for yes and 0 for no.}
#' \item{\code{stop_expected_success}}{
#' scalar. Did the trial stop for early success during imputation of patient
#' who had loss to follow up? 1 for yes and 0 for no.}
#'
#' }
#'
#' @export binomial_analysis
binomial_analysis <- function(
treatment,
outcome,
complete = NULL,
y0_treatment = NULL,
N0_treatment = NULL,
y0_control = NULL,
N0_control = NULL,
alternative = "greater",
N_impute = 10,
h0 = 0,
number_mcmc = 10000,
prob_ha = 0.95,
futility_prob = 0.10,
expected_success_prob = 0.90,
prior = c(1, 1),
discount_function = "identity",
fix_alpha = FALSE,
alpha_max = 1,
weibull_scale = 0.135,
weibull_shape = 3,
method = "fixed"
){
#if complete is NULL, assume the data is complete
if(is.null(complete)){
complete <- rep(1, length(outcome))
}
#reading the data
data_total <- data.frame(cbind(treatment, outcome, complete))
data_interim <- data_total %>%
mutate(futility = (complete == 0))
data <- data_interim %>%
filter(!futility)
prop <- data %>%
group_by(treatment) %>%
summarize(p_outcome = mean(outcome))
if(sum(data$treatment == 0) != 0){
y_c <- sum(data$outcome[data$treatment == 0])
N_c <- length(data$outcome[data$treatment == 0])
}
else{
y_c <- NULL
N_c <- NULL
}
# analyze the data using bayesDp
post <- bdpbinomial(y_t = sum(data$outcome[data$treatment == 1]),
N_t = length(data$outcome[data$treatment == 1]),
y_c = y_c,
N_c = N_c,
y0_t = y0_treatment,
N0_t = N0_treatment,
y0_c = y0_control,
N0_c = N0_control,
discount_function = discount_function,
number_mcmc = number_mcmc,
a0 = prior[1],
b0 = prior[2],
alpha_max = alpha_max,
fix_alpha = fix_alpha,
weibull_scale = weibull_scale,
weibull_shape = weibull_shape,
method = method)
# assigning stop_futility and expected success
stop_futility <- 0
stop_expected_success <- 0
expected_success_test <- 0
for(i in 1:N_impute){
data_control_success_impute <- data_interim %>%
filter(treatment == 0) %>%
mutate(outcome_impute = ifelse(futility,
rbinom(n(), 1, prop$p_outcome[1]),
outcome))
# imputing success for treatment group
data_treatment_success_impute <- data_interim %>%
filter(treatment == 1) %>%
mutate(outcome_impute = ifelse(futility,
rbinom(n(), 1, prop$p_outcome[2]),
outcome))
# combine the treatment and control imputed datasets
data_success_impute <- bind_rows(data_control_success_impute,
data_treatment_success_impute) %>%
mutate(outcome = outcome_impute) %>%
select(-outcome_impute)
# Create enrolled subject data frame for discount function analysis
data <- data_success_impute
# assigning input for control arm given it is a single or double arm
if(sum(data$treatment == 0) != 0){
y_c <- sum(data$outcome[data$treatment == 0])
N_c <- length(data$outcome[data$treatment == 0])
}
else{
y_c <- NULL
N_c <- NULL
}
# analyze complete+imputed data using discount funtion via binomial
post_imp <- bdpbinomial(y_t = sum(data$outcome[data$treatment == 1]),
N_t = length(data$outcome[data$treatment == 1]),
y_c = y_c,
N_c = N_c,
y0_t = y0_treatment,
N0_t = N0_treatment,
y0_c = y0_control,
N0_c = N0_control,
discount_function = discount_function,
number_mcmc = number_mcmc,
a0 = prior[1],
b0 = prior[2],
alpha_max = alpha_max,
fix_alpha = fix_alpha,
weibull_scale = weibull_scale,
weibull_shape = weibull_shape,
method = method)
if(sum(data$treatment == 0) != 0){
if(alternative == "two-sided"){
effect_imp <- post_imp$posterior_treatment$posterior - post_imp$posterior_control$posterior
success <- max(c(mean(effect_imp > h0), mean(-effect_imp > h0)))
}
else if(alternative == "greater"){
effect_imp <- post_imp$posterior_treatment$posterior - post_imp$posterior_control$posterior
success <- mean(effect_imp > h0)
}
else{
effect_imp <- post_imp$posterior_treatment$posterior - post_imp$posterior_control$posterior
success <- mean(-effect_imp > h0)
}
}
else{
effect_imp <- post_imp$final$posterior
if(alternative == "two-sided"){
success <- max(c(mean(effect_imp > h0), mean(effect_imp < h0)))
}
else if(alternative == "greater"){
success <- mean(effect_imp > h0)
}
else{
success <- mean(effect_imp < h0)
}
}
if(success > prob_ha){
expected_success_test <- expected_success_test + 1
}
}
if(expected_success_test / N_impute < futility_prob){
stop_futility <- 1
}
# Test if expected success criteria met
if(expected_success_test / N_impute > expected_success_prob ){
stop_expected_success <- 1
}
data_final <- data_interim %>%
filter(!futility)
if(sum(data_final$treatment == 0) != 0){
y_c <- sum(data_final$outcome[data_final$treatment == 0])
N_c <- length(data_final$outcome[data_final$treatment == 0])
}
else{
y_c <- NULL
N_c <- NULL
}
# Analyze complete data using discount funtion via binomial
post_final <- bdpbinomial(y_t = sum(data_final$outcome[data_final$treatment == 1]),
N_t = length(data_final$outcome[data_final$treatment == 1]),
y_c = y_c,
N_c = N_c,
y0_t = y0_treatment,
N0_t = N0_treatment,
y0_c = y0_control,
N0_c = N0_control,
number_mcmc = number_mcmc,
discount_function = discount_function,
a0 = prior[1],
b0 = prior[2],
alpha_max = alpha_max,
fix_alpha = fix_alpha,
weibull_scale = weibull_scale,
weibull_shape = weibull_shape,
method = method)
### Format and output results
# Posterior effect size: test vs control or treatment itself
if(sum(data_final$treatment == 0) != 0){
if(alternative == "two-sided"){
effect <- post_final$posterior_treatment$posterior - post_final$posterior_control$posterior
post_paa <- max(c(mean(effect > h0), mean(-effect > h0)))
}
else if(alternative == "greater"){
effect <- post_final$posterior_treatment$posterior - post_final$posterior_control$posterior
post_paa <- mean(effect > h0)
}
else{
effect <- post_final$posterior_treatment$posterior - post_final$posterior_control$posterior
post_paa <- mean(-effect > h0)
}
}
else{
effect <- post_final$final$posterior
if(alternative == "two-sided"){
post_paa <- max(c(mean(effect > h0), mean(effect < h0)))
}
else if(alternative == "greater"){
post_paa <- mean(effect > h0)
}
else{
post_paa <- mean(effect < h0)
}
}
N_treatment <- sum(data_final$treatment) # Total sample size analyzed - test group
N_control <- sum(!data_final$treatment) # Total sample size analyzed - control group
N_enrolled <- dim(data_total)[1]
#estimating prop
prop <- data %>%
group_by(treatment) %>%
summarize(p_outcome = mean(outcome))
## output
results_list <- list(
prob_of_accepting_alternative = prob_ha,
margin = h0, # margin for error
alternative = alternative, # alternative hypothesis
N_treatment = N_treatment,
N_control = N_control,
N_complete = N_treatment + N_control,
N_enrolled = N_enrolled, # Total sample size enrolled when trial stopped
post_prob_accept_alternative = post_paa, # Posterior probability that alternative hypothesis is true
est_final = mean(effect), # Posterior Mean of treatment effect
stop_futility = stop_futility, # Did the trial stop for futility
stop_expected_success = stop_expected_success # Did the trial stop for expected success
#MLE_est = MLE$coe[2], # Treatment effect useing MLE
#MLE_est_interim = MLE_int$coe[2] # Treatment effect useing MLE at interim analysis
)
return(results_list)
}
## quiets concerns of R CMD check re: the .'s that appear in pipelines
if(getRversion() >= "2.15.1") utils::globalVariables(c("complete", "outcome", "outcome_impute", "id",
"futility", "treatment",
"subject_impute_success",
"subject_impute_futility", "p_outcome"))
#' @title Data file for binomial analysis
#'
#' @description Wrapper function for data file in binomial analysis.
#'
#' @inheritParams binomial_analysis
#' @param .data NULL. stores the binomial data for analysis, please do not fill it in.
#'
#' @return a list with treatment, outcome and loss to follow up vector with binomial
#' outcome.
#'
#' @examples data_binomial(treatment = c(0, 1), outcome = c(1, 1), complete = c(1, 1))
#' @export data_binomial
data_binomial <- function(treatment, outcome, complete, .data = NULL){
.data$treatment <- treatment
.data$outcome <- outcome
.data$complete <- complete
.data
}
Try the bayesCT package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.