# R/feasibility_check.R In aaronjfisher/designOptim: Design optimizer

#### Documented in binsearchtolfeasibility_checkmax_power_feasiblemin_n_feasiblemin_n_multistage

# Functions for checking feasibility of a design,
# and searching for the closest design that is feasible.

###### Helper functions
######

# Interleave two vectors, v1 and v2.
# [INTERNAL FUNCTION]
# This code is based on code from Bogdan Romocea, at https://stat.ethz.ch/pipermail/r-help/2006-March/101023.html
# @param v1 one vector to interleave
# @param v2 a second vector to interleave
interleave <- function(v1,v2){
if(length(v1)==0) return(v2)
if(length(v2)==0) return(v1)

ord1 <- 2*(1:length(v1))-1
ord2 <- 2*(1:length(v2))
c(v1,v2)[order(c(ord1,ord2))]
}

# Reorder a vector so that mid points come first. This is similar in spirit to binary search.
# This recursive function is called exactly as many times as the length of \code{x}.
# [INTERNAL FUNCTION]
# @param x is a vector of integers returned by order(x)
orderByMids<-function(x){

ox <- order(x)

move <- ceiling(length(x)/2)
xm<- x[move]
x1<- x[ox<move]
x2<- x[ox>move]

if(length(x1) ==0 & length(x2)==0) return(xm)

if(length(x1)>0) x1<-orderByMids(x1)
if(length(x2)>0) x2<-orderByMids(x2)

return(c(xm,interleave(x1,x2)))

}

######
######

