R/survival.R
In thevaachandereng/BACT: Simulation and Analysis of Adaptive Bayesian Clinical Trials
Defines functions
data_survival
survival_outcome
gamma_prior
historical_survival
survival_analysis
survivalBACT
Documented in
data_survival
gamma_prior
historical_survival
survival_analysis
survivalBACT
survival_outcome
#' @title Time-to-event outcome for Bayesian Adaptive trials
#'
#' @description Simulation for time-to-event outcome for Bayesian Adaptive trial with
#' different inputs to control for power, sample size, type 1 error rate, etc.
#'
#' @param hazard_treatment vector. Constant hazard rates under the treatment arm.
#' @param hazard_control vector. Constant hazard rates under the control arm.
#' @inheritParams pw_exp_sim
#' @inheritParams normalBACT
#' @inheritParams survival_analysis
#'
#' @return A list of output for a single trial simulation.
#' \describe{
#' \item{\code{lambda_treatment}}{
#' vector. The input parameter of constant hazard rates in the
#' treatment group.}
#' \item{\code{cutpoint_treatment}}{
#' vector. The change-point vector when the constant hazard rate(s) changes for
#' the treatment group.}
#' \item{\code{lambda_control}}{
#' vector. The input parameter of constant hazard rates in the
#' control group.}
#' \item{\code{cutpoint_control}}{
#' vector. The change-point vector when the constant hazard rate(s) changes for
#' 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 treatment group for
#' each simulation.}
#' \item{\code{event_treatment}}{
#' scalar. The number of events in the treatment group for
#' each simulation.}
#' \item{\code{N_control}}{
#' scalar. The number of patients enrolled in the control group for
#' each simulation.}
#' \item{\code{event_control}}{
#' scalar. The number of events 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 runif
#' @importFrom dplyr mutate filter group_by bind_cols bind_rows select ungroup
#' @importFrom bayesDP bdpsurvival
#'
#' @export survivalBACT
survivalBACT <- function(
hazard_treatment,
cutpoint = NULL,
hazard_control = NULL,
N_total,
breaks = NULL,
time0 = NULL,
treatment0 = NULL,
event0 = NULL,
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", # 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 interim_look
if (!is.null(interim_look)) {
stopifnot(all(N_total > interim_look))
}
# Combining the histrotical data
if (!is.null(time0)) {
data0 <- data.frame(cbind(time = time0, event = event0, treatment = treatment0))
} else {
data0 <- NULL
}
# Checking if alternative is right
if (alternative != "two-sided" & alternative != "greater" & alternative != "less") {
stop("The input for alternative is wrong!")
}
# If cutpoint is not NULL, assign breaks to cutpoint
if (!is.null(cutpoint)) {
breaks <- cutpoint
}
# Make sure breaks is not more than EndofStudy
if (!is.null(breaks)) {
stopifnot(any(breaks < EndofStudy))
}
if (is.null(hazard_control) & h0 == 0) {
h0 <- 0.5
}
# 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(hazard_control)) {
group <- randomization(N_total = N_total, block = block, allocation = rand_ratio)
} else {
group <- rep(1, N_total)
}
time <- rep(NA, length = N_total)
event <- rep(NA, length = N_total)
# Simulate TTE outcome
if (!is.null(hazard_control)) {
sim_control <- pw_exp_sim(hazard = hazard_control,
n = sum(!group),
maxtime = EndofStudy,
cutpoint = cutpoint)
time[which(!group)] <- sim_control$time
event[which(!group)] <- sim_control$event
}
sim_treatment <- pw_exp_sim(hazard = hazard_treatment,
n = sum(group),
maxtime = EndofStudy,
cutpoint = cutpoint)
time[which(!!group)] <- sim_treatment$time
event[which(!!group)] <- sim_treatment$event
if (is.null(breaks)) {
breaks <- mean(time)
}
# 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(
time = time,
treatment = group,
event = event,
enrollment = enrollment,
id = 1:N_total,
loss_to_fu = loss_to_fu)
# Subjects lost are uniformly distributed
data_total$time[data_total$loss_to_fu] <- runif(n_loss_to_fu,
0, data_total$time[data_total$loss_to_fu])
data_total$event[data_total$loss_to_fu] <- rep(0, n_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 =
event * (enrollment[analysis_at_enrollnumber[i]] - enrollment <= EndofStudy & subject_enrolled) |
(1 - event) * (enrollment[analysis_at_enrollnumber[i]] - enrollment <= time & subject_enrolled) |
(subject_enrolled & loss_to_fu),
real_time =
ifelse(event * (enrollment[analysis_at_enrollnumber[i]] - enrollment <= EndofStudy & subject_enrolled) |
(1 - event) * (enrollment[analysis_at_enrollnumber[i]] - enrollment <= time & subject_enrolled),
enrollment[analysis_at_enrollnumber[i]] - enrollment + sd(time) / 10000, time),
real_event =
ifelse(event * (enrollment[analysis_at_enrollnumber[i]] - enrollment <= EndofStudy & subject_enrolled) |
(1 - event) * (enrollment[analysis_at_enrollnumber[i]] - enrollment <= time & subject_enrolled),
0, event))
data_interim <- data_interim %>%
mutate(time = real_time, event = real_event) %>%
select(-c(real_time, real_event))
# Carry out interim analysis on patients with complete data only
# - Set-up `new data` data frame
data <- data_interim %>%
filter(subject_enrolled) %>%
select(time, event, treatment)
# Analyze data using discount function via bdpsurvival
post <- suppressWarnings(
bdpsurvival(formula = Surv(time, event) ~ treatment,
data = data,
data0 = data0,
breaks = breaks,
a0 = prior[1],
b0 = prior[2],
surv_time = EndofStudy,
discount_function = discount_function,
alpha_max = alpha_max,
fix_alpha = fix_alpha,
number_mcmc = number_mcmc,
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) {
##########################################################################
### Expected success computations
##########################################################################
# Imputing the success for control group
if (!is.null(hazard_control)) {
control_impute <- data_interim %>%
filter(treatment == 0 & subject_impute_success)
impute_control <- pw_exp_impute(time = control_impute$time,
hazard = post$posterior_control$posterior_flat_hazard[i, ],
maxtime = EndofStudy,
cutpoint = post$args1$breaks)
data_control_success_impute <- data_interim %>%
filter(treatment == 0 & subject_impute_success) %>%
bind_cols(time_impute = impute_control$time,
event_impute = impute_control$event)
} else {
data_control_success_impute <- NULL
}
# Imputing success for treatment group
treatment_impute <- data_interim %>%
filter(treatment == 1 & subject_impute_success)
impute_treatment <- pw_exp_impute(time = treatment_impute$time,
hazard = post$posterior_treatment$posterior_flat_hazard[i, ],
maxtime = EndofStudy,
cutpoint = post$args1$breaks)
data_treatment_success_impute <- data_interim %>%
filter(treatment == 1 & subject_impute_success) %>%
bind_cols(time_impute = impute_treatment$time,
event_impute = impute_treatment$event)
# Combine the treatment and control imputed datasets
data_noimpute <- data_interim %>%
filter(!subject_impute_success) %>%
mutate(time_impute = time, event_impute = event)
# Combine the treatment and control imputed datasets
data_success_impute <- bind_rows(data_control_success_impute,
data_treatment_success_impute,
data_noimpute) %>%
mutate(time = time_impute, event = event_impute) %>%
select(-c(time_impute, event_impute))
# Create enrolled subject data frame for discount function analysis
data <- data_success_impute %>%
filter(subject_enrolled) %>%
ungroup() %>%
select(time, event, treatment)
# Analyze complete+imputed data using discount fucntion via bdpsurvival
post_imp <- suppressWarnings(
bdpsurvival(formula = Surv(time, event) ~ treatment,
data = data,
data0 = data0,
breaks = breaks,
a0 = prior[1],
b0 = prior[2],
surv_time = EndofStudy,
discount_function = discount_function,
alpha_max = alpha_max,
fix_alpha = fix_alpha,
number_mcmc = number_mcmc,
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) {
effect_imp <- post_imp$final$posterior_loghazard
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)
}
} else {
effect_imp <- post_imp$final$posterior_survival
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
# Imputing the control group
if (!is.null(hazard_control)) {
control_impute <- data_success_impute %>%
filter(treatment == 0 & subject_impute_futility)
impute_control <- pw_exp_sim(hazard = post$posterior_control$posterior_flat_hazard[i, ],
n = nrow(control_impute),
maxtime = EndofStudy,
cutpoint = post$args1$breaks)
data_control_futility_impute <- data_success_impute %>%
filter(treatment == 0 & subject_impute_futility) %>%
bind_cols(time_impute = impute_control$time,
event_impute = impute_control$event)
} else {
data_control_futility_impute <- NULL
}
# Imputing the treatment group
treatment_impute <- data_success_impute %>%
filter(treatment == 1 & subject_impute_futility)
impute_treatment <- pw_exp_sim(hazard = post$posterior_treatment$posterior_flat_hazard[i, ],
n = nrow(treatment_impute),
maxtime = EndofStudy,
cutpoint = post$args1$breaks)
data_treatment_futility_impute <- data_success_impute %>%
filter(treatment == 1 & subject_impute_futility) %>%
bind_cols(time_impute = impute_treatment$time,
event_impute = impute_treatment$event)
# Combine the treatment and control imputed datasets
data_noimpute_futility <- data_success_impute %>%
filter(!subject_impute_futility) %>%
mutate(time_impute = time, event_impute = event)
# Combine the treatment and control imputed datasets
data_futility_impute <- bind_rows(data_control_futility_impute,
data_treatment_futility_impute,
data_noimpute_futility) %>%
mutate(time = time_impute, event = event_impute) %>%
select(-c(time_impute, event_impute))
# Create enrolled subject data frame for discount function analysis
data <- data_futility_impute %>%
select(time, event, treatment)
# Analyze complete+imputed data using discount function via bdpsurvival
post_imp <- suppressWarnings(
bdpsurvival(formula = Surv(time, event) ~ treatment,
data = data,
data0 = data0,
breaks = breaks,
a0 = prior[1],
b0 = prior[2],
surv_time = EndofStudy,
discount_function = discount_function,
alpha_max = alpha_max,
fix_alpha = fix_alpha,
number_mcmc = number_mcmc,
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) {
effect_imp <- post_imp$final$posterior_loghazard
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)
}
} else {
effect_imp <- post_imp$final$posterior_survival
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
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 (sum(data$treatment == 0) != 0) {
if (alternative == "two-sided") {
effect_int <- post$final$posterior_loghazard
} else if (alternative == "greater") {
effect_int <- post$final$posterior_loghazard
} else {
effect_int <- post$final$posterior_loghazard
}
} else {
effect_int <- post_imp$posterior_treatment$posterior_survival
}
# 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) %>%
ungroup() %>%
select(time, event, treatment)
# Analyze complete data using discount function via bdpsurvival
post_final <- suppressWarnings(
bdpsurvival(formula = Surv(time, event) ~ treatment,
data = data_final,
data0 = data0,
breaks = breaks,
a0 = prior[1],
b0 = prior[2],
surv_time = EndofStudy,
discount_function = discount_function,
alpha_max = alpha_max,
fix_alpha = fix_alpha,
number_mcmc = number_mcmc,
weibull_scale = weibull_scale,
weibull_shape = weibull_shape,
method = method))
data_final <- data_total %>%
filter(id <= stage_trial_stopped)
# Format and output results
# Posterior effect size: test vs. control or treatment itself
if (sum(data_final$treatment == 0) != 0) {
effect <- post_final$final$posterior_loghazard
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)
}
} else {
effect <- post_final$final$posterior_survival
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
# Output
results_list <- list(
hazard_treatment = hazard_treatment, # Probability of treatment in binomial
hazard_control = hazard_control, # Probability of control in binomial
cutpoint = cutpoint,
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_enrolled = N_treatment + N_control,
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
)
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("time", "event", "event_impute", "real_time",
"time_impute", "id", "futility", "treatment",
"subject_impute_success", "real_event",
"subject_impute_futility"))
#' @title Analyzing Bayesian trial for time-to-event data
#'
#' @description Function to analyze Bayesian trial for time-to-event data which
#' allows early stopping and incorporation of historical data using the
#' discount function approach.
#'
#' @param time vector. exposure time for the subjects. It must be the same
#' length as the treatment variable.
#' @param event vector. The status indicator, normally 0=alive, 1=dead. Other
#' choices are TRUE/FALSE (TRUE = death) or 1/2 (2 = death). For censored
#' data, the status indicator is 0 = right censored, 1 = event at time.
#' Although unusual, the event indicator can be omitted, in which case all
#' subjects are assumed to have an event.
#' @param time0 vector. Historical exposure time for the subjects. It must be
#' the same length as the treatment variable.
#' @param event0 vector. Historical status indicator, normally 0=alive, 1=dead.
#' Other choices are TRUE/FALSE (TRUE = death) or 1/2 (2 = death). For
#' censored data, the status indicator is 0 = right censored, 1 = event at
#' time. Although unusual, the event indicator can be omitted, in which case
#' all subjects are assumed to have an event.
#' @param treatment0 vector. the historical treatment assignment for patients, 1
#' for treatment group and 0 for control group.
#' @param surv_time scalar. scalar. Survival time of interest for computing the
#' probability of survival for a single arm (OPC) trial. Default is overall
#' (pooled), i.e. current and historical, median survival time.
#' @param breaks vector. Breaks (interval starts) used to compose the breaks of
#' the piecewise exponential model. Do not include zero. Default breaks are
#' the quantiles of the input times.
#' @param prior vector. Prior values of the gamma rate, Gamma(a0, b0). The
#' default is set to Gamma(0.1, 0.1).
#' @inheritParams binomial_analysis
#' @inheritParams pw_exp_sim
#'
#' @return A list of output for the Bayesian trial for time-to-event:
#'
#' \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{alpha_max}}{
#' scalar. The alpha_max input.}
#' \item{\code{N_treatment}}{
#' scalar. The number of patients enrolled in the experimental group for
#' each simulation.}
#' \item{\code{event_treatment}}{
#' scalar. The number of events 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{event_control}}{
#' scalar. The number of events in the control group for
#' each simulation.}
#' \item{\code{N_enrolled}}{
#' scalar. 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 whose time passes the \code{surv_time}.}
#' \item{\code{alpha_discount}}{
#' vector. The alpha discount function used for control and treatment.}
#' \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.}
#' }
#'
#' @importFrom dplyr mutate bind_cols
#' @importFrom survival Surv
#' @importFrom dplyr mutate filter group_by bind_rows select n summarize
#' @importFrom bayesDP bdpsurvival
#' @importFrom stats median
#'
#' @export survival_analysis
survival_analysis <- function(
time,
treatment,
event = NULL,
time0 = NULL,
treatment0 = NULL,
event0 = NULL,
surv_time = NULL,
h0 = 0,
breaks = NULL,
alternative = "greater",
N_impute = 10,
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(event)) {
event <- rep(1, length(time))
}
# Combining the historical data
if (!is.null(time0)) {
data0 <- data.frame(cbind(time = time0, event = event0, treatment = treatment0))
} else {
data0 <- NULL
}
if (is.null(surv_time)) {
surv_time <- median(c(time, time0))
}
# Reading the data
data_total <- data.frame(cbind(time, event, treatment))
# Added interim data for incomplete data
data_interim <- data_total %>%
mutate(futility = (time < surv_time & event == 0))
# Analyze the data using bdpsurvival
post <- suppressWarnings(
bdpsurvival(formula = Surv(time, event) ~ treatment,
data = data_total,
data0 = data0,
breaks = breaks,
a0 = prior[1],
b0 = prior[2],
surv_time = surv_time,
discount_function = discount_function,
alpha_max = alpha_max,
fix_alpha = fix_alpha,
number_mcmc = number_mcmc,
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) {
##########################################################################
### Expected success computations
##########################################################################
control_impute <- data_interim %>%
filter(treatment == 0 & futility)
impute_control <- pw_exp_impute(time = control_impute$time,
hazard = post$posterior_control$posterior_flat_hazard[i, ],
maxtime = surv_time,
cutpoint = post$args1$breaks)
data_control_success_impute <- data_interim %>%
filter(treatment == 0 & futility) %>%
bind_cols(time_impute = impute_control$time,
event_impute = impute_control$event)
treatment_impute <- data_interim %>%
filter(treatment == 1 & futility)
impute_treatment <- pw_exp_impute(time = treatment_impute$time,
hazard = post$posterior_treatment$posterior_flat_hazard[i, ],
maxtime = surv_time,
cutpoint = post$args1$breaks)
data_treatment_success_impute <- data_interim %>%
filter(treatment == 1 & futility) %>%
bind_cols(time_impute = impute_treatment$time,
event_impute = impute_treatment$event)
data_noimpute <- data_interim %>%
filter(!futility) %>%
mutate(time_impute = time, event_impute = event)
# Combine the treatment and control imputed datasets
data_success_impute <- bind_rows(data_control_success_impute,
data_treatment_success_impute,
data_noimpute) %>%
mutate(time = time_impute, event = event_impute) %>%
select(-c(time_impute, event_impute))
# Analyze complete+imputed data using discount function via bdpsurvival
post_imp <- suppressWarnings(
bdpsurvival(formula = Surv(time, event) ~ treatment,
data = data_success_impute,
data0 = data0,
breaks = breaks,
a0 = prior[1],
b0 = prior[2],
surv_time = surv_time,
discount_function = discount_function,
alpha_max = alpha_max,
fix_alpha = fix_alpha,
number_mcmc = number_mcmc,
weibull_scale = weibull_scale,
weibull_shape = weibull_shape,
method = method))
# Posterior effect size: test vs. control or treatment itself
if (sum(data_success_impute$treatment == 0) != 0) {
effect_imp <- post_imp$final$posterior_loghazard
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)
}
} else {
effect <- post_imp$posterior_treatment$posterior_survival
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
}
##########################################################################
### Futility computations
##########################################################################
# TBA
data_final <- data_total
# Analyze complete data using discount function via binomial
post_final <- suppressWarnings(
bdpsurvival(formula = Surv(time, event) ~ treatment,
data = data_final,
data0 = data0,
breaks = breaks,
a0 = prior[1],
b0 = prior[2],
surv_time = surv_time,
discount_function = discount_function,
alpha_max = alpha_max,
fix_alpha = fix_alpha,
number_mcmc = number_mcmc,
weibull_scale = weibull_scale,
weibull_shape = weibull_shape,
method = method)
)
##############################################################################
### Summary at the interim analysis
##############################################################################
# Posterior effect size: test vs. control or treatment itself
if (sum(data_final$treatment == 0) != 0) {
effect <- post_final$final$posterior_loghazard
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)
}
} else {
effect <- post_final$posterior_treatment$posterior_survival
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 <- N_treatment + N_control
N_complete <- sum(!data_interim$futility)
alpha_discount <- c(post_final$posterior_control$alpha_discount,
post_final$posterior_treatment$alpha_discount)
event_treatment <- sum(data_final$event[!!data_final$treatment])
event_control <- sum(data_final$event[!data_final$treatment])
# Output
results_list <- list(
prob_of_accepting_alternative = prob_ha,
margin = h0, # Margin for error
alternative = alternative, # Alternative hypothesis
alpha_max = alpha_max,
N_treatment = N_treatment,
event_treatment = event_treatment,
N_control = N_control,
event_control = event_control,
N_enrolled = N_treatment + N_control,
N_complete = N_complete,
alpha_discount = alpha_discount,
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?
)
return(results_list)
}
## Quiets concerns of R CMD check re: the .'s that appear in pipelines
if (getRversion() >= "2.15.1") utils::globalVariables(c("time", "event", "event_impute",
"time_impute", "id", "futility", "treatment",
"subject_impute_success",
"subject_impute_futility"))
#' @title Historical data for survival analysis
#'
#' @description Wrapper function for historical data from time-to-event outcome.
#'
#' @inheritParams survival_analysis
#' @param .data NULL. Stores the historical time, treatment and event. Should
#' not be edited by the user.
#'
#' @return A list with historical data for time-to-event outcome with the
#' discount function.
#'
#' @examples
#' historical_survival(time = rexp(10, 0.01),
#' treatment = rep(10, 1),
#' event = rep(10, 1))
#'
#' @export historical_survival
historical_survival <- function(time = NULL,
treatment = NULL,
event = NULL,
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",
.data = NULL
) {
.data$time0 <- time
.data$treatment0 <- treatment
.data$event0 <- event
.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 Gamma prior for for control and treatment group
#'
#' @description Wrapper function for gamma prior \code{Gamma(a0, b0)}.
#'
#' @param a0 numeric. The shape parameter in the gamma distribution
#' (\code{beta(a0, b0)}).
#' @param b0 numeric. The scale parameter in the beta distribution
#' (\code{beta(a0, b0)}).
#' @param .data NULL. Stores the gamma prior rate. Should not be edited by the
#' user.
#'
#' @return A list with vector of gamma rate for the gamma prior for treatment
#' and control group.
#'
#' @examples gamma_prior(a0 = 0.1, b0 = 0.1)
#'
#' @export gamma_prior
gamma_prior <- function(a0 = 0.1, b0 = 0.1, .data = NULL) {
.data$prior <- c(a0, b0)
.data
}
#' @title Piecewise constant hazard rates and the cutpoint for control and
#' treatment group
#'
#' @description Wrapper function for the piecewise constant hazard rates and the
#' cutpoint for control and treatment group.
#'
#' @inheritParams survivalBACT
#' @param .data NULL. Stores the hazard rates and cutpoint. Should not be edited
#' by the user.
#'
#' @return A list with hazard rates and cutpoint for control and treatment
#' group.
#'
#' @examples survival_outcome(hazard_treatment = 0.06,
#' hazard_control = 0.08,
#' cutpoint = NULL)
#' @export survival_outcome
survival_outcome <- function(hazard_treatment = NULL,
cutpoint = NULL,
hazard_control = NULL,
.data = NULL) {
.data$hazard_treatment <- hazard_treatment
.data$cutpoint <- cutpoint
.data$hazard_control <- hazard_control
.data
}
#' @title Data file for survival analysis
#'
#' @description Wrapper function for data file in survival analysis.
#'
#' @inheritParams survival_analysis
#' @param .data NULL. Stores the survival data for analysis. Should not be
#' edited by the user.
#'
#' @return A list with time, treatment, and event with time-to-event outcome.
#'
#' @examples data_survival(time = c(6.2, 8.2, 8.0, 2.3),
#' treatment = c(0, 1, 0, 1),
#' event = c(1, 1, 1, 1))
#'
#' @export data_survival
data_survival <- function(time, treatment, event, .data = NULL) {
.data$time <- time
.data$treatment <- treatment
.data$event <- event
.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_survival survival_outcome gamma_prior historical_survival survival_analysis survivalBACT
Documented in data_survival gamma_prior historical_survival survival_analysis survivalBACT survival_outcome
#' @title Time-to-event outcome for Bayesian Adaptive trials
#'
#' @description Simulation for time-to-event outcome for Bayesian Adaptive trial with
#' different inputs to control for power, sample size, type 1 error rate, etc.
#'
#' @param hazard_treatment vector. Constant hazard rates under the treatment arm.
#' @param hazard_control vector. Constant hazard rates under the control arm.
#' @inheritParams pw_exp_sim
#' @inheritParams normalBACT
#' @inheritParams survival_analysis
#'
#' @return A list of output for a single trial simulation.
#' \describe{
#' \item{\code{lambda_treatment}}{
#' vector. The input parameter of constant hazard rates in the
#' treatment group.}
#' \item{\code{cutpoint_treatment}}{
#' vector. The change-point vector when the constant hazard rate(s) changes for
#' the treatment group.}
#' \item{\code{lambda_control}}{
#' vector. The input parameter of constant hazard rates in the
#' control group.}
#' \item{\code{cutpoint_control}}{
#' vector. The change-point vector when the constant hazard rate(s) changes for
#' 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 treatment group for
#' each simulation.}
#' \item{\code{event_treatment}}{
#' scalar. The number of events in the treatment group for
#' each simulation.}
#' \item{\code{N_control}}{
#' scalar. The number of patients enrolled in the control group for
#' each simulation.}
#' \item{\code{event_control}}{
#' scalar. The number of events 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 runif
#' @importFrom dplyr mutate filter group_by bind_cols bind_rows select ungroup
#' @importFrom bayesDP bdpsurvival
#'
#' @export survivalBACT
survivalBACT <- function(
hazard_treatment,
cutpoint = NULL,
hazard_control = NULL,
N_total,
breaks = NULL,
time0 = NULL,
treatment0 = NULL,
event0 = NULL,
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", # 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 interim_look
if (!is.null(interim_look)) {
stopifnot(all(N_total > interim_look))
}
# Combining the histrotical data
if (!is.null(time0)) {
data0 <- data.frame(cbind(time = time0, event = event0, treatment = treatment0))
} else {
data0 <- NULL
}
# Checking if alternative is right
if (alternative != "two-sided" & alternative != "greater" & alternative != "less") {
stop("The input for alternative is wrong!")
}
# If cutpoint is not NULL, assign breaks to cutpoint
if (!is.null(cutpoint)) {
breaks <- cutpoint
}
# Make sure breaks is not more than EndofStudy
if (!is.null(breaks)) {
stopifnot(any(breaks < EndofStudy))
}
if (is.null(hazard_control) & h0 == 0) {
h0 <- 0.5
}
# 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(hazard_control)) {
group <- randomization(N_total = N_total, block = block, allocation = rand_ratio)
} else {
group <- rep(1, N_total)
}
time <- rep(NA, length = N_total)
event <- rep(NA, length = N_total)
# Simulate TTE outcome
if (!is.null(hazard_control)) {
sim_control <- pw_exp_sim(hazard = hazard_control,
n = sum(!group),
maxtime = EndofStudy,
cutpoint = cutpoint)
time[which(!group)] <- sim_control$time
event[which(!group)] <- sim_control$event
}
sim_treatment <- pw_exp_sim(hazard = hazard_treatment,
n = sum(group),
maxtime = EndofStudy,
cutpoint = cutpoint)
time[which(!!group)] <- sim_treatment$time
event[which(!!group)] <- sim_treatment$event
if (is.null(breaks)) {
breaks <- mean(time)
}
# 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(
time = time,
treatment = group,
event = event,
enrollment = enrollment,
id = 1:N_total,
loss_to_fu = loss_to_fu)
# Subjects lost are uniformly distributed
data_total$time[data_total$loss_to_fu] <- runif(n_loss_to_fu,
0, data_total$time[data_total$loss_to_fu])
data_total$event[data_total$loss_to_fu] <- rep(0, n_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 =
event * (enrollment[analysis_at_enrollnumber[i]] - enrollment <= EndofStudy & subject_enrolled) |
(1 - event) * (enrollment[analysis_at_enrollnumber[i]] - enrollment <= time & subject_enrolled) |
(subject_enrolled & loss_to_fu),
real_time =
ifelse(event * (enrollment[analysis_at_enrollnumber[i]] - enrollment <= EndofStudy & subject_enrolled) |
(1 - event) * (enrollment[analysis_at_enrollnumber[i]] - enrollment <= time & subject_enrolled),
enrollment[analysis_at_enrollnumber[i]] - enrollment + sd(time) / 10000, time),
real_event =
ifelse(event * (enrollment[analysis_at_enrollnumber[i]] - enrollment <= EndofStudy & subject_enrolled) |
(1 - event) * (enrollment[analysis_at_enrollnumber[i]] - enrollment <= time & subject_enrolled),
0, event))
data_interim <- data_interim %>%
mutate(time = real_time, event = real_event) %>%
select(-c(real_time, real_event))
# Carry out interim analysis on patients with complete data only
# - Set-up `new data` data frame
data <- data_interim %>%
filter(subject_enrolled) %>%
select(time, event, treatment)
# Analyze data using discount function via bdpsurvival
post <- suppressWarnings(
bdpsurvival(formula = Surv(time, event) ~ treatment,
data = data,
data0 = data0,
breaks = breaks,
a0 = prior[1],
b0 = prior[2],
surv_time = EndofStudy,
discount_function = discount_function,
alpha_max = alpha_max,
fix_alpha = fix_alpha,
number_mcmc = number_mcmc,
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) {
##########################################################################
### Expected success computations
##########################################################################
# Imputing the success for control group
if (!is.null(hazard_control)) {
control_impute <- data_interim %>%
filter(treatment == 0 & subject_impute_success)
impute_control <- pw_exp_impute(time = control_impute$time,
hazard = post$posterior_control$posterior_flat_hazard[i, ],
maxtime = EndofStudy,
cutpoint = post$args1$breaks)
data_control_success_impute <- data_interim %>%
filter(treatment == 0 & subject_impute_success) %>%
bind_cols(time_impute = impute_control$time,
event_impute = impute_control$event)
} else {
data_control_success_impute <- NULL
}
# Imputing success for treatment group
treatment_impute <- data_interim %>%
filter(treatment == 1 & subject_impute_success)
impute_treatment <- pw_exp_impute(time = treatment_impute$time,
hazard = post$posterior_treatment$posterior_flat_hazard[i, ],
maxtime = EndofStudy,
cutpoint = post$args1$breaks)
data_treatment_success_impute <- data_interim %>%
filter(treatment == 1 & subject_impute_success) %>%
bind_cols(time_impute = impute_treatment$time,
event_impute = impute_treatment$event)
# Combine the treatment and control imputed datasets
data_noimpute <- data_interim %>%
filter(!subject_impute_success) %>%
mutate(time_impute = time, event_impute = event)
# Combine the treatment and control imputed datasets
data_success_impute <- bind_rows(data_control_success_impute,
data_treatment_success_impute,
data_noimpute) %>%
mutate(time = time_impute, event = event_impute) %>%
select(-c(time_impute, event_impute))
# Create enrolled subject data frame for discount function analysis
data <- data_success_impute %>%
filter(subject_enrolled) %>%
ungroup() %>%
select(time, event, treatment)
# Analyze complete+imputed data using discount fucntion via bdpsurvival
post_imp <- suppressWarnings(
bdpsurvival(formula = Surv(time, event) ~ treatment,
data = data,
data0 = data0,
breaks = breaks,
a0 = prior[1],
b0 = prior[2],
surv_time = EndofStudy,
discount_function = discount_function,
alpha_max = alpha_max,
fix_alpha = fix_alpha,
number_mcmc = number_mcmc,
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) {
effect_imp <- post_imp$final$posterior_loghazard
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)
}
} else {
effect_imp <- post_imp$final$posterior_survival
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
# Imputing the control group
if (!is.null(hazard_control)) {
control_impute <- data_success_impute %>%
filter(treatment == 0 & subject_impute_futility)
impute_control <- pw_exp_sim(hazard = post$posterior_control$posterior_flat_hazard[i, ],
n = nrow(control_impute),
maxtime = EndofStudy,
cutpoint = post$args1$breaks)
data_control_futility_impute <- data_success_impute %>%
filter(treatment == 0 & subject_impute_futility) %>%
bind_cols(time_impute = impute_control$time,
event_impute = impute_control$event)
} else {
data_control_futility_impute <- NULL
}
# Imputing the treatment group
treatment_impute <- data_success_impute %>%
filter(treatment == 1 & subject_impute_futility)
impute_treatment <- pw_exp_sim(hazard = post$posterior_treatment$posterior_flat_hazard[i, ],
n = nrow(treatment_impute),
maxtime = EndofStudy,
cutpoint = post$args1$breaks)
data_treatment_futility_impute <- data_success_impute %>%
filter(treatment == 1 & subject_impute_futility) %>%
bind_cols(time_impute = impute_treatment$time,
event_impute = impute_treatment$event)
# Combine the treatment and control imputed datasets
data_noimpute_futility <- data_success_impute %>%
filter(!subject_impute_futility) %>%
mutate(time_impute = time, event_impute = event)
# Combine the treatment and control imputed datasets
data_futility_impute <- bind_rows(data_control_futility_impute,
data_treatment_futility_impute,
data_noimpute_futility) %>%
mutate(time = time_impute, event = event_impute) %>%
select(-c(time_impute, event_impute))
# Create enrolled subject data frame for discount function analysis
data <- data_futility_impute %>%
select(time, event, treatment)
# Analyze complete+imputed data using discount function via bdpsurvival
post_imp <- suppressWarnings(
bdpsurvival(formula = Surv(time, event) ~ treatment,
data = data,
data0 = data0,
breaks = breaks,
a0 = prior[1],
b0 = prior[2],
surv_time = EndofStudy,
discount_function = discount_function,
alpha_max = alpha_max,
fix_alpha = fix_alpha,
number_mcmc = number_mcmc,
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) {
effect_imp <- post_imp$final$posterior_loghazard
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)
}
} else {
effect_imp <- post_imp$final$posterior_survival
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
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 (sum(data$treatment == 0) != 0) {
if (alternative == "two-sided") {
effect_int <- post$final$posterior_loghazard
} else if (alternative == "greater") {
effect_int <- post$final$posterior_loghazard
} else {
effect_int <- post$final$posterior_loghazard
}
} else {
effect_int <- post_imp$posterior_treatment$posterior_survival
}
# 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) %>%
ungroup() %>%
select(time, event, treatment)
# Analyze complete data using discount function via bdpsurvival
post_final <- suppressWarnings(
bdpsurvival(formula = Surv(time, event) ~ treatment,
data = data_final,
data0 = data0,
breaks = breaks,
a0 = prior[1],
b0 = prior[2],
surv_time = EndofStudy,
discount_function = discount_function,
alpha_max = alpha_max,
fix_alpha = fix_alpha,
number_mcmc = number_mcmc,
weibull_scale = weibull_scale,
weibull_shape = weibull_shape,
method = method))
data_final <- data_total %>%
filter(id <= stage_trial_stopped)
# Format and output results
# Posterior effect size: test vs. control or treatment itself
if (sum(data_final$treatment == 0) != 0) {
effect <- post_final$final$posterior_loghazard
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)
}
} else {
effect <- post_final$final$posterior_survival
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
# Output
results_list <- list(
hazard_treatment = hazard_treatment, # Probability of treatment in binomial
hazard_control = hazard_control, # Probability of control in binomial
cutpoint = cutpoint,
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_enrolled = N_treatment + N_control,
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
)
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("time", "event", "event_impute", "real_time",
"time_impute", "id", "futility", "treatment",
"subject_impute_success", "real_event",
"subject_impute_futility"))
#' @title Analyzing Bayesian trial for time-to-event data
#'
#' @description Function to analyze Bayesian trial for time-to-event data which
#' allows early stopping and incorporation of historical data using the
#' discount function approach.
#'
#' @param time vector. exposure time for the subjects. It must be the same
#' length as the treatment variable.
#' @param event vector. The status indicator, normally 0=alive, 1=dead. Other
#' choices are TRUE/FALSE (TRUE = death) or 1/2 (2 = death). For censored
#' data, the status indicator is 0 = right censored, 1 = event at time.
#' Although unusual, the event indicator can be omitted, in which case all
#' subjects are assumed to have an event.
#' @param time0 vector. Historical exposure time for the subjects. It must be
#' the same length as the treatment variable.
#' @param event0 vector. Historical status indicator, normally 0=alive, 1=dead.
#' Other choices are TRUE/FALSE (TRUE = death) or 1/2 (2 = death). For
#' censored data, the status indicator is 0 = right censored, 1 = event at
#' time. Although unusual, the event indicator can be omitted, in which case
#' all subjects are assumed to have an event.
#' @param treatment0 vector. the historical treatment assignment for patients, 1
#' for treatment group and 0 for control group.
#' @param surv_time scalar. scalar. Survival time of interest for computing the
#' probability of survival for a single arm (OPC) trial. Default is overall
#' (pooled), i.e. current and historical, median survival time.
#' @param breaks vector. Breaks (interval starts) used to compose the breaks of
#' the piecewise exponential model. Do not include zero. Default breaks are
#' the quantiles of the input times.
#' @param prior vector. Prior values of the gamma rate, Gamma(a0, b0). The
#' default is set to Gamma(0.1, 0.1).
#' @inheritParams binomial_analysis
#' @inheritParams pw_exp_sim
#'
#' @return A list of output for the Bayesian trial for time-to-event:
#'
#' \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{alpha_max}}{
#' scalar. The alpha_max input.}
#' \item{\code{N_treatment}}{
#' scalar. The number of patients enrolled in the experimental group for
#' each simulation.}
#' \item{\code{event_treatment}}{
#' scalar. The number of events 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{event_control}}{
#' scalar. The number of events in the control group for
#' each simulation.}
#' \item{\code{N_enrolled}}{
#' scalar. 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 whose time passes the \code{surv_time}.}
#' \item{\code{alpha_discount}}{
#' vector. The alpha discount function used for control and treatment.}
#' \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.}
#' }
#'
#' @importFrom dplyr mutate bind_cols
#' @importFrom survival Surv
#' @importFrom dplyr mutate filter group_by bind_rows select n summarize
#' @importFrom bayesDP bdpsurvival
#' @importFrom stats median
#'
#' @export survival_analysis
survival_analysis <- function(
time,
treatment,
event = NULL,
time0 = NULL,
treatment0 = NULL,
event0 = NULL,
surv_time = NULL,
h0 = 0,
breaks = NULL,
alternative = "greater",
N_impute = 10,
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(event)) {
event <- rep(1, length(time))
}
# Combining the historical data
if (!is.null(time0)) {
data0 <- data.frame(cbind(time = time0, event = event0, treatment = treatment0))
} else {
data0 <- NULL
}
if (is.null(surv_time)) {
surv_time <- median(c(time, time0))
}
# Reading the data
data_total <- data.frame(cbind(time, event, treatment))
# Added interim data for incomplete data
data_interim <- data_total %>%
mutate(futility = (time < surv_time & event == 0))
# Analyze the data using bdpsurvival
post <- suppressWarnings(
bdpsurvival(formula = Surv(time, event) ~ treatment,
data = data_total,
data0 = data0,
breaks = breaks,
a0 = prior[1],
b0 = prior[2],
surv_time = surv_time,
discount_function = discount_function,
alpha_max = alpha_max,
fix_alpha = fix_alpha,
number_mcmc = number_mcmc,
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) {
##########################################################################
### Expected success computations
##########################################################################
control_impute <- data_interim %>%
filter(treatment == 0 & futility)
impute_control <- pw_exp_impute(time = control_impute$time,
hazard = post$posterior_control$posterior_flat_hazard[i, ],
maxtime = surv_time,
cutpoint = post$args1$breaks)
data_control_success_impute <- data_interim %>%
filter(treatment == 0 & futility) %>%
bind_cols(time_impute = impute_control$time,
event_impute = impute_control$event)
treatment_impute <- data_interim %>%
filter(treatment == 1 & futility)
impute_treatment <- pw_exp_impute(time = treatment_impute$time,
hazard = post$posterior_treatment$posterior_flat_hazard[i, ],
maxtime = surv_time,
cutpoint = post$args1$breaks)
data_treatment_success_impute <- data_interim %>%
filter(treatment == 1 & futility) %>%
bind_cols(time_impute = impute_treatment$time,
event_impute = impute_treatment$event)
data_noimpute <- data_interim %>%
filter(!futility) %>%
mutate(time_impute = time, event_impute = event)
# Combine the treatment and control imputed datasets
data_success_impute <- bind_rows(data_control_success_impute,
data_treatment_success_impute,
data_noimpute) %>%
mutate(time = time_impute, event = event_impute) %>%
select(-c(time_impute, event_impute))
# Analyze complete+imputed data using discount function via bdpsurvival
post_imp <- suppressWarnings(
bdpsurvival(formula = Surv(time, event) ~ treatment,
data = data_success_impute,
data0 = data0,
breaks = breaks,
a0 = prior[1],
b0 = prior[2],
surv_time = surv_time,
discount_function = discount_function,
alpha_max = alpha_max,
fix_alpha = fix_alpha,
number_mcmc = number_mcmc,
weibull_scale = weibull_scale,
weibull_shape = weibull_shape,
method = method))
# Posterior effect size: test vs. control or treatment itself
if (sum(data_success_impute$treatment == 0) != 0) {
effect_imp <- post_imp$final$posterior_loghazard
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)
}
} else {
effect <- post_imp$posterior_treatment$posterior_survival
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
}
##########################################################################
### Futility computations
##########################################################################
# TBA
data_final <- data_total
# Analyze complete data using discount function via binomial
post_final <- suppressWarnings(
bdpsurvival(formula = Surv(time, event) ~ treatment,
data = data_final,
data0 = data0,
breaks = breaks,
a0 = prior[1],
b0 = prior[2],
surv_time = surv_time,
discount_function = discount_function,
alpha_max = alpha_max,
fix_alpha = fix_alpha,
number_mcmc = number_mcmc,
weibull_scale = weibull_scale,
weibull_shape = weibull_shape,
method = method)
)
##############################################################################
### Summary at the interim analysis
##############################################################################
# Posterior effect size: test vs. control or treatment itself
if (sum(data_final$treatment == 0) != 0) {
effect <- post_final$final$posterior_loghazard
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)
}
} else {
effect <- post_final$posterior_treatment$posterior_survival
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 <- N_treatment + N_control
N_complete <- sum(!data_interim$futility)
alpha_discount <- c(post_final$posterior_control$alpha_discount,
post_final$posterior_treatment$alpha_discount)
event_treatment <- sum(data_final$event[!!data_final$treatment])
event_control <- sum(data_final$event[!data_final$treatment])
# Output
results_list <- list(
prob_of_accepting_alternative = prob_ha,
margin = h0, # Margin for error
alternative = alternative, # Alternative hypothesis
alpha_max = alpha_max,
N_treatment = N_treatment,
event_treatment = event_treatment,
N_control = N_control,
event_control = event_control,
N_enrolled = N_treatment + N_control,
N_complete = N_complete,
alpha_discount = alpha_discount,
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?
)
return(results_list)
}
## Quiets concerns of R CMD check re: the .'s that appear in pipelines
if (getRversion() >= "2.15.1") utils::globalVariables(c("time", "event", "event_impute",
"time_impute", "id", "futility", "treatment",
"subject_impute_success",
"subject_impute_futility"))
#' @title Historical data for survival analysis
#'
#' @description Wrapper function for historical data from time-to-event outcome.
#'
#' @inheritParams survival_analysis
#' @param .data NULL. Stores the historical time, treatment and event. Should
#' not be edited by the user.
#'
#' @return A list with historical data for time-to-event outcome with the
#' discount function.
#'
#' @examples
#' historical_survival(time = rexp(10, 0.01),
#' treatment = rep(10, 1),
#' event = rep(10, 1))
#'
#' @export historical_survival
historical_survival <- function(time = NULL,
treatment = NULL,
event = NULL,
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",
.data = NULL
) {
.data$time0 <- time
.data$treatment0 <- treatment
.data$event0 <- event
.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 Gamma prior for for control and treatment group
#'
#' @description Wrapper function for gamma prior \code{Gamma(a0, b0)}.
#'
#' @param a0 numeric. The shape parameter in the gamma distribution
#' (\code{beta(a0, b0)}).
#' @param b0 numeric. The scale parameter in the beta distribution
#' (\code{beta(a0, b0)}).
#' @param .data NULL. Stores the gamma prior rate. Should not be edited by the
#' user.
#'
#' @return A list with vector of gamma rate for the gamma prior for treatment
#' and control group.
#'
#' @examples gamma_prior(a0 = 0.1, b0 = 0.1)
#'
#' @export gamma_prior
gamma_prior <- function(a0 = 0.1, b0 = 0.1, .data = NULL) {
.data$prior <- c(a0, b0)
.data
}
#' @title Piecewise constant hazard rates and the cutpoint for control and
#' treatment group
#'
#' @description Wrapper function for the piecewise constant hazard rates and the
#' cutpoint for control and treatment group.
#'
#' @inheritParams survivalBACT
#' @param .data NULL. Stores the hazard rates and cutpoint. Should not be edited
#' by the user.
#'
#' @return A list with hazard rates and cutpoint for control and treatment
#' group.
#'
#' @examples survival_outcome(hazard_treatment = 0.06,
#' hazard_control = 0.08,
#' cutpoint = NULL)
#' @export survival_outcome
survival_outcome <- function(hazard_treatment = NULL,
cutpoint = NULL,
hazard_control = NULL,
.data = NULL) {
.data$hazard_treatment <- hazard_treatment
.data$cutpoint <- cutpoint
.data$hazard_control <- hazard_control
.data
}
#' @title Data file for survival analysis
#'
#' @description Wrapper function for data file in survival analysis.
#'
#' @inheritParams survival_analysis
#' @param .data NULL. Stores the survival data for analysis. Should not be
#' edited by the user.
#'
#' @return A list with time, treatment, and event with time-to-event outcome.
#'
#' @examples data_survival(time = c(6.2, 8.2, 8.0, 2.3),
#' treatment = c(0, 1, 0, 1),
#' event = c(1, 1, 1, 1))
#'
#' @export data_survival
data_survival <- function(time, treatment, event, .data = NULL) {
.data$time <- time
.data$treatment <- treatment
.data$event <- event
.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.