R/Command_Line_Version_of_Adaptive_Group_Sequential_Design.R

In interAdapt: interAdapt

#' Table of default parameters
#'
#' A table of the default parameters entered to \code{\link{compute_design_performance}}.
#'
#' @details Each row of this table corresponds to one input to \code{\link{compute_design_performance}}, and contains: (1) the name of the argument passed to \code{\link{compute_design_performance}}, (2) a label describing the input in more detail, (3) the minimum value the input can take, (4) the maximum value the input can take, and (5) the input's default value.
#'
#'
#' @name default_parameter_table
#' @docType data
NULL

#' Compute Design Performance
#'
#' @description
#' Generates decision rules for group
#' sequential trial designs with adaptive enrollment criteria. Tables are
#' also generated which compare the performance of these designs to the
#' performance of standard group sequential designs. We use the notation \code{AD} to refer to the design with adaptive enrollment, \code{SC} to refer to a standard group sequential design enrolling from the combined population, and \code{SS} to refer to a standard group sequential design enrolling from only the subpopulation where there is greater prior evidence of a positive treatment effect. Further details are provided in (Fisher et al. 2014).
#' 
#' 
#' 
#' 
#' @param p_1 The proportion of the population in subpopulation 1,
#' which is the subpopulation having stronger prior evidence of a positive treatment effect.
#' 
#' @param p_10 The probability of a successful outcome for subpopulation 1 under assignment
#' to the control arm. This is used in estimating power and expected sample size.
#' 
#' @param p_20 The probability of a successful outcome for subpopulation 2 under assignment
#' to the control arm. This is used in estimating power and expected sample size.
#' 
#' @param p_11 The probability of a successful outcome for subpopulation 1 under assignment
#' to the treatment arm. Note that the user does not specify the probability of success under treatment for subpopulation 2 (\eqn{p_2t}).
#' Instead, \code{compute_design_performance} considers a range of possible values of \eqn{p_2t} (see the \code{lower_bound_treatment_effect_subpop_2} and \code{upper_bound_treatment_effect_subpop_2} arguments).
#' 
#' @param per_stage_sample_size_combined_AD The
#' number of participants enrolled per stage in the adaptive design while
#' both subpopulations are being enrolled.
#' 
#' @param per_stage_sample_size_subpop_1_AD The number of participants required for each stage
#' in the adaptive design after enrollment for subpopulation 2 has been
#' stopped.
#' 
#' @param FWER The familywise Type I error rate (\eqn{\alpha}) for all designs (see Fisher et al. 2014).
#' 
#' @param FWER_H0C_proportion
#' Proportion of \eqn{\alpha} allocated to \eqn{H_{0C}} for the adaptive design. Here, \eqn{H_{0C}} refers to the null hypothesis of no treatment effect in the combined population.
#' 
#' @param Delta Used as the exponent in defining the efficacy and
#' futility boundaries, as described (Fisher et al. 2014).
#' 
#' @param iter The number of simulated trials used to
#' estimate the power, expected sample size, and expected trial duration.
#' 
#' @param time_limit
#' Time limit for simulation in seconds. If the simulation
#' exceeds the time limit, calculations will stop and the user will get an
#' error message ("reached CPU time limit"). See \code{\link[base]{setTimeLimit}}.
#' To avoid this, the number of iterations (\code{iter}) can be reduced or
#' the time limit can be increased.
#' 
#' @param total_number_stages
#' Total number of stages
#' used in each design  (\eqn{K}).  The maximum allowable number of stages is 20.
#' 
#' @param last_stage_subpop_2_enrolled_AD
#' The last stage subpopulation 2 is enrolled, under the adaptive design. We refer to this stage number as \eqn{k*}.
#' 
#' @param enrollment_rate_combined_population The assumed
#' enrollment rate per year for the combined population.  This impacts the
#' expected duration of each trial design. Active enrollments from
#' the two subpopulations are assumed to be independent.  The enrollment rates
#' for subpopulations 1 and 2 are assumed proportional, based on \code{p_1}.
#' This implies that each stage of the adaptive design up to and including stage \code{k*} takes the same amount of time to complete, regardless of whether or not enrollment stops for subpopulation 2.  Each stage after \code{k*} will also take the same amount of time to complete. 
#' 
#' @param per_stage_sample_size_combined_SC
#' The number of participants enrolled in each stage of the standard group sequential design enrolling the
#' combined population (\code{SC}).
#' 
#' @param per_stage_sample_size_SS
#' The number of participants enrolled in each stage for standard group sequential design enrolling only
#' subpopulation 1 (\code{SS}).
#' 
#' @param subpop_2_stop_boundary_constant_AD
#' Stopping boundary proportionality constant for subpopulation 2 enrollment
#' in the adaptive design.
#' 
#' @param H01_futility_boundary_constant_AD
#' Futility boundary proportionality constant for \eqn{H_{01}} in the adaptive design. This is used to calculate the futility boundaries
#' (\eqn{l_{1,k}}) for the z-statistics calculated in subpopulation 1
#' (\eqn{Z_{1,k}}) as defined in (Fisher et al. 2014).
#' 
#' @param H0C_futility_boundary_constant_SC
#' Futility boundary proportionality constant for \eqn{H_{0C}} in the standard design always enrolling from the combined population.
#' 
#' @param H01_futility_boundary_constant_SS
#' Futility boundary proportionality constant for \eqn{H_{01}} in the standard design only enrolling from subpopulation 1.
#' 
#' @param lower_bound_treatment_effect_subpop_2
#' Simulations
#' are performed under a range of treatment effect sizes for
#' subpopulation 2 (i.e. \code{p_{2t}-p_{2c}}).  This parameter sets the lower bound for this range.
#' This effectively sets the lower bound for the probability of success under treatment for subpopulation 2 (\code{p_21}), since \code{p_20} is set
#' by the user.
#' 
#' @param upper_bound_treatment_effect_subpop_2 Simulations
#' are performed under a range of treatment effect sizes for
#' subpopulation 2 (i.e. \code{p_{2t}-p_{2c}}). This parameter sets the upper bound for this range.
#' 
#' @param CSV
#' Rather than manually entering the arguments above, this allows for the arguments to be entered via a tabular csv file. The \code{CSV} argument should contain a character vector or list of csv filenames. The
#' table must minimally include the columns "inputId" and "value" (as in the
#' \code{\link[interAdapt]{default_parameter_table}}).
#' 
#' 
#' 
#' @details
#' This function is meant to be applied when there is prior
#' evidence that a treatment might work better in a one subpopulation
#' than in another. In this context, a trial with an adaptive enrollment
#' criteria would determine whether or not to continue enrolling patients
#' from each subpopulation based on interim analyses of whether each
#' subpopulation is benefiting. In order for the type I error and the
#' power of the trial to be calculable, the decision rules for changing
#' enrollment must be set before the trial starts. This function generates
#' decision rules for group sequential trial designs with adaptive
#' enrollment criteria, and compares the performance of these designs
#' against standard group sequential designs with fixed enrollment criteria. Performance is compared in
#' terms of power, expected sample size, and expected trial duration.
#' 
#' 
#' @return
#' A list with 5 elements:
#' 
#' 		\item{performance_comparison}{A table comparing the performance of the three trials, in terms of power, expected sample size, and expected duration. See examples.}
#' 		\item{AD_design}{Efficacy and futility boundaries for the group sequential design with adaptive enrollment}
#' 		\item{SC_design}{Efficacy and futility boundaries for the standard group sequential design enrolling from the combined population}
#' 		\item{SS_design}{Efficacy and futility boundaries for the standard group sequential design enrolling subpopulation 1 only}
#' 		\item{input_parameters}{List of input argument values}
#' 
#' 
#' @references
#' Aaron Fisher, Harris Jaffee, and Michael Rosenblum. interAdapt -- An
#' Interactive Tool for Designing and Evaluating Randomized Trials with
#' Adaptive Enrollment Criteria. Working Paper, 2014.
#' http://arxiv.org/abs/1404.0734
#' 
#' @export
#' @import mvtnorm
#' 
#' @keywords Adaptive clinical trial, Adaptive enrollment, Group sequential designs
#' @examples
#' 
#' #Store function output
#' o<-compute_design_performance()
#' 
#' names(o)
#' 
#' #Decision rules for trial designs
#' print(o$SS[[1]])
#' print(o$SC[[1]])
#' print(o$AD[[1]])
#' 
#' #Plot decision rules
#' par(mfrow=c(1,3))
#' matplot(t(o$SS[[1]][2:3,]),type='o',
#'		main='Standard trial - subpop 1',
#'		xlab='stage',ylab='Z-statistic',
#'		col='blue',pch=2:3,lty=3)
#' legend('topright',c('H01 Efficacy','H01 Futility')
#'		,col='blue',pch=2:3,lty=3)
#' matplot(t(o$SC[[1]][4:5,]),type='o',
#'		main='Standard trial - combined pop',xlab='stage',
#'		ylab='Z-statistic',col='red',pch=0:1,lty=3)
#' legend('topright',c('H0C Efficacy','H0C Futility'),
#'		col='red',pch=0:1,lty=3)
#' matplot(t(o$AD[[1]][4:7,]),type='o',
#'		main='Adaptive Enrollment',xlab='stage',
#'		ylab='Z-statistic',col=c('red','red','blue','blue'),
#'		pch=0:3,lty=3)
#' legend('topright',c('H0C Efficacy','H0C Futility',
#'		'H01 Efficacy','H01 Futility'),col=c('red',
#'		'red','blue','blue'),pch=0:3,lty=3)
#' 
#' #Check performance
#' o$performance[[1]]
#' 
#' 
#' #Plot performance over a range of treatment effects for subpop2
#' col1<-c('black','black','black','green','blue')
#' 
#' 
#' perform_names<-rownames(o$performance[[1]])
#' 
#' #index for parts of the table corresponding to power
#' p_ind<- grep('Power',perform_names)
#' #index for parts of the table corresponding to sample size
#' s_ind<- grep('Sample Size',perform_names)
#' #index for parts of the table corresponding to trial duration
#' d_ind<- grep('DUR',perform_names)
#' 
#' 
#' par(mfrow=c(1,3))
#' 
#' lty1<-c(2,3,1,4,5)
#' matplot(x=t(o$performance[[1]][1,]),t(o$performance[[1]][p_ind,]),
#'		type='l',lty=lty1,col=col1[1:5],xlab='Subpop.2 Tx. Effect',
#'		ylab='Power',main='Power')
#' legend('bottomleft',perform_names[p_ind],col=col1[1:length(p_ind)],
#'		lty=lty1)
#' 
#' 
#' matplot(x=t(o$performance[[1]][1,]),t(o$performance[[1]][s_ind,]),
#'		type='l',lty=1:3,col=col1[3:5],xlab='Subpop.2 Tx. Effect',
#'		ylab='Expected Sample Size',main='Expected Sample Size')
#' legend('topright',substr(perform_names[d_ind],1,2),col=col1[3:5],
#'		lty=1:3)
#' 
#' 
#' matplot(x=t(o$performance[[1]][1,]),t(o$performance[[1]][d_ind,]),
#'		type='l',lty=1:3,col=col1[3:5],xlab='Subpop.2 Tx. Effect',
#'		ylab="Expected Duration",main='Expected Duration')
#' legend('topright',substr(perform_names[d_ind],1,2),col=col1[3:5],
#'		lty=1:3)
#' 
#'
compute_design_performance <- function(

## Subpopulation 1 proportion (Range: 0 to 1)
p_1 = 0.33,

## Note: throughout, we denote the treatment arm by A=1 and control arm by A=0. 
## We represent p_{1c} (probability of successful outcome Y=1 under assignment to control, for subpopulation s) by ps0, for each s = 1, 2.
## Similarly, we represent p_{1t} (probability of successful outcome Y=1 under assignment to treatment, for subpopulation s) by ps1, for each s = 1, 2.

## Probability outcome = 1 under control:
## for Subpopulation 1 (Range: 0 to 1)
p_10 = 0.25,
## for Subpopulation 2 (Range: 0 to 1)
p_20 = 0.20,

## Prob. outcome = 1 under treatment, at alternative:
## for Subpopulation 1 (Range: 0 to 1)
p_11 = p_10 + 0.125,
## The user does not input the corresponding value for Subpopulation 2, since this quantity is varied (it is the quantity on the horizontal axis in the plots in the Performance section of the user interface)

## Adaptive Design Per-stage Sample Sizes
per_stage_sample_size_combined_AD = 280, #(Range: 0 to 1000)
per_stage_sample_size_subpop_1_AD = 148, #(Range: 0 to 1000)

## Alpha allocation
# Desired familywise type I error rate (one-sided) (Range: 0 to 1)
FWER = 0.025,
# Proportion of FWER to test of H0C (Range: 0 to 1)
FWER_H0C_proportion = 0.09,

## End of list of User Controlled Parameters that are set with Sliders

## Additional Parameters (expected to be less frequently changed than parameters set by sliders) to be input using textboxes under PARAMETERS Tab

## Parameters used in all designs (adaptive and standard)
# Group Sequential Boundary Parameter (exponent) (Range: -0.5 to 0.5)
Delta = -0.5,
# Number simulated trials used to compute power, expected sample size, and expected trial duration for plots and tables
iter = 10000, # Range: 1 to 500,000
# Time limit for primary function 
time_limit = 45, #(seconds, range from 5 to 60)
# Number stages in trial design
total_number_stages = 5, # Range 1:20

## Parameters used only by adaptive design
last_stage_subpop_2_enrolled_AD = 3,  #(Range: 1 to total_number_stages)
# Stopping boundary proportionality constant for subpopulation 2 in adaptive design (z-statistic scale)

# Enrollment rate for combined population (patients per year)
enrollment_rate_combined_population = 420,

# Per stage sample size for standard design enrolling combined population
per_stage_sample_size_combined_SC = 106, #(Range: 0 to 1000)

# Per stage sample size for standard design enrolling on subpopulation 1 
per_stage_sample_size_SS = 100, #(Range: 0 to 1000)

subpop_2_stop_boundary_constant_AD = 0, #(Range: -10 to 10)

# Futility boundary proportionality constant for H01 for adaptive design (z-statistic scale)
H01_futility_boundary_constant_AD = 0, #(Range: -10 to 10)
## Parameters used only by standard designs
# Futility boundary proportionality constant for standard design enrolling combined population (z-statistic scale)
H0C_futility_boundary_constant_SC = -0.1, #(Range: -10 to 10)
# Futiltiy boundary proportionality constant for standard design enrolling subpopulation 1 only
H01_futility_boundary_constant_SS = -0.1, ##(Range: -10 to 10)

# Horizontal axis range in plots, in terms of mean treatment effect on risk difference scale:
lower_bound_treatment_effect_subpop_2 = -0.2,  # range (-1,1)
upper_bound_treatment_effect_subpop_2 = 0.2,  # range (-1,1)

CSV){



#### check for invalid inputs


#First check if all are in the range required by params.csv
all_default_inputs <- read.csv(paste0(path.package("interAdapt"),"/csv/params.csv"), header=TRUE, as.is=TRUE)
for(i in 1:dim(all_default_inputs)[1]){

	var_value_i<-eval(parse(text=all_default_inputs[i,'inputId']))
	#two warnings:
	if( var_value_i>all_default_inputs[i,'max'])  {
		assign(all_default_inputs[i,"inputId"],value=all_default_inputs[i,"max"])
		warning(paste0(all_default_inputs[i,"inputId"], ' is outside the allowed range, and has been set to ',all_default_inputs[i,"max"]))
	}
	if( var_value_i<all_default_inputs[i,'min'])  {
		assign(all_default_inputs[i,"inputId"],value=all_default_inputs[i,"min"])
		warning(paste0(all_default_inputs[i,"inputId"], ' is outside the allowed range, and has been set to ',all_default_inputs[i,"min"]))
	}
}
#Next check some variable specific requirements
if(total_number_stages<last_stage_subpop_2_enrolled_AD){
    last_stage_subpop_2_enrolled_AD<-total_number_stages
    warning(paste0("The last stage subpopulation 2 is enrolled must be less than the total number of stages. Here the last stage in which subpopulation 2 is enrolled has been set to ",total_number_stages,", the total number of stages"))
}



#### Import user's csv files
if (!missing(CSV)) {
    error <- 0
    if(inherits(CSV, "connection"))
	CSV <- list(CSV)
    for (file in CSV) {
	# parameter table format:
	#	inputId,label,min,max,value,step[,animate]
	#	  ***          ^   ^   ***
	# only inputId and value are required
	# range check done if min/max are given
	# label, step, and animate are ignored
	local <- is.character(file)
	if (local)
		cat("reading", file, "...")
	else
		print(file)
        t<-read.csv(file, header=TRUE, as.is=TRUE)
	if (local)
		cat(" Done\n")
	values <- t[,'value']
	Names <- names(t)
	names <- t[,'inputId']
	for (name in names)
		if (!exists(name))
			stop("not a parameter:", name)
	if ("min" %in% Names && "max" %in% Names) {
		m <- t[,'min']
		M <- t[,'max']
		if(!all(m <= values & values <= M)) {
			error <- 1
			cat("out of range:\n")
			for (i in 1:nrow(t))
				if(m[i] > values[i] || values[i] > M[i])
					cat(names[i], "\t", values[i], "\n")
		}
	}
        print(data.frame(value=values, row.names=names))
        for(i in 1:nrow(t))
                assign(names[i], values[i])
    }
    if(error)
	stop()
}
#### End Import user's csv files

## List of global variables not directly accessible by user; these are functions of above variables, but are updated every time table_constructor function is called
H0C_efficacy_boundary_proportionality_constant_standard_design <- 2.04
H01_efficacy_boundary_proportionality_constant_standard_design <- 2.04
H0C_efficacy_boundary_proportionality_constant_adaptive_design <- 2.54
H01_efficacy_boundary_proportionality_constant_adaptive_design <- 2.12
subpop_1_efficacy_boundaries_adaptive_design <- H01_efficacy_boundary_proportionality_constant_adaptive_design*((1:total_number_stages)/total_number_stages)^Delta
# List of treatment effect values (on risk difference scale) at which power, expected sample size, and expected duration will be evaluated 
risk_difference_list <- seq(lower_bound_treatment_effect_subpop_2,upper_bound_treatment_effect_subpop_2,length=7)

# Construct stopping boundaries for adaptive design based on user-input proportionality constants:
subpopulation_2_stopping_boundaries_adaptive_design <- c(subpop_2_stop_boundary_constant_AD*(((1:(last_stage_subpop_2_enrolled_AD-1))/(last_stage_subpop_2_enrolled_AD-1))^Delta),rep(Inf,total_number_stages-last_stage_subpop_2_enrolled_AD+1))
# Compute subpopulation 1 cumulative sample size vector for adaptive design
if(total_number_stages>last_stage_subpop_2_enrolled_AD){
	subpop_1_sample_size_vector <- c((1:last_stage_subpop_2_enrolled_AD)*per_stage_sample_size_combined_AD*p_1,(last_stage_subpop_2_enrolled_AD*per_stage_sample_size_combined_AD*p_1) + (1:(total_number_stages-last_stage_subpop_2_enrolled_AD))*per_stage_sample_size_subpop_1_AD)
} else {
	subpop_1_sample_size_vector <- c((1:last_stage_subpop_2_enrolled_AD)*per_stage_sample_size_combined_AD*p_1)
}
subpop_1_futility_boundaries_adaptive_design <- c(H01_futility_boundary_constant_AD*(subpop_1_sample_size_vector[1:total_number_stages-1]/subpop_1_sample_size_vector[total_number_stages-1])^Delta,H01_efficacy_boundary_proportionality_constant_adaptive_design)

# Construct stopping boundaries for standard designs based on user-input proportionality constants:
combined_pop_futility_boundaries_standard_design_H0C <-c(H0C_futility_boundary_constant_SC*(((1:(total_number_stages-1))/(total_number_stages-1))^Delta),Inf)
subpop_1_futility_boundaries_standard_design_H0C <-c(rep(Inf,total_number_stages-1),Inf)

combined_pop_futility_boundaries_standard_design_H01 <-c(H01_futility_boundary_constant_SS*(((1:(total_number_stages-1))/(total_number_stages-1))^Delta),Inf)
subpop_1_futility_boundaries_standard_design_H01 <-c(rep(Inf,total_number_stages-1),Inf)

# Function to compute power for given design and data generating distribution.
# The data generating distribution is defined by the signal to noise ratio (SNR) for subpopulation 1 patients: SNR_subpop_1
# and the signal to noise ratio (SNR) for subpopulation 2 patients: SNR_subpop_2. 
# SNR for a given population s (s=1 or 2) is defined as the mean difference (p_{st}-p_{sc}) divided by the corresponding standard deviation. Formulas for these quantities are given below.
# The inputs to the function are:
#  design_type: may be "adaptive" or "standard"
#  p1: proportion of population in subpopulation 1 (denoted pi_1 in the documentation)
#  total_number_stages: total number of stages in the trial (denoted K in the documentation)
#  standard_design_futility_boundaries: only used for standard designs
#  subpop_1_futility_boundaries: vector (l_{1,1},...l_{1,K}), where each l_{1,k} is the subpopulation 1 futility boundary defined in the documentation
#  subpop_2_futility_boundaries: vector (l_{2,1},...l_{2,K}), where each l_{2,k} is the subpopulation 2 futility boundary defined in the documentation
#  n1: per-stage sample size for combined population during each stage at or before last_stage_subpop_2_enrolled_AD (i.e., for each stage k <= k*)
#  n2: per-stage sample size for subpopulation 1 during each stage from last_stage_subpop_2_enrolled_AD through last stage (i.e., for each stage k > k*)
#  e_AD_C: the proportionality constant for efficacy boundaries for the combined population in the adaptive design, denoted e_{AD,C} in documentation
#  e_AD_1: the proportionality constant for efficacy boundaries for subpopulation 1 in the adaptive design, denoted e_{AD,1} in documentation
#  In the following, r1 and r2 are the probabilities of being randomized to treatment, for subpopulation 1 and subpopulation 2, respectively; these are set throughout to equal 1/2, so as to achieve 1:1 randomization to treatment vs. control. 
#  SNR_subpop_1: (p11-p10)/sqrt(p11*(1-p11)/r1+p10*(1-p10)/(1-r1)) 
#  SNR_subpop_2:	(p21-p20)/sqrt(p21*(1-p21)/r2+p20*(1-p20)/(1-r2)) 
#  outcome_variance_subpop_1: p11*(1-p11)/r1+p10*(1-p10)/(1-r1)
#  outcome_variance_subpop_2: p21*(1-p21)/r2+p20*(1-p20)/(1-r2)

get_power <- function(design_type="adaptive",p1,total_number_stages=5,k,standard_design_futility_boundaries=rep(-Inf,total_number_stages),subpop_1_futility_boundaries=rep(-Inf,total_number_stages),subpop_2_futility_boundaries=rep(-Inf,total_number_stages),n1,n2,e_AD_C,e_AD_1,SNR_subpop_1,SNR_subpop_2,outcome_variance_subpop_1,outcome_variance_subpop_2){
	p2 <- (1-p1) # proportion of population in subpopulation 2
	# Enrollment rate subpop. 1 (patients per year)
	enrollment_rate_subpop_1 <- p_1*enrollment_rate_combined_population
	# Enrollment rate subpop. 2 (patients per year)
	enrollment_rate_subpop_2 <- (1-p_1)*enrollment_rate_combined_population 

	# These are sample sizes of number who have final outcomes available at each interim analysis
	
	combined_population_stagewise_sample_sizes <- c(rep(n1,k),rep(n2,total_number_stages-k))
	subpop_1_stagewise_sample_sizes <- c(rep(p1*n1,k),rep(n2,total_number_stages-k))
	subpop_2_stagewise_sample_sizes <- combined_population_stagewise_sample_sizes - subpop_1_stagewise_sample_sizes

	# Build up vectors of cumulative sample sizes for combined population and each subpopulation, respectively, based on stagewise sample sizes 
	cum_sample_sizes_combined_population <- rep(0,total_number_stages)
	cum_sample_sizes_subpop_1 <- rep(0,total_number_stages)
	cum_sample_sizes_subpop_2 <- rep(0,total_number_stages)
	for(stage in 1:total_number_stages)
	{
		if(stage==1){cum_sample_sizes_subpop_1[stage] <- subpop_1_stagewise_sample_sizes[1]} else {cum_sample_sizes_subpop_1[stage] <- cum_sample_sizes_subpop_1[stage-1] + subpop_1_stagewise_sample_sizes[stage]}
		if(stage==1){cum_sample_sizes_subpop_2[stage] <- subpop_2_stagewise_sample_sizes[1]} else {cum_sample_sizes_subpop_2[stage] <- cum_sample_sizes_subpop_2[stage-1] + subpop_2_stagewise_sample_sizes[stage]}
	}
	cum_sample_sizes_combined_population <- cum_sample_sizes_subpop_2 + cum_sample_sizes_subpop_1
 
	combined_population_efficacy_boundaries <- e_AD_C*c((cum_sample_sizes_combined_population[1:k]/cum_sample_sizes_combined_population[k])^Delta,rep(Inf,total_number_stages-k))
	subpop_1_efficacy_boundaries <- e_AD_1*(cum_sample_sizes_subpop_1/cum_sample_sizes_subpop_1[total_number_stages])^Delta
        
## Get list of sample sizes corresponding to each interim analysis
	all_relevant_subpop_1_sample_sizes <- sort(unique(c(cum_sample_sizes_subpop_1)))
	all_relevant_subpop_2_sample_sizes <- sort(unique(c(cum_sample_sizes_subpop_2)))

## generate z-statistic increments
	Z_subpop_1_increment <- array(0,c(length(all_relevant_subpop_1_sample_sizes),iter)) 
	Z_subpop_1_increment[1,] <- rnorm(iter)+SNR_subpop_1*sqrt(all_relevant_subpop_1_sample_sizes[1])
	if(length(all_relevant_subpop_1_sample_sizes)>1)
	{	for(i in 2:length(all_relevant_subpop_1_sample_sizes))
		{
			Z_subpop_1_increment[i,] <- rnorm(iter)+SNR_subpop_1*sqrt(all_relevant_subpop_1_sample_sizes[i]-all_relevant_subpop_1_sample_sizes[i-1])
		}
	}
	Z_subpop_2_increment <- array(0,c(length(all_relevant_subpop_2_sample_sizes),iter)) 
	Z_subpop_2_increment[1,] <- rnorm(iter)+SNR_subpop_2*sqrt(all_relevant_subpop_2_sample_sizes[1])
	if(length(all_relevant_subpop_2_sample_sizes)>1)
	{
		for(i in 2:length(all_relevant_subpop_2_sample_sizes))
		{
			Z_subpop_2_increment[i,] <- rnorm(iter)+SNR_subpop_2*sqrt(all_relevant_subpop_2_sample_sizes[i]-all_relevant_subpop_2_sample_sizes[i-1])
		}
	}
	
## generate partial sums of increments
## Construct cumulative z-statistics:
	# First for subpop_1 
	Z_subpop_1_partial_weighted_sum_of_increments <- Z_subpop_1_increment
	if(length(all_relevant_subpop_1_sample_sizes)>1)
	{
		for(i in 2:length(all_relevant_subpop_1_sample_sizes))
		{
			Z_subpop_1_partial_weighted_sum_of_increments[i,] <- 
		((sqrt(all_relevant_subpop_1_sample_sizes[i-1]/all_relevant_subpop_1_sample_sizes[i])*Z_subpop_1_partial_weighted_sum_of_increments[i-1,])		
			+ (sqrt((all_relevant_subpop_1_sample_sizes[i]-all_relevant_subpop_1_sample_sizes[i-1])/all_relevant_subpop_1_sample_sizes[i])*Z_subpop_1_increment[i,]))
		}
	}
	Z_subpop_1_cumulative <- array(0,c(total_number_stages,iter))
	for(i in 1:total_number_stages){
		index <- which(all_relevant_subpop_1_sample_sizes==cum_sample_sizes_subpop_1[i])
		Z_subpop_1_cumulative[i,] <- Z_subpop_1_partial_weighted_sum_of_increments[index,]
	}
	# For subpopulation 2
	Z_subpop_2_partial_weighted_sum_of_increments <- Z_subpop_2_increment
	if(length(all_relevant_subpop_2_sample_sizes)>1)
	{
		for(i in 2:length(all_relevant_subpop_2_sample_sizes))
		{
			Z_subpop_2_partial_weighted_sum_of_increments[i,] <- 
		((sqrt(all_relevant_subpop_2_sample_sizes[i-1]/all_relevant_subpop_2_sample_sizes[i])*Z_subpop_2_partial_weighted_sum_of_increments[i-1,])		
			+ (sqrt((all_relevant_subpop_2_sample_sizes[i]-all_relevant_subpop_2_sample_sizes[i-1])/all_relevant_subpop_2_sample_sizes[i])*Z_subpop_2_increment[i,]))
		}
	}
	Z_subpop_2_cumulative <- array(0,c(total_number_stages,iter))
	for(i in 1:total_number_stages){
		index <- which(all_relevant_subpop_2_sample_sizes==cum_sample_sizes_subpop_2[i])
		Z_subpop_2_cumulative[i,] <- Z_subpop_2_partial_weighted_sum_of_increments[index,]
	}
	# Define combined_population population z-statistics
	variance_component1 <- (p1^2)*outcome_variance_subpop_1/cum_sample_sizes_subpop_1
	if(p2!=0){variance_component2 <- (p2^2)*outcome_variance_subpop_2/cum_sample_sizes_subpop_2}else{variance_component2 <- 0*variance_component1}
	correlation_Z_subpop_1_with_Z_combined_population <- sqrt(variance_component1/(variance_component1+variance_component2)) 
	correlation_Z_subpop_2_with_Z_combined_population <- sqrt(variance_component2/(variance_component1+variance_component2))
	Z_combined_population_cumulative <- (correlation_Z_subpop_1_with_Z_combined_population*Z_subpop_1_cumulative + correlation_Z_subpop_2_with_Z_combined_population*Z_subpop_2_cumulative)
	
## Determine outcomes of each simulated trial
if(design_type=="adaptive"){
    # record if efficacy boundary ever crossed, for each of H0C and H01:
 	ever_cross_H0C_efficacy_boundary <- rep(0,iter)
	ever_cross_H01_efficacy_boundary <- rep(0,iter)
	# indicator of stopping all enrollment, and of stopping only subpopulation 2, respectively:
	all_stopped <- rep(0,iter)
    subpop_2_stopped <- rep(0,iter)
    # indicators of rejecting null hypotheses:
	reject_H0C <- rep(0,iter)
    reject_H01 <- rep(0,iter)
    # record stage (just) after which trial stops
	final_stage_subpop_1_enrolled_up_through <- rep(total_number_stages,iter)
    final_stage_subpop_2_enrolled_up_through <- rep(total_number_stages,iter)
	for(stage in 1:total_number_stages)
	{
		  #below, k represents k^* from paper
          if(stage <= k){ever_cross_H0C_efficacy_boundary <- ifelse(Z_combined_population_cumulative[stage,]>combined_population_efficacy_boundaries[stage],1,ever_cross_H0C_efficacy_boundary)} # since always stop H0C testing after stage k
          ever_cross_H01_efficacy_boundary <- ifelse(Z_subpop_1_cumulative[stage,]>subpop_1_efficacy_boundaries[stage],1,ever_cross_H01_efficacy_boundary)
		# Step 1 of algorithm: Determine if any new events where H0C rejected for efficacy:
         if(stage <= k){reject_H0C <- ifelse((!all_stopped) & (!subpop_2_stopped) & Z_combined_population_cumulative[stage,]>combined_population_efficacy_boundaries[stage],1,reject_H0C)}

         reject_H01 <- ifelse((!all_stopped) & Z_subpop_1_cumulative[stage,]>subpop_1_efficacy_boundaries[stage],1,reject_H01)
          
         all_stopped <- ifelse(reject_H0C | reject_H01 | (Z_subpop_1_cumulative[stage,]< subpop_1_futility_boundaries[stage]),1,all_stopped)          

         subpop_2_stopped <- ifelse(all_stopped | (Z_subpop_2_cumulative[stage,] < subpop_2_futility_boundaries[stage]),1,subpop_2_stopped)

        # force subpop 2 stop at stage k if not yet stopped already
        if(stage==k){subpop_2_stopped <- rep(1,iter)}
		
		# record at what stage each subpop. stopped
        final_stage_subpop_1_enrolled_up_through <- ifelse((final_stage_subpop_1_enrolled_up_through==total_number_stages) & (all_stopped==1),stage,final_stage_subpop_1_enrolled_up_through)
        final_stage_subpop_2_enrolled_up_through <- ifelse((final_stage_subpop_2_enrolled_up_through==total_number_stages) & (subpop_2_stopped==1),stage,final_stage_subpop_2_enrolled_up_through)
	}
return(c(
mean(cum_sample_sizes_subpop_1[final_stage_subpop_1_enrolled_up_through]+cum_sample_sizes_subpop_2[final_stage_subpop_2_enrolled_up_through]), # expected sample size
mean(cum_sample_sizes_subpop_1[final_stage_subpop_1_enrolled_up_through])/(p1*enrollment_rate_combined_population), # expected duration
mean(reject_H0C), # power to reject H0C
mean(reject_H01), # power to reject H01
mean(reject_H0C | reject_H01) # power to reject H01 or H0C
))
} else if(design_type=="standard"){  
    # record if efficacy boundary ever crossed:
 	ever_cross_H0C_efficacy_boundary <- rep(0,iter)
	# indicator of stopping all enrollment, and of stopping only subpopulation 2, respectively:
	all_stopped <- rep(0,iter)
    # indicator of rejecting null hypotheses:
	reject_H0C <- rep(0,iter)
    # record stage (just) after which trial stops
	final_stage_enrolled_up_through <- rep(total_number_stages,iter)
	for(stage in 1:total_number_stages)
	{
		  #below, k represents k^* from paper
          ever_cross_H0C_efficacy_boundary <- ifelse(Z_combined_population_cumulative[stage,]>combined_population_efficacy_boundaries[stage],1,ever_cross_H0C_efficacy_boundary)
		# Step 1 of algorithm: Determine if any new events where H0C rejected for efficacy:
         reject_H0C <- ifelse((!all_stopped) & Z_combined_population_cumulative[stage,]>combined_population_efficacy_boundaries[stage],1,reject_H0C)          
         all_stopped <- ifelse(reject_H0C | (Z_combined_population_cumulative[stage,]< standard_design_futility_boundaries[stage]),1,all_stopped)          		
		# record at what stage each subpop. stopped
        final_stage_enrolled_up_through <- ifelse((final_stage_enrolled_up_through==total_number_stages) & (all_stopped==1),stage,final_stage_enrolled_up_through)
	}
return(c(
mean(cum_sample_sizes_subpop_1[final_stage_enrolled_up_through]+cum_sample_sizes_subpop_2[final_stage_enrolled_up_through]), # expected sample size
mean(reject_H0C) # power to reject H0C
))
}}

#Binary search to find smallest proportionality constants e_AD_C and e_AD_1 such that worst-case familywise Type I error is at most alpha (set by user) for adaptive design
get_adaptive_efficacy_boundaries <- function(alpha_FWER,alpha_H0C,outcome_variance_subpop_1,
outcome_variance_subpop_2)
{
	p1 <- p_1
	p2 <- 1-p1
	k_star <- last_stage_subpop_2_enrolled_AD
ss<- rep(per_stage_sample_size_combined_AD,k_star)
if(k_star>1)
{
	for(i in 2:k_star){
	ss[i] <- ss[i-1]+per_stage_sample_size_combined_AD
	}
}
	cov_matrix <- diag(k_star)
for(i in 1:k_star){for(j in 1:k_star) cov_matrix[i,j] <- sqrt(min(ss[i],ss[j])/max(ss[i],ss[j]))}

boundary_vector_with_unit_proportionality_constant <- ((1:k_star)/k_star)^Delta

OF_prop_constant_upper_bnd <- 10
OF_prop_constant_lower_bnd <- 0
while(OF_prop_constant_upper_bnd-OF_prop_constant_lower_bnd > 0.000001)
{
	OF_prop_constant_midpt <- mean(c(OF_prop_constant_lower_bnd,OF_prop_constant_upper_bnd))
	type_I_error <- 1-(pmvnorm(lower=rep(-Inf,k_star),upper=OF_prop_constant_midpt*boundary_vector_with_unit_proportionality_constant,mean=rep(0,k_star),sigma=cov_matrix))
	if(type_I_error < alpha_H0C) OF_prop_constant_upper_bnd <- OF_prop_constant_midpt else OF_prop_constant_lower_bnd <- OF_prop_constant_midpt
}
H0C_prop_const <- OF_prop_constant_midpt
OF_final_boundaries_H0C <- OF_prop_constant_midpt*boundary_vector_with_unit_proportionality_constant

# Generate candidate efficacy boundaries for subpopulation 1
# First construct cumulative sample sizes for subpopulation 1
cum_sample_sizes_subpop_1 <- c(p_1*ss,rep(0,total_number_stages-k_star))
if(total_number_stages-k_star >0)
{
	for(i in (k_star+1):total_number_stages){
	cum_sample_sizes_subpop_1[i] <-  cum_sample_sizes_subpop_1[i-1]+ per_stage_sample_size_subpop_1_AD
	}
}
# Construct covariance matrix of Z_{C,1},dots,Z_{C,k_star},Z_{S,1},dots,Z_{S,total_number_stages}
cov_matrix_full <- cov_matrix_full_check <- diag(k_star+total_number_stages)
# upper left set to cov_matrix for Z_{C,j}:
cov_matrix_full[1:k_star,1:k_star] <- cov_matrix
# build lower right for Z_{S,j} 
for(i in 1:total_number_stages){for(j in 1:total_number_stages) cov_matrix_full[i+k_star,j+k_star] <- sqrt(min(cum_sample_sizes_subpop_1[i],cum_sample_sizes_subpop_1[j])/max(cum_sample_sizes_subpop_1[i],cum_sample_sizes_subpop_1[j]))}

cum_sample_sizes_subpop_2 <- c(ss*(1-p_1),rep(ss[k_star]*(1-p_1),total_number_stages-k_star))

# build upper right and lower left parts of covariance matrix Z_{C,i}Z_{S,j}
for(i in 1:k_star){
	for(j in 1:total_number_stages){
	cov_matrix_full[j+k_star,i] <-cov_matrix_full[i,j+k_star] <- sqrt((min(cum_sample_sizes_subpop_1[i],cum_sample_sizes_subpop_1[j])/max(cum_sample_sizes_subpop_1[i],cum_sample_sizes_subpop_1[j]))*(p1*outcome_variance_subpop_1/(p1*outcome_variance_subpop_1+p2*outcome_variance_subpop_2)))}}
# Do binary search to set efficacy threshold for subpopulation 1.
boundary_vector_with_unit_proportionality_constant_subpop_1 <- (cum_sample_sizes_subpop_1/cum_sample_sizes_subpop_1[total_number_stages])^Delta
#((1:total_number_stages)/total_number_stages)^Delta

OF_prop_constant_upper_bnd <- 10
OF_prop_constant_lower_bnd <- 0
while(OF_prop_constant_upper_bnd-OF_prop_constant_lower_bnd > 0.000001)
{
	OF_prop_constant_midpt <- mean(c(OF_prop_constant_lower_bnd,OF_prop_constant_upper_bnd))
	type_I_error <- 1-(pmvnorm(lower=rep(-Inf,k_star+total_number_stages),upper=c(OF_final_boundaries_H0C,OF_prop_constant_midpt*boundary_vector_with_unit_proportionality_constant_subpop_1),mean=rep(0,k_star+total_number_stages),sigma=cov_matrix_full))
	
	if(type_I_error < alpha_FWER) OF_prop_constant_upper_bnd <- OF_prop_constant_midpt else OF_prop_constant_lower_bnd <- OF_prop_constant_midpt
}
H01_prop_const <- OF_prop_constant_midpt
OF_final_boundaries_H01 <- OF_prop_constant_midpt*boundary_vector_with_unit_proportionality_constant_subpop_1
return(c(H0C_prop_const,H01_prop_const,boundary_vector_with_unit_proportionality_constant_subpop_1))
}

##
## Construct table of values to be displayed in user interface
##
table_constructor <- function(){

setTimeLimit(time_limit) # stops computation if taking greater than time_limit

error_counter <- 0 # checks if errors encountered
k <- last_stage_subpop_2_enrolled_AD
p1 <- p_1
p2 <- (1-p1)
p11 <- p_11
p10 <- p_10
p20 <- p_20

# Probability randomized to control Arm
# for Subpop. 1 (Range: 0 to 1)
r1 <- 1/2 
# for Subpop. 2 (Range: 0 to 1)
r2 <- 1/2

futility_boundaries_standard_design_H0C <-c(H0C_futility_boundary_constant_SC*(((1:(total_number_stages-1))/(total_number_stages-1))^Delta),Inf)

futility_boundaries_standard_design_H01 <-c(H01_futility_boundary_constant_SS*(((1:(total_number_stages-1))/(total_number_stages-1))^Delta),Inf)

#Placeholders: these are not used in algorithm
subpop_1_futility_boundaries_standard_design_H0C <- rep(-Inf,total_number_stages)
subpop_1_futility_boundaries_standard_design_H01 <- rep(-Inf,total_number_stages)
subpop_2_futility_boundaries_standard_design_H0C <- rep(-Inf,total_number_stages)
subpop_2_futility_boundaries_standard_design_H01 <- rep(-Inf,total_number_stages)
combined_population_stopping_boundaries_adaptive_design <- rep(-Inf,total_number_stages)

outcome_variance_subpop_1_under_null <- p10*(1-p10)/r1+p10*(1-p10)/(1-r1)
outcome_variance_subpop_2_under_null <- p20*(1-p20)/r2+p20*(1-p20)/(1-r2)
prop_consts <- get_adaptive_efficacy_boundaries(FWER,FWER_H0C_proportion*FWER,outcome_variance_subpop_1_under_null,outcome_variance_subpop_2_under_null)
H0C_efficacy_boundary_proportionality_constant_adaptive_design <<- prop_consts[1] 
H01_efficacy_boundary_proportionality_constant_adaptive_design <<- prop_consts[2] 
adaptive_design_boundary_vector_with_unit_proportionality_constant_subpop_1 <- prop_consts[3:(3+total_number_stages-1)]

## Compute efficacy boundary for standard design enrolling combined population and testing only H0C, such that worst-case Type I error is at most alpha
ss <- 1:total_number_stages
cov_matrix <- diag(total_number_stages)
for(i in 1:total_number_stages){for(j in 1:total_number_stages) cov_matrix[i,j] <- sqrt(min(ss[i],ss[j])/max(ss[i],ss[j]))}
boundary_vector_with_unit_proportionality_constant <- ((1:total_number_stages)/total_number_stages)^Delta
OF_prop_constant_upper_bnd <- 10
OF_prop_constant_lower_bnd <- 0
while(OF_prop_constant_upper_bnd-OF_prop_constant_lower_bnd > 0.000001)
{
	OF_prop_constant_midpt <- mean(c(OF_prop_constant_lower_bnd,OF_prop_constant_upper_bnd))
	type_I_error <- 1-(pmvnorm(lower=rep(-Inf,total_number_stages),upper=c(OF_prop_constant_midpt*boundary_vector_with_unit_proportionality_constant),mean=rep(0,total_number_stages),sigma=cov_matrix))

	if(type_I_error < FWER) OF_prop_constant_upper_bnd <- OF_prop_constant_midpt else OF_prop_constant_lower_bnd <- OF_prop_constant_midpt
}
H0C_efficacy_boundary_proportionality_constant_standard_design <<- OF_prop_constant_midpt
H01_efficacy_boundary_proportionality_constant_standard_design <<- OF_prop_constant_midpt

H0C_efficacy_boundaries <- H0C_efficacy_boundary_proportionality_constant_adaptive_design*c(((1:k)/k)^Delta,rep(Inf,total_number_stages-k))

subpop_1_efficacy_boundaries_adaptive_design <<- H01_efficacy_boundary_proportionality_constant_adaptive_design*adaptive_design_boundary_vector_with_unit_proportionality_constant_subpop_1

subpopulation_2_stopping_boundaries_adaptive_design <<- c(subpop_2_stop_boundary_constant_AD*(((1:(last_stage_subpop_2_enrolled_AD-1))/(last_stage_subpop_2_enrolled_AD-1))^Delta),rep(Inf,total_number_stages-last_stage_subpop_2_enrolled_AD+1))
# Compute subpop. 1 cumulative sample size vector
if(total_number_stages>last_stage_subpop_2_enrolled_AD){
	subpop_1_sample_size_vector <- c((1:last_stage_subpop_2_enrolled_AD)*per_stage_sample_size_combined_AD*p_1,(last_stage_subpop_2_enrolled_AD*per_stage_sample_size_combined_AD*p_1) + (1:(total_number_stages-last_stage_subpop_2_enrolled_AD))*per_stage_sample_size_subpop_1_AD)
} else {
	subpop_1_sample_size_vector <- c((1:last_stage_subpop_2_enrolled_AD)*per_stage_sample_size_combined_AD*p_1)
}
subpop_1_futility_boundaries_adaptive_design <<- c(H01_futility_boundary_constant_AD*(subpop_1_sample_size_vector[1:total_number_stages-1]/subpop_1_sample_size_vector[total_number_stages-1])^Delta,H01_efficacy_boundary_proportionality_constant_adaptive_design)

print_ss_and_boundaries_flag <- 1
SNR_subpop_1 <- (p11-p10)/sqrt(p11*(1-p11)/r1+p10*(1-p10)/(1-r1))
outcome_variance_subpop_1 <- p11*(1-p11)/r1+p10*(1-p10)/(1-r1)

risk_difference_list <<- sort(unique(c(seq(max(c(min(c(lower_bound_treatment_effect_subpop_2,upper_bound_treatment_effect_subpop_2,0)),-p20)),min(c(max(c(lower_bound_treatment_effect_subpop_2,upper_bound_treatment_effect_subpop_2,0)),1-p20)),length=10))))

standard_combined_population_df <- array(0,c(length(risk_difference_list),4))
standard_subpop_1_only_df <- array(0,c(length(risk_difference_list),3))
adaptive_df <- array(0,c(length(risk_difference_list),5))
overrun_df <- array(0,c(length(risk_difference_list),3))
counter_combined_population <- 1
counter_subpop_1 <- 1
counter_adaptive <- 1
 
for(percent_benefit_subpop_2 in rev(risk_difference_list))
{
	p21 <- p20 + percent_benefit_subpop_2 
	#next two lines ensure p21 is in the range  (.001 , .999)
	p21<-max(p21,.001)
	p21<-min(p21,.999)
	SNR_subpop_2 <- (p21-p20)/sqrt(p21*(1-p21)/r2+p20*(1-p20)/(1-r2)) 
    outcome_variance_subpop_2 <- p21*(1-p21)/r2+p20*(1-p20)/(1-r2)

power_vec <- get_power(design_type="adaptive",p1=p_1,total_number_stages,k=last_stage_subpop_2_enrolled_AD,combined_population_stopping_boundaries_adaptive_design,subpop_1_futility_boundaries_adaptive_design,subpopulation_2_stopping_boundaries_adaptive_design,n1=per_stage_sample_size_combined_AD,n2=per_stage_sample_size_subpop_1_AD,e_AD_C=H0C_efficacy_boundary_proportionality_constant_adaptive_design,e_AD_1=H01_efficacy_boundary_proportionality_constant_adaptive_design,SNR_subpop_1=SNR_subpop_1,SNR_subpop_2=SNR_subpop_2,outcome_variance_subpop_1,outcome_variance_subpop_2)

adaptive_df[counter_adaptive,] <-  power_vec # Expected sample size, Expected Duration, Power to Reject H0C, Power to reject H01, Power to reject at least one of H0C or H01
overrun_df[counter_adaptive,1] <- 0
counter_adaptive <- counter_adaptive + 1

power_vec <- get_power(design_type="standard",p1=p_1,total_number_stages=total_number_stages,k=total_number_stages,futility_boundaries_standard_design_H0C,subpop_1_futility_boundaries=subpop_1_futility_boundaries_standard_design_H0C,subpop_2_futility_boundaries=subpop_2_futility_boundaries_standard_design_H0C,n1=per_stage_sample_size_combined_SC,n2=0,e_AD_C=H0C_efficacy_boundary_proportionality_constant_standard_design,e_AD_1=Inf,SNR_subpop_1=SNR_subpop_1,SNR_subpop_2=SNR_subpop_2,outcome_variance_subpop_1,outcome_variance_subpop_2)
standard_combined_population_df[counter_combined_population,] <- c(power_vec[1],power_vec[1]/enrollment_rate_combined_population,power_vec[2],0)
overrun_df[counter_combined_population,2] <- 0
counter_combined_population <- counter_combined_population +1

power_vec <- get_power(design_type="standard",p1=1,total_number_stages=total_number_stages,k=total_number_stages,futility_boundaries_standard_design_H01,subpop_1_futility_boundaries=subpop_1_futility_boundaries_standard_design_H01,subpop_2_futility_boundaries=subpop_2_futility_boundaries_standard_design_H01,n1=per_stage_sample_size_SS,n2=0,e_AD_C=H01_efficacy_boundary_proportionality_constant_standard_design,e_AD_1=Inf,SNR_subpop_1=SNR_subpop_1,SNR_subpop_2=SNR_subpop_2,outcome_variance_subpop_1,outcome_variance_subpop_2)
standard_subpop_1_only_df[counter_subpop_1,] <- c(power_vec[1],power_vec[1]/(p1*enrollment_rate_combined_population),power_vec[2])
overrun_df[counter_subpop_1,3] <- 0

counter_subpop_1 <- counter_subpop_1 +1

}
return(data.frame(cbind(rev(risk_difference_list),adaptive_df,standard_combined_population_df,standard_subpop_1_only_df,overrun_df)))
}

## Construct Power Curve Plot for display in user interface
power_curve_plot <- function()
{
plot(0,type="n",xlim=c(min(risk_difference_list),max(risk_difference_list)),ylim=c(0,1),main="Power versus Average Treatment Effect in Subpopulation 2",xlab="Avg. Treatment Effect on Risk Difference Scale in Subpopulation 2",ylab="Power")

lines(x=rev(risk_difference_list),y=table1[,6],lty=1,col=1,lwd=3)
# H0C adaptive
lines(x=rev(risk_difference_list),y=table1[,4],lty=2,col=1,lwd=3)
# H01 adaptive
lines(x=rev(risk_difference_list),y=table1[,5],lty=3,col=1,lwd=3)

# H0C standard
lines(x=rev(risk_difference_list),y=table1[,9],lty=4,col=3,lwd=3)

# H01 standard
lines(x=rev(risk_difference_list),y=table1[,13],lty=5,col=4,lwd=3)
ltext<-rep(NA,5)
ltext[1]<-expression(paste("Adaptive, Power to Reject at Least one of H"[0][C]," or H"[0][1]))
ltext[2]<-expression(paste("Adaptive, Power to Reject H"[0][C]))
ltext[3]<-expression(paste("Adaptive, Power to Reject H"[0][1]))
ltext[4]<-expression(paste("Standard Design Total Pop., Power to Reject H"[0][C],""))
ltext[5]<-expression(paste("Standard Design Subpop. 1 Only, Power to Reject H"[0][1],""))
legend("bottomright",legend=ltext,lty=c(1,2,3,4,5),col=c(1,1,1,3,4),lwd=c(3,3,3,3,3))
}

## Expected Sample Size Plot
expected_sample_size_plot <- function()
{
min_ess <- 0
max_ess <- max(c(table1[,2],table1[,7],table1[,11]))
plot(0,type="n",xlim=c(min(risk_difference_list),max(risk_difference_list)),ylim=c(min_ess,max_ess),main="Expected Sample Size versus Average Treatment Effect in Subpopulation 2",xlab="Avg. Treatment Effect on Risk Difference Scale in Subpopulation 2",ylab="Expected Sample Size")

# adaptive
lines(x=rev(risk_difference_list),y=table1[,2],lty=1,col=1,lwd=3)

# H0C standard
lines(x=rev(risk_difference_list),y=table1[,7],lty=2,col=3,lwd=3)

# H01 standard
lines(x=rev(risk_difference_list),y=table1[,11],lty=3,col=4,lwd=3)
legend("bottomright",legend=c("Adaptive Design","Standard Design Total Pop.","Standard Design Subpop. 1 Only"),lty=c(1,2,3),col=c(1,3,4),lwd=c(3,3,3))

}

## Expected Duration Plot
expected_duration_plot <- function()
{
min_dur <- 0
max_dur <- max(c(table1[,3],table1[,8],table1[,12]))
plot(0,type="n",xlim=c(min(risk_difference_list),max(risk_difference_list)),ylim=c(min_dur,max_dur),main="Expected Duration versus Average Treatment Effect in Subpopulation 2",xlab="Avg. Treatment Effect on Risk Difference Scale in Subpopulation 2",ylab="Expected Duration in Years")

# adaptive
lines(x=rev(risk_difference_list),y=table1[,3],lty=1,col=1,lwd=3)

# H0C standard
lines(x=rev(risk_difference_list),y=table1[,8],lty=2,col=3,lwd=3)

# H01 standard
lines(x=rev(risk_difference_list),y=table1[,12],lty=3,col=4,lwd=3)
legend("bottomright",legend=c("Adaptive Design","Standard Design Total Pop.","Standard Design Subpop. 1 Only"),lty=c(1,2,3),col=c(1,3,4),lwd=c(3,3,3))

}

## Construct diagram of efficacy and futility boundaries for adaptive design
boundary_adapt_plot <-function()
{
adapt_boundary_mat<- t(adaptive_design_sample_sizes_and_boundaries_table()[[1]][c("H0C Efficacy Boundaries u(C,k) for z-statistics Z(C,k)", "Boundaries l(2,k) for Z(2,k) to Stop Subpop. 2 Enrollment", "H01 Efficacy Boundaries u(1,k) for Z(1,k)", "Boundaries l(1,k) for Z(1,k) to Stop All Enrollment"), ])
adapt_boundary_mat[adapt_boundary_mat==Inf]<-NA
fancyTitle<-expression(atop('Decision Boundaries for Sequential Test of Combined Population','Null Hypothesis ( H'[0][C]~') and Subpopulation 1 Null Hypothesis ( H'[0][1]~')'))
matplot(adapt_boundary_mat,type='o',main=fancyTitle,pch=c(0,1,2,3),col=c('red','red','blue','blue'),lty=2,cex=1.5, xlab='Stage',ylab='Boundaries on Z-score scale')
ltext<-rep(NA,4)
ltext[1]<-expression(paste("H"[0][C]," Efficacy Boundaries"))
ltext[2]<-expression(paste("Boundaries l(2,k) for Z(2,k) to Stop Subpop. 2"))
ltext[3]<-expression(paste('H'[0][1],' Efficacy Boundaries'))
ltext[4]<-expression(paste("Boundaries l(1,k) for Z(1,k) to Stop All Enrollment"))
legend('topright',ltext,pch=c(0,1,2,3),col=c('red','red','blue','blue'),lty=2,cex=1.25) 
}

## Construct diagram of efficacy and futility boundaries for design enrolling combined population
boundary_standard_H0C_plot <-function()
{
H0C_boundary_mat<- t(standard_H0C_design_sample_sizes_and_boundaries_table()[[1]][c("H0C Efficacy Boundaries for z-statistics Z(C,k)", "H0C Futility Boundaries for z-statistics Z(C,k)"), ])
fancyTitle<-expression(atop('Decision Boundaries for Sequential Test of', 'Combined Population Null Hypothesis ( H'[0][C]~')'))
matplot(H0C_boundary_mat,type='o',main=fancyTitle,lty=2,pch=c(0,1),col='red', xlab='Stage',ylab='Boundaries on Z-score scale',cex=1.5)
ltext<-rep(NA,2)
ltext[1]<-expression(paste('H'[0][C],' Efficacy Boundaries'))
ltext[2]<-expression(paste('H'[0][C],' Futility Boundaries'))
legend('topright',ltext,lty=2,pch=c(0,1),col='red',cex=1.25)
}

## Construct diagram of efficacy and futility boundaries for design enrolling only subpopulation 1
boundary_standard_H01_plot <-function()
{
H01_boundary_mat<- t(standard_H01_design_sample_sizes_and_boundaries_table()[[1]][c("H01 Efficacy Boundaries for z-statistics Z(1,k)", "H01 Futility Boundaries for z-statistics Z(1,k)"), ])
fancyTitle<-expression(atop('Decision Boundaries for Sequential Test of', 'Combined Population Null Hypothesis ( H'[0][1]~')'))
matplot(H01_boundary_mat,type='o',main=fancyTitle,lty=2,pch=c(0,1),col='blue', xlab='Stage',ylab='Boundaries on Z-score scale',cex=1.5)
ltext<-rep(NA,2)
ltext[1]<-expression(paste('H'[0][1],' Efficacy Boundaries'))
ltext[2]<-expression(paste('H'[0][1],' Futility Boundaries'))
legend('topright',ltext,lty=2,pch=c(0,1),col='blue',cex=1.25)
}

## Construct table displaying performance of each design, including the following: power, expected sample size, and expected trial duration
performance_table <- function(table1)
{
output_df <- cbind(as.matrix(table1))
output_df_formatted <- cbind(output_df[,1],output_df[,2],output_df[,3],100*output_df[,4],100*output_df[,5],100*output_df[,6],output_df[,7],output_df[,8],100*output_df[,9],
#100*output_df[,10],
output_df[,11],output_df[,12],100*output_df[,13])
colnames(output_df_formatted) <- c("Subpop.2 Tx. Effect","AD:Sample Size","AD:DUR","AD:Power H0C","AD:Power H01","AD:Power H0C or H01","SC:Sample Size","SC:DUR","SC:Power H0C","SS:Sample Size","SS:DUR","SS:Power H01")
return(list(output_df_formatted,digits=c(0,2,0,1,0,0,0,0,1,0,0,1,0),caption=paste0("Comparison of avg sample size, avg duration (DUR), and power (as a percent), for the following designs: the Adaptive Design (AD), the Standard Design Enrolling Combined Population (SC), and the Standard Design Enrolling Subpop. 1 Only (SS). All designs strongly control the familywise Type I error rate at ",FWER,".")))
}

## Format the performance table
transpose_performance_table<-function(ptable){
	#make a matrix of digit values for each cell
	#take off the 1st entry (=0) from the digits vector (for the head col names)
	#and add a zero column for the new row names
	dpt<-dim(ptable[[1]])
	digitMatPre<-matrix(ptable$digits[-1],nrow=dpt[2],ncol=dpt[1],byrow=FALSE)
	digitMat<-cbind(0,digitMatPre) #add a column of zeros for the row names
	orderedTab<-t(ptable[[1]])[,(dpt[1]:1)] #reverse the col order
	outTab<-data.frame(orderedTab)
	#We need to say include.colnames=FALSE in our custom renderTable function
	return(list( outTab, digits=digitMat, caption=ptable$caption ) )
}

## Construct table displaying efficacy and futility boundaries for adaptive design
adaptive_design_sample_sizes_and_boundaries_table <- function()
{
k <- last_stage_subpop_2_enrolled_AD
p1 <- p_1
p2 <- 1-p1

H0C_efficacy_boundaries <- H0C_efficacy_boundary_proportionality_constant_adaptive_design*c(((1:k)/k)^Delta,rep(NA,total_number_stages-k))

subpopulation_2_stopping_boundaries_adaptive_design <<- c(subpop_2_stop_boundary_constant_AD*(((1:(last_stage_subpop_2_enrolled_AD-1))/(last_stage_subpop_2_enrolled_AD-1))^Delta),Inf,rep(NA,total_number_stages-last_stage_subpop_2_enrolled_AD))
# Compute subpop. 1 cumulative sample size vector
if(total_number_stages>last_stage_subpop_2_enrolled_AD){
	subpop_1_sample_size_vector <- c((1:last_stage_subpop_2_enrolled_AD)*per_stage_sample_size_combined_AD*p_1,(last_stage_subpop_2_enrolled_AD*per_stage_sample_size_combined_AD*p_1) + (1:(total_number_stages-last_stage_subpop_2_enrolled_AD))*per_stage_sample_size_subpop_1_AD)
} else {
	subpop_1_sample_size_vector <- c((1:last_stage_subpop_2_enrolled_AD)*per_stage_sample_size_combined_AD*p_1)
}
subpop_1_futility_boundaries_adaptive_design <<- c(H01_futility_boundary_constant_AD*(subpop_1_sample_size_vector[1:total_number_stages-1]/subpop_1_sample_size_vector[total_number_stages-1])^Delta,subpop_1_efficacy_boundaries_adaptive_design[total_number_stages])

if(k<total_number_stages){
row1 <- c(p1*per_stage_sample_size_combined_AD*(1:k),p1*per_stage_sample_size_combined_AD*k+per_stage_sample_size_subpop_1_AD*(1:(total_number_stages-k)))

row2 <- c(p2*per_stage_sample_size_combined_AD*(1:k),rep(p2*per_stage_sample_size_combined_AD*k,total_number_stages-k))

row3 <- c(per_stage_sample_size_combined_AD*(1:k),per_stage_sample_size_combined_AD*k+per_stage_sample_size_subpop_1_AD*(1:(total_number_stages-k)))
}else{
row1 <- c(p1*per_stage_sample_size_combined_AD*(1:k))

row2 <- c(p2*per_stage_sample_size_combined_AD*(1:k))

row3 <- c(per_stage_sample_size_combined_AD*(1:k))
}

H0C_efficacy <-  H0C_efficacy_boundaries

subpopulation_2_stopping_boundaries_adaptive_design_copy <- subpopulation_2_stopping_boundaries_adaptive_design
H01_efficacy <- subpop_1_efficacy_boundaries_adaptive_design
H01_futility <- subpop_1_futility_boundaries_adaptive_design

output_df <- rbind(row1,row2,row3,H0C_efficacy,subpopulation_2_stopping_boundaries_adaptive_design_copy,H01_efficacy,H01_futility)
row.names(output_df) <- c("Cumulative Sample Size Subpop. 1","Cumulative Sample Size Subpop. 2","Cumulative Sample Size Combined Pop.","H0C Efficacy Boundaries u(C,k) for z-statistics Z(C,k)","Boundaries l(2,k) for Z(2,k) to Stop Subpop. 2 Enrollment","H01 Efficacy Boundaries u(1,k) for Z(1,k)","Boundaries l(1,k) for Z(1,k) to Stop All Enrollment")
colnames(output_df) <- 1:total_number_stages
dig_array <- array(0,c(7,(total_number_stages+1)))

dig_array[4:7,2:(total_number_stages+1)] <- ifelse(is.na(output_df[4:7,]),0,ifelse(output_df[4:7,]==0,0,2))
return(list(output_df,digits=dig_array,caption="Cumulative Sample Sizes and Decision Boundaries for Adaptive Design. Each column corresponds to a stage. All thresholds are given on the z-statistic scale."))

}

## Construct table displaying efficacy and futility boundaries for standard design enrolling combined population
standard_H0C_design_sample_sizes_and_boundaries_table <- function()
{
k<-total_number_stages
p1 <- p_1
p2 <- 1-p1

H0C_efficacy_boundaries <- H0C_efficacy_boundary_proportionality_constant_standard_design*c(((1:k)/k)^Delta)

H0C_futility <- futility_boundaries_standard_design_H0C <-c(H0C_futility_boundary_constant_SC*(((1:(total_number_stages-1))/(total_number_stages-1))^Delta),H0C_efficacy_boundaries[total_number_stages])

row1 <- c(p1*per_stage_sample_size_combined_SC*(1:k))

row2 <- c(p2*per_stage_sample_size_combined_SC*(1:k))

row3 <- c(per_stage_sample_size_combined_SC*(1:k))

H0C_efficacy <-  H0C_efficacy_boundaries

output_df <- rbind(row1,row2,row3,H0C_efficacy,H0C_futility)
row.names(output_df) <- c("Cumulative Sample Size Subpop. 1","Cumulative Sample Size Subpop. 2","Cumulative Sample Size Combined Pop.","H0C Efficacy Boundaries for z-statistics Z(C,k)","H0C Futility Boundaries for z-statistics Z(C,k)")
colnames(output_df) <- 1:total_number_stages
dig_array <- array(0,c(5,(total_number_stages+1)))

dig_array[4:5,2:(total_number_stages+1)] <- ifelse(is.na(output_df[4:5,]),0,ifelse(output_df[4:5,]==0,0,2))
return(list(output_df,digits=dig_array,caption="Cumulative Sample Sizes and Decision Boundaries for Standard Design Enrolling Combined Population. Each column corresponds to a stage. All thresholds are given on the z-statistic scale."))


}

## Construct table displaying efficacy and futility boundaries for standard design enrolling subpopulation 1 only
standard_H01_design_sample_sizes_and_boundaries_table <- function()
{
k<-total_number_stages
p1 <- p_1
p2 <- 1-p1

H01_efficacy_boundaries <- H01_efficacy_boundary_proportionality_constant_standard_design*c(((1:k)/k)^Delta)

H01_futility <- futility_boundaries_standard_design_H01 <-c(H01_futility_boundary_constant_SS*(((1:(total_number_stages-1))/(total_number_stages-1))^Delta),H01_efficacy_boundaries[total_number_stages])


row1 <- c(per_stage_sample_size_SS*(1:k))

H01_efficacy <-  H01_efficacy_boundaries

output_df <- rbind(row1,H01_efficacy,H01_futility)
row.names(output_df) <- c("Cumulative Sample Size","H01 Efficacy Boundaries for z-statistics Z(1,k)","H01 Futility Boundaries for z-statistics Z(1,k)")
colnames(output_df) <- 1:total_number_stages
dig_array <- array(0,c(3,(total_number_stages+1)))
#dig_array[4:7,] <- array(2,c(4,(total_number_stages+1)))
dig_array[2:3,2:(total_number_stages+1)] <- ifelse(is.na(output_df[2:3,]),0,ifelse(output_df[2:3,]==0,0,2))
return(list(output_df,digits=dig_array,caption="Cumulative Sample Sizes and Decision Boundaries for Standard Design enrolling only Subpopulation 1. Each column corresponds to a stage. All thresholds are given on the z-statistic scale."))
}

table1 <- table_constructor()

ptable <- performance_table(table1)

Performance_Comparison_Among_Designs <- transpose_performance_table(ptable)

Performance_Comparison_Among_Designs$digits <- NULL

Performance_Comparison_Among_Designs$caption <- "Comparison of avg
sample size, avg duration (DUR), and power (as a percent), for the
following designs: the Adaptive Design (AD), the Standard Design
Enrolling Combined Population (SC), and the Standard Design Enrolling
Subpop. 1 Only (SS). All designs strongly control the familywise Type
I error rate at level set by user."

AD_Design_Per_Stage_Sample_Sizes_and_Boundaries <-
adaptive_design_sample_sizes_and_boundaries_table()

AD_Design_Per_Stage_Sample_Sizes_and_Boundaries$digits <- NULL

SC_Design_Per_Stage_Sample_Sizes_and_Boundaries <-
standard_H0C_design_sample_sizes_and_boundaries_table()

SC_Design_Per_Stage_Sample_Sizes_and_Boundaries$digits <- NULL

SS_Design_Per_Stage_Sample_Sizes_and_Boundaries <-
standard_H01_design_sample_sizes_and_boundaries_table()

SS_Design_Per_Stage_Sample_Sizes_and_Boundaries$digits <- NULL

inputIds <- read.csv(paste0(path.package("interAdapt"),"/csv/params.csv"), header=TRUE, as.is=TRUE)[,"inputId"]
Input_Parameters <- mget(inputIds)
names(Input_Parameters) <- inputIds

return(list(performance_comparison=Performance_Comparison_Among_Designs,AD_design=AD_Design_Per_Stage_Sample_Sizes_and_Boundaries,SC_design=SC_Design_Per_Stage_Sample_Sizes_and_Boundaries,SS_design=SS_Design_Per_Stage_Sample_Sizes_and_Boundaries,input_parameters=Input_Parameters))

}