#' Check feasibility of power constraints
#'
#' Functions to check feasibility of a design, adjust the design in order to make it feasible, or adjust the constraints so that they can be achieved. Here, "feasibility" means that there exists a 1-stage design at the supplied sample size that meets the supplied power and Type I error constraints.
#'
#' \code{feasibility_check} checks whether it is possible to meet a set of power constraints while maintaining a given familywise Type I error rate (FWER) and total sample size.
#'
#' \code{min_n_feasible} implements a binary search to find the smallest sample size that meets the power constraints and FWER specified. In each iteration of the search, \code{min_n_feasible} calls \code{feasibility_check}.
#'
#' Given a sample size and set of power cases, \code{max_power_feasible} finds the number \eqn{m} between 0 and 1 such that if all of the minimum power thresholds are multiplied by \eqn{m}, then these power constraints will be satisfied at the supplied value of \code{n_total}. In other words, it finds the factor by which the power constraints must be relaxed.
#'
#' @param FWER the required familywise Type I error rate for the trial
#' @param p1 population proportion in subpopulation 1
#' @param trial_method the type of trial to run. 'cov' for covariance based, 'MB' for Maurer Bretz, and 'covMB' for a combination approach.
#' @param r1 probability of being randomized to treatment in subpopulation 1
#' @param r2 probability of being randomized to treatment in subpopulation 2
#' @param n_total the total sample size for the 1-stage trial
#' @param cases a list of power constraints, of the same format as those passed to \code{\link{optimizeTrial}}
#' @param npoints_sqrt \code{feasibility_check} determines feasibility by searching over a grid of points for the alpha to be allocated between the two trials. Determines the number of grid points to search over. The search will be conducted over a triangle of points, with \code{(npoints_sqrt^2)/2} points.
#' @param trial_args either an empty list (for trial_method=='cov'), or a list containing graph edges for alpha reallocation (for trial_method=='MB').
#' @return \code{feasibility_check} returns a named vector, with the first element being 1 if the trial is feasible, and zero otherwise. The remaining elements of this vector tell the alpha allocated to each hypothesis in one particular setup that results in a feasible trial. \cr \code{min_n_feasible} returns the smallest sample size that meets the power constraints and FWER specified, as well as the output from \code{feasibility_check} at that sample size. \cr \code{max_power_feasible} returns the multiplier \eqn{m} and a list of modified cases in which each power threshold has been multiplied by \eqn{m}.
#' @export
#'
feasibility_check <- function(FWER, p1, trial_method, r1=0.5, r2=1-r1, n_total,cases,npoints_sqrt=10,trial_args=list()){

seq1 <- seq(10^-10,FWER,length=npoints_sqrt)
seq1 <- seq1[orderByMids(order(seq1))]
#reordering the grid search in a smart way can let us stop the search sooner if we find a feasible design.
if(npoints_sqrt==1) seq1 <- FWER / 3 #equal allocation to each hypothesis

for(alpha_1 in seq1){
seq2<- seq1[seq1<=FWER-alpha_1]

for(alpha_2 in seq2){

# Reset all_constraints_satisfied to TRUE
all_constraints_satisfied <- TRUE

for(case in cases){
# Set means and variances according to user input for power constraint:
#!! In future we can make this code cleaner by combining many of these into trial_args
# (specifically p1, r1, r2, FWER), but for now we're keeping it just to avoid
# compatibility errors.

trial_args$iter <- 10000 #!! Hard-coded trial_args$FWER <- FWER
trial_args$num_stages<-1 trial_args$n_total<-n_total
trial_args$n_per_stage<-c(n_total) trial_args$p1<-p1
trial_args$r1<-r1 trial_args$r2<-r2
trial_args$H01_eff_allocated<-alpha_1 trial_args$H02_eff_allocated<-alpha_2
trial_args$H0C_eff_allocated<-(FWER-(alpha_1+alpha_2)) trial_args$H01_futility_boundaries<- -Inf # these are essentially ignored anyway
trial_args$H02_futility_boundaries<- -Inf trial_args$H0C_futility_boundaries<- -Inf
trial_args$time_limit<-200 trial_args$enrollment_rate_combined<-100
trial_args$delay<-0#irrelevant since it's 1 stage trial_args$errtol<-0.01
trial_args$maxpts<-10000 trial_args$abseps<-.01

full_case_args <- c(trial_args, case)

if(trial_method=='cov') buildFun <- buildTrial
if(trial_method=='MB') buildFun <- buildTrial_Maurer_Bretz_2013

case$performance <- do.call(buildFun,full_case_args)$performance

# check that power constraints are satisfied
# if not, stop searching at this alpha configuration, and declare that we've failed to satisfy the requirements.
# If no alpha configurations succeed, this will be the last time we change all_constraints_satisfied
case_constraints_satisfied <- check_power_case(case)
if(!case_constraints_satisfied) {
all_constraints_satisfied <- FALSE
break #no need to check the remaining cases
}
}

#After we've searched over all stages, if we've met all constraints, stop the function.
if(all_constraints_satisfied){
return(c(
'feasible'=1,
'H01_eff_allocated'=alpha_1,
'H02_eff_allocated'=alpha_2,
'H0C_eff_allocated'=(FWER-alpha_1-alpha_2)
))
}
#print(trial_output$performance) } } #After we've searched over all configurations with no success, stop the trial. return(c('feasible'=0)) } # internal function check_power_case <-function(case){ power_ind <- c('Pow_H01','Pow_H02','Pow_H0C') #need to make sure we access powers in the same order for required and achieved power. power_mins <- unlist(case$power_mins[power_ind])
power_achieved <- case$performance[power_ind] if(any(power_achieved < power_mins)) { constraints_satisfied <- FALSE }else{ constraints_satisfied <- TRUE } return(constraints_satisfied) } #' binary search with arbitrary resolution #' #' The \code{\link{binsearch}} function in the \code{gtools} package searches only over integers. This is a wrapper for \code{\link{binsearch}} that also allows searching over the grid of points with distance \code{tol} between them. A target must also be entered (see \code{\link{binsearch}}) #' @param fun a function that determines the output over which we search (passed to \code{\link{binsearch}}) #' @param tol resolution of the grid over which to search for \code{target} (see \code{\link{binsearch}}) #' @param range a range over which to search for the input to \code{fun} #' @param ... passed to \code{\link{binsearch}} #' @export #' @import gtools #' @seealso \code{\link{binsearch}} #' @examples #' # best solution here is at x0, #' # which we can find with increasing precision #' x0 <- 10.241 #' binsearchtol( function(x) x-x0, target=0, range=c(0,2*x0) , tol=1.00) #' binsearchtol( function(x) x-x0, target=0, range=c(0,2*x0) , tol=0.10) #' binsearchtol( function(x) x-x0, target=0, range=c(0,2*x0) , tol=0.05) #' binsearchtol <- function(fun, tol=1, range, ...){ funTol <- function(x) fun(x*tol) soln <- binsearch(fun=funTol, range=range/tol, ...) soln$where <- soln$where*tol soln } #' Find the smallest total sample size such that a multi-stage trial meets desired constraints #' #' Takes as input a list of arguments (args) that define an adaptive trial (see \code{\link{buildTrial}} or \code{\link{optimizeTrial}}). This function adjusts the \code{n_total} argument in order to find the smallest maximum sample size that meets the power constraints specified in the \code{cases} argument. #' #' This function requires that the objective function contain a 'base' element, and a 'power_diffs' element that is nonnegative when power constraints are met. For example, see \code{\link{min_E_SS_power_constraints}}. #' #' @param args a list containing a subset of the arguments for the functions \code{\link{getEffBounds}} and \code{\link{simTrial}} (or comparable functions, see \code{trial_method} argument). This should include a FWER constraint. #' @param cases A list of power constraints, in the same format as those sent to \code{\link{optimizeTrial}} #' @param trial_method either 'cov' or 'MB' for Maurer-Bretz (2013). #' @param objective_fun see \code{\link{optimizeTrial}} #' @param min_n The smallest sample size to consider #' @param max_n The largest sample size to consider #' @param step_n The step size to consider when carrying out the binary search. For example, if \code{step_n = 10}, \code{min_n=0}, and \code{max_n=100}, the function will find the smallest n_total satisfying the supplied constraints, and that is also a multiple of 10. #' @param showiter passed to \code{\link{binsearch}} #' @export #' @references #' Maurer, W. and F. Bretz (2013). Multiple testing in group sequential trials using graphical approaches. \emph{Statistics in Biopharmaceutical Research.} #' @seealso \code{\link{min_n_feasible}}, \code{\link{feasibility_check}} #' @return A list containing #' \item{n}{The smallest feasible n_total} #' \item{soln}{Output from \code{\link{binsearch}}}. min_n_multistage<-function(args, cases, trial_method, objective_fun, min_n=1, max_n=min_n*1000, step_n=10, showiter=FALSE){ fun4search<-function(x){ args_x <- args args_x$n_total <- x
args_x <- cleanArgList(args_x,max_K=args_x$n_total) perf_list<-get_case_perf_obj(base_args=args_x, cases = cases, skip_penalty=FALSE, objective_fun=objective_fun, trial_method=trial_method) # get_perf_mat(perf_list) # output here generally increasing after power constraint is met # and constant before that constraints_satisfied <- all(perf_list$obj$power_diffs>=0) if( constraints_satisfied) return(perf_list$obj$base) if(!constraints_satisfied) return(0) #if statements required, or we get Inf * 0. } soln<-binsearchtol(fun=fun4search, range=c(min_n,max_n), tol=step_n, target=.5, showiter=showiter ) return(list( 'n'=soln$where[length(soln$where)], 'soln'=soln )) } #' @rdname feasibility_check #' @export #' @param min_n smallest sample size to consider for a 1-stage trial #' @param max_n largest sample size to consider #' @param step_n the step size for sample size. For example, if \code{step=20} then we will search in increments of 20 people. #' @param ... passed to \code{feasibility_check} min_n_feasible <- function(min_n=100, max_n = min_n * 10,step_n=5,showiter=FALSE, trial_method,...){ # In order to use binsearch, we need to insert a monotonic function # We use the sample size multiplied by whether power constraints are met # This is zero for n_total too small, positive for n_total big enough, and monotonically increasing. # Our target is 0.5, which will highlight the point at which our # monotonic function just begins to be positive. last_fc <- 0 #used for it's names later on fc_tracker <- data.frame() fun4search <- function(x){ fc<-feasibility_check(n_total=x, trial_method=trial_method, ...) if(fc[1]==1){ last_fc <<- fc fc_tracker <<- rbind(c('n_total'=x,fc),fc_tracker) } x*(fc[1]==1) } soln<-binsearchtol(fun=fun4search, range=c(min_n,max_n), tol=step_n, target=.5, showiter=showiter ) names(fc_tracker)<-c('n_total',names(last_fc)) out_n <- soln$value[2]

if(all(soln$value==0)){ warning('No solution found, max_n is too small. Returning max_n.') return(max_n) } if(soln$flag=='Lower Boundary'){
warning('Solution found at min_n. Consider decreasing min_n, as a smaller trial might also be feasible.')
out_n <- soln$value[1] } out <- (fc_tracker[fc_tracker$n_total==out_n,])[1,]
return(unlist(out))
}

#' @rdname feasibility_check
#' @export
#' @param step_multiplier the step size
#' @param showiter passed to \code{\link[gtools]{binsearch}}, determines whether to show progress.
max_power_feasible <- function(n_total,cases, p1, trial_method,step_multiplier=.01,showiter=FALSE,FWER=0.025,npoints_sqrt=10,...){

#The same approach is taken here as in min_n_feasible

#check inputs
if(step_multiplier<0 | step_multiplier>1) stop('step_multiplier must be between 0 and 1')

#Check initial/default n_total performance
fc1<-feasibility_check(n_total=n_total, p1, trial_method=trial_method,cases=cases, FWER=FWER,npoints_sqrt=npoints_sqrt)
if(fc1[1]==1){
warning('No search is necessary, initial n_total can achieve these power constraints. Returning original cases.')
return(list(
'multiplier'=1,cases=cases
))
}

#function to decrease min power required.
deflate_mins<-function(z,multiplier){
lapply(z,function(q) q*multiplier)
}
deflate_cases<-function(cases,multiplier){
cases_2<-cases
for(p in 1:length(cases)){
cases_2[[p]]$power_mins <- deflate_mins(cases_2[[p]]$power_mins,multiplier=multiplier)
}
cases_2
}

#Search for largest multiplier meeting power constraints
soln<-binsearchtol(fun=function(x){

if(x==1){

# recycle fc1 here to reduce redundant computation and potential
# for inconsistencies due to random monte carlo error
fc<-fc1

}else{

#multiply all power constraints by x
cases_2<-deflate_cases(cases, multiplier=x)
#check feasibility of resulting constraints
fc<-feasibility_check(n_total=n_total, FWER=FWER, p1=p1, trial_method=trial_method,cases=cases_2, ...)

}

out <- ((1/(x+1))+10)* ( fc[1] == 1)
out
#divide by x+1 to avoid dividing by zero
#For x <= optimal multiplier, out is decreasing (towards target=5) in x but bounded above 10
#For x >  optimal multiplier, out is zero
#We search for a target at 5
},
range=c(0,1),
tol=step_multiplier,
target=5.0,
showiter=showiter
)

m <- soln\$where[1]
cases_final <- deflate_cases(cases,multiplier=m)
return(list(
multiplier = m,
cases=cases_final
))
}

aaronjfisher/designOptim documentation built on April 3, 2018, 8:25 p.m.