R/binomial.R
In thevaachandereng/BACT: 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 trials 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 treatment (default 1:1)
prop_loss_to_followup = 0.10, # Proportion of loss in data
alternative = "greater", # 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!")
}
# Checking if N_impute is <= number_mcmc
if (N_impute > number_mcmc) {
stop("The number of imputations must not be greater than the number of MCMC draws!")
}
# Assigning 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 imputation
# - subject_impute_futility: subject has no data present in the current look;
# needs imputation
data_interim <- data_total %>%
mutate(
subject_enrolled = id <= analysis_at_enrollnumber[i],
subject_impute_futility = !subject_enrolled,
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 function 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
# Sub-sample values from posterior to use for imputation stage
id_impute <- sample(1:number_mcmc, N_impute)
p_treatment_imp <- post$posterior_treatment$posterior[id_impute]
p_control_imp <- post$posterior_control$posterior[id_impute]
for (j in 1:N_impute) {
##########################################################################
### Expected success computations
##########################################################################
# 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_imp[j]),
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_imp[j]),
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 function 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_imp[j]),
Y))
data_treatment_futility_impute <- data_success_impute %>%
filter(treatment == 1) %>%
mutate(Y_impute = ifelse(subject_impute_futility,
rbinom(n(), 1, p_treatment_imp[j]),
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 function via bdpbinomial
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
}
}
# 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, ])
} else {
# Assigning stage trial stopped given no interim look
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 function via bdpbinomial
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 using MLE
#MLE_est_interim = MLE_int$coe[2] # Treatment effect using 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
return(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. Should not
#' be edited by the user.
#'
#' @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. Should not
#' be edited by the user.
#'
#' @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 proportion of control and treatment. Should not
#' be edited by the user.
#'
#' @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 a Bayesian trial for binomial counts
#'
#' @description Function to analyze a Bayesian trial for binomial count data
#' which allows for 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.
#' @param N_max_treatment integer. Maximum allowable sample size for the
#' treatment arm (including the currently enrolled subjects). Default is NULL,
#' meaning we are already at the final analysis.
#' @param N_max_control integer. Maximum allowable sample size for the control
#' arm (including the currently enrolled subjects). Default is NULL, meaning
#' we are already at the final analysis.
#' @inheritParams normalBACT
#' @inheritParams binomialBACT
#'
#' @details If the enrollment size is at the final sample size, i.e. the maximum
#' allowable sample size for the trial, then this function is not of practical
#' use since there is no opportunity to stop trial enrollment. In such a case,
#' it is expected that the follow-up will be conducted per the study protocol
#' and a final analysis made.
#'
#' @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}}{
#' integer. 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}}{
#' integer. The number of patients enrolled in the experimental group for
#' each simulation.}
#' \item{\code{N_control}}{
#' integer. 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}}{
#' integer. The number of patients who completed the trial and had no
#' loss to follow-up.}
#' \item{\code{N_max_treatment}}{
#' integer. Maximum allowable sample size for the treatment arm
#' (including the currently enrolled subjects).}
#' \item{\code{N_max_control}}{
#' integer. Maximum allowable sample size for the control arm
#' (including the currently enrolled subjects).}
#' \item{\code{post_prob_accept_alternative}}{
#' scalar. 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}}{
#' integer. 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}}{
#' integer. 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,
p_treatment = 0,
N_max_treatment = NULL,
N_max_control = 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))
}
# Final analysis or interim?
if (is.null(N_max_treatment) & is.null(N_max_control)) {
final_analysis <- TRUE
N_max_treatment <- sum(treatment == 1)
N_max_control <- sum(treatment == 0)
} else if (length(complete) < (N_max_treatment + N_max_control)) {
final_analysis <- FALSE
# Number of subjects still to enroll
N_horizon <- N_max_treatment + N_max_control - length(complete)
} else if (length(complete) == (N_max_treatment + N_max_control)) {
final_analysis <- TRUE
} else if (length(complete) > (N_max_treatment + N_max_control)) {
stop("Number of enrolled subjects exceeds maximum defined in analysis plan!")
}
if (final_analysis) {
message("Are you at the final sample size? There is no opportunity to stop enrollment.\n
Consider setting N_max_treatment and/or N_max_control")
}
# Reading the data
data_interim <- data.frame(cbind(treatment, outcome, complete))
data_interim$subject_enrolled <- TRUE
# Add additional subjects who will be enrolled under current design:
# = N_max_treatment + N_max_control - number of subjects currently enrolled
if (!final_analysis) {
data_horizon <- data.frame(
"treatment" = c(rep(1, N_max_treatment), rep(0, N_max_control)),
"outcome" = rep(NA, N_horizon),
"complete" = rep(0, N_horizon),
"subject_enrolled" = rep(FALSE, N_horizon)
)
data_interim <- rbind(data_interim, data_horizon)
}
# 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 imputation
# - subject_impute_futility: subject has no data present in the current look;
# needs imputation
data_interim <- data_interim %>%
mutate(subject_impute_futility = !subject_enrolled,
subject_impute_success = (complete == 0))
data <- data_interim %>%
filter(subject_enrolled,
!subject_impute_success)
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 data using discount function via bdpbinomial
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 stop_expected_success
stop_futility <- 0
stop_expected_success <- 0
futility_test <- 0
expected_success_test <- 0
# Sub-sample values from posterior to use for imputation stage
id_impute <- sample(1:number_mcmc, N_impute)
p_treatment_imp <- post$posterior_treatment$posterior[id_impute]
p_control_imp <- post$posterior_control$posterior[id_impute]
for(i in 1:N_impute) {
##########################################################################
### Expected success computations
##########################################################################
# Imputing success for control group
data_control_success_impute <- data_interim %>%
filter(treatment == 0) %>%
mutate(outcome_impute = ifelse(subject_impute_success & subject_enrolled,
rbinom(n(), 1, p_control_imp[i]),
outcome))
# Imputing success for treatment group
data_treatment_success_impute <- data_interim %>%
filter(treatment == 1) %>%
mutate(outcome_impute = ifelse(subject_impute_success & subject_enrolled,
rbinom(n(), 1, p_treatment_imp[i]),
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 %>%
filter(subject_enrolled)
# 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 function via bdpbinomial
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)
# Estimation of the posterior effect for difference between test and
# control. If expected success, add 1 to the counter.
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
}
##########################################################################
### Futility computations
##########################################################################
# For patients not enrolled, impute the outcome
data_control_futility_impute <- data_success_impute %>%
filter(treatment == 0) %>%
mutate(outcome_impute = ifelse(subject_impute_futility,
rbinom(n(), 1, p_control_imp[i]),
outcome))
data_treatment_futility_impute <- data_success_impute %>%
filter(treatment == 1) %>%
mutate(outcome_impute = ifelse(subject_impute_futility,
rbinom(n(), 1, p_treatment_imp[i]),
outcome))
# Combine the treatment and control imputed datasets
data_futility_impute <- bind_rows(data_control_futility_impute,
data_treatment_futility_impute) %>%
mutate(outcome = outcome_impute) %>%
select(-outcome_impute)
# Create enrolled subject data frame for discount function analysis
data <- data_futility_impute
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 function via bdpbinomial
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)
# Estimation of the posterior effect for difference between test and control
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)
}
}
# Increase futility counter by 1 if P(effect_imp < h0) > ha
if (success > prob_ha) {
futility_test <- futility_test + 1
}
}
# Test if expected success criteria met
if (expected_success_test / N_impute > expected_success_prob) {
stop_expected_success <- 1
}
# Test if futility success criteria is met
if (futility_test / N_impute < futility_prob) {
stop_futility <- 1
}
##############################################################################
### Summary at the interim analysis
##############################################################################
data <- data_interim %>%
filter(subject_enrolled,
!subject_impute_success)
# Posterior effect size: test vs. control or treatment itself
if (sum(data$treatment == 0) != 0) {
if (alternative == "two-sided") {
effect <- post$posterior_treatment$posterior - post$posterior_control$posterior
post_paa <- max(c(mean(effect > h0), mean(-effect > h0)))
} else if (alternative == "greater") {
effect <- post$posterior_treatment$posterior - post$posterior_control$posterior
post_paa <- mean(effect > h0)
} else {
effect <- post$posterior_treatment$posterior - post$posterior_control$posterior
post_paa <- mean(-effect > h0)
}
} else {
effect <- post$final$posterior
if (alternative == "two-sided") {
success <- max(c(mean(effect > h0), mean(effect < h0)))
} else if (alternative == "greater") {
success <- mean(effect > h0)
} else {
success <- mean(effect < h0)
}
}
N_treatment <- sum(data$treatment) # Total sample size analyzed / evaluable - test group
N_control <- sum(!data$treatment) # Total sample size analyzed / evaluable - control group
N_enrolled <- sum(data_interim$subject_enrolled)
# 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
N_max_treatment = N_max_treatment,
N_max_control = N_max_control,
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 using MLE
#MLE_est_interim = MLE_int$coe[2] # Treatment effect using 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",
"futility", "treatment",
"subject_impute_success",
"subject_impute_futility"))
#' @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. Should not be
#' edited by user.
#'
#' @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
}
thevaachandereng/BACT documentation built on July 24, 2020, 2:35 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 trials 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 treatment (default 1:1)
prop_loss_to_followup = 0.10, # Proportion of loss in data
alternative = "greater", # 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!")
}
# Checking if N_impute is <= number_mcmc
if (N_impute > number_mcmc) {
stop("The number of imputations must not be greater than the number of MCMC draws!")
}
# Assigning 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 imputation
# - subject_impute_futility: subject has no data present in the current look;
# needs imputation
data_interim <- data_total %>%
mutate(
subject_enrolled = id <= analysis_at_enrollnumber[i],
subject_impute_futility = !subject_enrolled,
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 function 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
# Sub-sample values from posterior to use for imputation stage
id_impute <- sample(1:number_mcmc, N_impute)
p_treatment_imp <- post$posterior_treatment$posterior[id_impute]
p_control_imp <- post$posterior_control$posterior[id_impute]
for (j in 1:N_impute) {
##########################################################################
### Expected success computations
##########################################################################
# 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_imp[j]),
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_imp[j]),
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 function 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_imp[j]),
Y))
data_treatment_futility_impute <- data_success_impute %>%
filter(treatment == 1) %>%
mutate(Y_impute = ifelse(subject_impute_futility,
rbinom(n(), 1, p_treatment_imp[j]),
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 function via bdpbinomial
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
}
}
# 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, ])
} else {
# Assigning stage trial stopped given no interim look
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 function via bdpbinomial
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 using MLE
#MLE_est_interim = MLE_int$coe[2] # Treatment effect using 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
return(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. Should not
#' be edited by the user.
#'
#' @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. Should not
#' be edited by the user.
#'
#' @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 proportion of control and treatment. Should not
#' be edited by the user.
#'
#' @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 a Bayesian trial for binomial counts
#'
#' @description Function to analyze a Bayesian trial for binomial count data
#' which allows for 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.
#' @param N_max_treatment integer. Maximum allowable sample size for the
#' treatment arm (including the currently enrolled subjects). Default is NULL,
#' meaning we are already at the final analysis.
#' @param N_max_control integer. Maximum allowable sample size for the control
#' arm (including the currently enrolled subjects). Default is NULL, meaning
#' we are already at the final analysis.
#' @inheritParams normalBACT
#' @inheritParams binomialBACT
#'
#' @details If the enrollment size is at the final sample size, i.e. the maximum
#' allowable sample size for the trial, then this function is not of practical
#' use since there is no opportunity to stop trial enrollment. In such a case,
#' it is expected that the follow-up will be conducted per the study protocol
#' and a final analysis made.
#'
#' @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}}{
#' integer. 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}}{
#' integer. The number of patients enrolled in the experimental group for
#' each simulation.}
#' \item{\code{N_control}}{
#' integer. 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}}{
#' integer. The number of patients who completed the trial and had no
#' loss to follow-up.}
#' \item{\code{N_max_treatment}}{
#' integer. Maximum allowable sample size for the treatment arm
#' (including the currently enrolled subjects).}
#' \item{\code{N_max_control}}{
#' integer. Maximum allowable sample size for the control arm
#' (including the currently enrolled subjects).}
#' \item{\code{post_prob_accept_alternative}}{
#' scalar. 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}}{
#' integer. 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}}{
#' integer. 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,
p_treatment = 0,
N_max_treatment = NULL,
N_max_control = 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))
}
# Final analysis or interim?
if (is.null(N_max_treatment) & is.null(N_max_control)) {
final_analysis <- TRUE
N_max_treatment <- sum(treatment == 1)
N_max_control <- sum(treatment == 0)
} else if (length(complete) < (N_max_treatment + N_max_control)) {
final_analysis <- FALSE
# Number of subjects still to enroll
N_horizon <- N_max_treatment + N_max_control - length(complete)
} else if (length(complete) == (N_max_treatment + N_max_control)) {
final_analysis <- TRUE
} else if (length(complete) > (N_max_treatment + N_max_control)) {
stop("Number of enrolled subjects exceeds maximum defined in analysis plan!")
}
if (final_analysis) {
message("Are you at the final sample size? There is no opportunity to stop enrollment.\n
Consider setting N_max_treatment and/or N_max_control")
}
# Reading the data
data_interim <- data.frame(cbind(treatment, outcome, complete))
data_interim$subject_enrolled <- TRUE
# Add additional subjects who will be enrolled under current design:
# = N_max_treatment + N_max_control - number of subjects currently enrolled
if (!final_analysis) {
data_horizon <- data.frame(
"treatment" = c(rep(1, N_max_treatment), rep(0, N_max_control)),
"outcome" = rep(NA, N_horizon),
"complete" = rep(0, N_horizon),
"subject_enrolled" = rep(FALSE, N_horizon)
)
data_interim <- rbind(data_interim, data_horizon)
}
# 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 imputation
# - subject_impute_futility: subject has no data present in the current look;
# needs imputation
data_interim <- data_interim %>%
mutate(subject_impute_futility = !subject_enrolled,
subject_impute_success = (complete == 0))
data <- data_interim %>%
filter(subject_enrolled,
!subject_impute_success)
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 data using discount function via bdpbinomial
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 stop_expected_success
stop_futility <- 0
stop_expected_success <- 0
futility_test <- 0
expected_success_test <- 0
# Sub-sample values from posterior to use for imputation stage
id_impute <- sample(1:number_mcmc, N_impute)
p_treatment_imp <- post$posterior_treatment$posterior[id_impute]
p_control_imp <- post$posterior_control$posterior[id_impute]
for(i in 1:N_impute) {
##########################################################################
### Expected success computations
##########################################################################
# Imputing success for control group
data_control_success_impute <- data_interim %>%
filter(treatment == 0) %>%
mutate(outcome_impute = ifelse(subject_impute_success & subject_enrolled,
rbinom(n(), 1, p_control_imp[i]),
outcome))
# Imputing success for treatment group
data_treatment_success_impute <- data_interim %>%
filter(treatment == 1) %>%
mutate(outcome_impute = ifelse(subject_impute_success & subject_enrolled,
rbinom(n(), 1, p_treatment_imp[i]),
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 %>%
filter(subject_enrolled)
# 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 function via bdpbinomial
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)
# Estimation of the posterior effect for difference between test and
# control. If expected success, add 1 to the counter.
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
}
##########################################################################
### Futility computations
##########################################################################
# For patients not enrolled, impute the outcome
data_control_futility_impute <- data_success_impute %>%
filter(treatment == 0) %>%
mutate(outcome_impute = ifelse(subject_impute_futility,
rbinom(n(), 1, p_control_imp[i]),
outcome))
data_treatment_futility_impute <- data_success_impute %>%
filter(treatment == 1) %>%
mutate(outcome_impute = ifelse(subject_impute_futility,
rbinom(n(), 1, p_treatment_imp[i]),
outcome))
# Combine the treatment and control imputed datasets
data_futility_impute <- bind_rows(data_control_futility_impute,
data_treatment_futility_impute) %>%
mutate(outcome = outcome_impute) %>%
select(-outcome_impute)
# Create enrolled subject data frame for discount function analysis
data <- data_futility_impute
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 function via bdpbinomial
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)
# Estimation of the posterior effect for difference between test and control
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)
}
}
# Increase futility counter by 1 if P(effect_imp < h0) > ha
if (success > prob_ha) {
futility_test <- futility_test + 1
}
}
# Test if expected success criteria met
if (expected_success_test / N_impute > expected_success_prob) {
stop_expected_success <- 1
}
# Test if futility success criteria is met
if (futility_test / N_impute < futility_prob) {
stop_futility <- 1
}
##############################################################################
### Summary at the interim analysis
##############################################################################
data <- data_interim %>%
filter(subject_enrolled,
!subject_impute_success)
# Posterior effect size: test vs. control or treatment itself
if (sum(data$treatment == 0) != 0) {
if (alternative == "two-sided") {
effect <- post$posterior_treatment$posterior - post$posterior_control$posterior
post_paa <- max(c(mean(effect > h0), mean(-effect > h0)))
} else if (alternative == "greater") {
effect <- post$posterior_treatment$posterior - post$posterior_control$posterior
post_paa <- mean(effect > h0)
} else {
effect <- post$posterior_treatment$posterior - post$posterior_control$posterior
post_paa <- mean(-effect > h0)
}
} else {
effect <- post$final$posterior
if (alternative == "two-sided") {
success <- max(c(mean(effect > h0), mean(effect < h0)))
} else if (alternative == "greater") {
success <- mean(effect > h0)
} else {
success <- mean(effect < h0)
}
}
N_treatment <- sum(data$treatment) # Total sample size analyzed / evaluable - test group
N_control <- sum(!data$treatment) # Total sample size analyzed / evaluable - control group
N_enrolled <- sum(data_interim$subject_enrolled)
# 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
N_max_treatment = N_max_treatment,
N_max_control = N_max_control,
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 using MLE
#MLE_est_interim = MLE_int$coe[2] # Treatment effect using 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",
"futility", "treatment",
"subject_impute_success",
"subject_impute_futility"))
#' @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. Should not be
#' edited by user.
#'
#' @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
}
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.