R/zzz.R

In OneArmPhaseTwoStudy: Planning, Conducting, and Analysing Single-Arm Phase II Studies

  loadModule("simon", TRUE)
  loadModule("sub1", TRUE)


#' Sets the parameters for a given "simon"-object.
#'
#' @param s a "simon"-object which is generated by the function \code{\link{setupSimon}}.
#' @param alpha the maximal type I error rate.
#' @param beta the maximal type II error rate.
#' @param p0 the response probability under the null hypothesis.
#' @param p1 the response probability under the alternative hypothesis.
#'
#' @seealso \code{\link{setupSimon}}
#'
#' @examples
#' #Create "simon"-object.
#' simon <- setupSimon()
#' #Change parameters.
#' setSimonParams(simon, alpha = 0.1, beta = 0.2, p0 = 0.25, p1 = 0.45)
#' #Calculate designs for the given "simon"-object.
#' designs <- getSolutions(simon)$Solutions
#' designs
setSimonParams <- function(s, alpha = 0.05, beta = 0.05, p0 = 0.1, p1 = 0.3)
{
  stopifnot(class(alpha) == "numeric", alpha > 0, alpha < 1, class(beta) == "numeric", beta > 0, beta < 1)
  stopifnot(class(p0) == "numeric", p0 >= 0, p0 < 1, class(p1) == "numeric", p1 > 0, p1 < 1)
  stopifnot(class(s)[1] == "Rcpp_SimonDesign")

  s$setAlpha(alpha);
  s$setBeta(beta);
  s$setP0(p0);
  s$setP1(p1);

}

#' Creates a "simon"-object.
#'
#' Creates a "simon"-object which can be used in the function \code{\link{getSolutions}} to identify possible designs.
#'
#' @param alpha the maximal type I error rate.
#' @param beta the maximal type II error rate.
#' @param p0 the response probability under the null hypothesis.
#' @param p1 the response probability under the alternative hypothesis.
#'
#' @examples
#' #Create a "simon"-object
#' simon <- setupSimon()
#' #Calculate designs for the given "simon"-object.
#' designs <- getSolutions(simon)$Solutions
#' designs
setupSimon <-function(alpha = 0.05, beta = 0.05, p0 = 0.1, p1 = 0.3)
{
  stopifnot(class(alpha) == "numeric", alpha > 0, alpha < 1, class(beta) == "numeric", beta > 0, beta < 1)
  stopifnot(class(p0) == "numeric", p0 >= 0, p0 < 1, class(p1) == "numeric", p1 > 0, p1 < 1)

  #simon <- Module("simon", "OneArmPhaseTwoStudy", mustStart = TRUE);
  #s = new(simon$SimonDesign);
  s = new(SimonDesign)
  #simon_mod <- loadModule("simon")
  #s = new(simon_mod)
  setSimonParams(s,alpha,beta,p0,p1);
  s$aproximateMaxN();
  s;
}

#' Sets the parameters for a given "sub1"-object.
#'
#' @param sub1 a "sub1"-object which is generated by the function \code{\link{setupSub1Design}}.
#' @param alpha the maximal type I error rate.
#' @param beta the maximal type II error rate.
#' @param pc0 the response probability for the subset endpoint under the null hypothesis.
#' @param pt0 the response probability for the superset endpoint under the null hypothesis.
#' @param pc1 the response probability for the subset endpoint under the alternative hypothesis.
#' @param pt1 the response probability for the superset endpoint under the alternative hypothesis.
#'
#' @seealso \code{\link{setupSub1Design}}
#'
#' @examples
#' #Create "sub1"-object.
#' sub1 <- setupSub1Design()
#' #Change parameters.
#' setSub1Params(sub1, beta = 0.2, pc0 = 0.5, pt0 = 0.6)
#' #Calculate designs for the given "sub1"-object.
#' designs <- getSolutionsSub1(sub1)$Solutions
#' designs
setSub1Params <- function(sub1, alpha = 0.1, beta = 0.1, pc0 = 0.6, pt0 = 0.7, pc1 = 0.8, pt1 = 0.9 )
{
  stopifnot(class(alpha) == "numeric", alpha > 0, alpha < 1, class(beta) == "numeric", beta > 0, beta < 1)
  stopifnot(class(pc0) == "numeric", pc0 > 0, pc0 < 1, class(pt0) == "numeric", pt0 > 0, pt0 < 1)
  stopifnot(class(pc1) == "numeric", pc1 > 0, pc1 < 1, class(pt1) == "numeric", pt1 > 0, pt1 < 1)
  stopifnot(class(sub1)[1] == "Rcpp_Sub1Design")

  sub1$setAlpha(alpha);
  sub1$setBeta(beta);
  sub1$setPc0(pc0);
  sub1$setPt0(pt0);
  sub1$setPc1(pc1);
  sub1$setPt1(pt1);

}

#' Creates a "sub1"-object.
#'
#' Creates a "sub1"-object which can be used in the function \code{getSolutionsSub1} to identify possible designs.
#'
#' @param alpha the maximal type I error rate.
#' @param beta the maximal type II error rate.
#' @param pc0 the response probability for the subset endpoint under the null hypothesis.
#' @param pt0 the response probability for the superset endpoint under the null hypothesis.
#' @param pc1 the response probability for the subset endpoint under the alternative hypothesis.
#' @param pt1 the response probability for the superset endpoint under the alternative hypothesis.
#'
#' @examples
#' #Create "sub1"-object.
#' sub1 <- setupSub1Design()
#' #Calculate designs for the given "sub1"-object.
#' designs <- getSolutionsSub1(sub1)$Solutions
#' designs
setupSub1Design <-function(alpha = 0.1, beta = 0.1, pc0 = 0.6, pt0 = 0.7, pc1 = 0.8, pt1 = 0.9)
{
  stopifnot(class(alpha) == "numeric", alpha > 0, alpha < 1, class(beta) == "numeric", beta > 0, beta < 1)
  stopifnot(class(pc0) == "numeric", pc0 > 0, pc0 < 1, class(pt0) == "numeric", pt0 > 0, pt0 < 1)
  stopifnot(class(pc1) == "numeric", pc1 > 0, pc1 < 1, class(pt1) == "numeric", pt1 > 0, pt1 < 1)

  #sub1 <- Module("sub1", "OneArmPhaseTwoStudy", mustStart = TRUE);
  #sub1 = new(sub1$Sub1Design);
  sub1 = new(Sub1Design)
  #sub1_mod <-loadModule("sub1")
  #sub1 = new(sub1_mod)
  setSub1Params(sub1, alpha, beta, pc0, pt0, pc1, pt1);
  sub1$aproximateMaxN();
  sub1;
}

#' Helper function for \code{\link{getSolutions}} and \code{\link{getSolutionsSub1}}.
#'
#' Transfers the found designs of a "design"-object to a dataframe. This function is a helper-function used in
#' \code{\link{getSolutions}} and \code{\link{getSolutionsSub1}} to structure the found designs for a given parameter set specified
#' in a "sub1"-object or "simon"-object.
#' The return value is a list containing a dataframe for the found designs and multiple dataframes containing the results for (non-)stochastic curtailment, if present.
#'
#' @param designObject either a "sub1"-object or a "simon"-object containing the designs for a given parameter set specified inside the "designObject".
#' @param useCurtailment boolean value dedicating the use of (non-)stochastic curtailment. This parameter indicates if that information should also gathered.
toDataframe <-function(designObject, useCurtailment =F)
{
  stopifnot(class(designObject)[1] == "Rcpp_SimonDesign" | class(designObject)[1] == "Rcpp_Sub1Design")
  stopifnot(class(useCurtailment) == "logical")

  tmp = designObject$getResultsForR();
  df = data.frame();
  curList = list();
  if(length(tmp)>0)
  {
    df = data.frame()
    if(length(tmp)>=1)
    {
      for(i in 1:length(tmp))
      {
        df = rbind(df, tmp[[i]]$Solution);
        if(useCurtailment)
        {
          temp = tmp[[i]]$Curtailment;
          if(length(temp) > 0)
          {
            tmpDF = temp[[1]]$Details;
            stoppingRulesList <- list();
            stoppingRulesList[[paste("Stoppingrules_for_Row:", 1, sep="")]] <- temp[[1]]$Stoppingrules;
            if(length(temp) > 1)
            {
              for(j in 2: length(temp))
              {
                tmpDF = rbind(tmpDF, temp[[j]]$Details)
                stoppingRulesList[[paste("Stoppingrules_for_Row:", j, sep="")]] <- temp[[j]]$Stoppingrules;
              }
            }
            curList[[paste("CurResultsForID:",df$ID[i],sep="")]] <- tmpDF;
            curList[[paste("StoppingrulesForID:", df$ID[i], sep="")]] <- stoppingRulesList;
          }
        }
      }
    }
    if(class(designObject)[1] == "Rcpp_SimonDesign"){
      df[,c("enP0", "enP1")] <- round(df[,c("enP0", "enP1")], digits = 2)
      df[,c("petP0", "petP1", "Alpha", "Beta", "Admiss_Start", "Admiss_End")] <-
        round(df[,c("petP0", "petP1", "Alpha", "Beta", "Admiss_Start", "Admiss_End")], digits = 4)
    }else if(class(designObject)[1] == "Rcpp_Sub1Design"){
      df[,"enP0"] <- round(df[,"enP0"], digits = 2)
      df[,c("petP0", "Alpha", "Beta", "Admiss_Start", "Admiss_End")] <-
        round(df[,c("petP0", "Alpha", "Beta", "Admiss_Start", "Admiss_End")], digits = 4)
    }
  }
  solutionresults = list(Solutions =df, Curtailment_Results = curList);
  solutionresults
}

#' Returns designs for a given "simon"-object (see \code{\link{setupSimon}})
#'
#' getSolutions uses a "simon"-object to calculate two-stage designs as they were described by Simon.
#'
#' @param simon a "simon"-object which will be used to calculate designs.
#' @param useCurtailment boolean value determining whether (non-)stochastic curtailment is used.
#' @param curtail_All boolean value; if true the effect of (non-)stochastic curtailment will be calculated for different cut points in 0.05 steps starting with the value of the parameter "cut".
#' @param cut sets the "cut point" used to calculate the effect of (non-)stochastic curtailment. A study is stopped if the conditional power falls below the value of "cut".
#' @param replications number of simulations to estimate the effect of (non-)stochastic curtailment.
#' @param upperBorder maximal possible value for n. If set to sero (default) the programm will aproximate a upper border automaticly.
#'
#' @seealso \code{\link{setupSimon}}
#'
#' @references Simon, R. (1989): Optimal two-stage designs for phase II clinical trials. Controlled Clinical Trials 10,1-10.
#' @references Kunz C.U., Kieser M (2012): Curtailment in single-arm two-stage phase II oncology trials. Biometrical Journal 54, 445-456
#'
#' @examples
#' # Example 1: Using the default values
#' designs <- getSolutions()
#' designs <- designs$Solutions
#' designs
#' \donttest{
#' # Example 2: Setting up a "simon"-object, then calculate designs
#' simon <- setupSimon(alpha = 0.1, beta = 0.2, p0 = 0.3, p1 = 0.5)
#' designs <- getSolutions(simon)$Solutions
#' designs
#'
#' # Esample 3: Calculating designs and simulating the influence of
#' # stochastic curtailment for each design.
#' simon <- setupSimon(alpha = 0.1, beta = 0.2, p0 = 0.3, p1 = 0.5)
#' designs <- getSolutions(simon, useCurtailment = TRUE, curtail_All = TRUE, cut = 0.3)
#' #List containing the found designs, the influence of stochastic curtailment
#' # and the regarding stopping rules.
#' designs
#' }
getSolutions <-function(simon = setupSimon(), useCurtailment = FALSE, curtail_All = FALSE, cut = 0.0, replications = 10000, upperBorder = 0)
{
  stopifnot(class(useCurtailment) == "logical", class(curtail_All) == "logical")
  stopifnot(class(cut) == "numeric", cut >= 0, cut <= 1, class(replications) == "numeric", replications > 0, class(upperBorder) == "numeric", upperBorder >= 0)
  stopifnot(class(simon)[1] == "Rcpp_SimonDesign")

  if(upperBorder > 0)
  {
    simon$setMaxN(upperBorder)
  }else{
    simon$aproximateMaxN()
  }
  simon$calculateStudySolutions();
  if(useCurtailment)
  {
    numSol = simon$getSolutionCount()
    for(i in 0:(numSol -1))
    {
      simon$calculateSC(i,cut, replications, curtail_All)
    }
  }
  toDataframe(simon, useCurtailment)
}



#' Calculates designs for a given "sub1"-object.
#'
#' By iterating over all possible values for "r1", "n1", "r", "s" and "n" designs for a given "sub1"-object are found.
#' Proceeding to the second stage of the study more than "r1" responses among the first "n1" patients in the subset endpoint are needed.
#' Rejecting the null hypothesis more than "r" responses in the subset endpoint
#' or more than "s" responses in the superset endpoint among "n" patients are needed.
#'
#' @param sub1 a "sub1"-object which will be used to calculate fitting designs
#' @param skipS boolean value; skips the iteration over "s" at certian points to improve calculation speed (finds less designs)
#' @param skipR boolean value; skips the iteration over "r" at certian points to improve calculation speed (finds less designs)
#' @param skipN1 boolean value; skips the iteration over "n1" at certian points to improve calculation speed (finds less designs and it is impossible to determine the optimalization criteria of the found designs)
#' @param lowerBorder sets a minimal value for "n" (number of patients to be recruited)
#' @param upperBorder sets a maximal value for "n" (number of patients to be recruited)
#' @param useCurtailment determines if the effect of (non-)stochastic curtailment should also be calculated for the found designs
#' @param curtailAll boolean value; if true the effect of (non-)stochastic curtailment will be calculated for different cut points in 0.05 steps starting with the value of the parameter "cut".
#' @param cut sets the "cut point" used to calculate the effect of (non-)stochastic curtailment. A study is stopped if the conditional power falls below the value of "cut".
#' @param replications number of simulations to estimate the effect of (non-)stochastic curtailment.
#'
#' @references Kunz C.U., Kieser M (2012): Curtailment in single-arm two-stage phase II oncology trials. Biometrical Journal 54, 445-456
#'
#' @seealso \code{\link{setupSub1Design}}
#'
#' @examples
#' # Example 1: Using the default values
#' sub1 <- setupSub1Design()
#' getSolutionsSub1(sub1)
#' \donttest{
#' # Example 2: Setting up a "sub1"-object, then calculating designs
#' sub1 <- setupSub1Design(alpha = 0.1, beta = 0.2, pc0 = 0.3, pt0 = 0.4)
#' designs <- getSolutionsSub1(sub1)$Solutions
#' designs
#'
#' # Example 2: Calculating designs and simulating the influence of stochastic curtailment
#' # for each design.
#' sub1 <- setupSub1Design(alpha = 0.1, beta = 0.2, pc0 = 0.3, pt0 = 0.4)
#' designs <- getSolutionsSub1(sub1, useCurtailment = TRUE, curtailAll = TRUE, cut = 0.3)
#' #Contains the found designs, the influence of stochastic curtailment
#' #and the regarding stopping rules .
#' designs
#' }
getSolutionsSub1 <-function(sub1 = setupSub1Design(), skipS = TRUE, skipR = TRUE, skipN1 = TRUE, lowerBorder = 0, upperBorder = 0, useCurtailment = FALSE, curtailAll = FALSE, cut = 0.00, replications = 1000)
{
  stopifnot(class(skipS) == "logical", class(skipR) == "logical", class(skipN1) == "logical")
  stopifnot(class(useCurtailment) == "logical", class(curtailAll) == "logical")
  stopifnot(class(cut) == "numeric", cut >= 0, cut <= 1, class(replications) == "numeric", replications > 0)
  stopifnot(class(lowerBorder) == "numeric", lowerBorder >= 0, class(upperBorder) == "numeric", upperBorder >= 0, upperBorder >= lowerBorder)
  stopifnot(class(sub1)[1] == "Rcpp_Sub1Design")

  #helper function
  ctp <- function(r1, n1, r, s, n, pc0, pt0, alpha){
    sim <- setupSimon();
    ep1_alpha <- sim$calcAlpha(n1, r1, n, r, pc0)
    if(ep1_alpha <= alpha){
      ep2_pv <-binom.test(s+1, n, pt0, alternative = "g")$p.value
      if(ep2_pv <= alpha){
        return(T)
      }
    }else{
      return(F)
    }
    return(F)
  }


  sub1$calculateStudySolutions(skipS, skipR, skipN1, lowerBorder, upperBorder);
  if(useCurtailment)
  {
    numSol = sub1$getSolutionCount()
    for(i in 0:(numSol -1))
    {
      sub1$calculateSC(i, cut, replications, curtailAll)
    }
  }
  result <- toDataframe(sub1, useCurtailment)
  res_sol <- result$Solutions
  ctp_res <- mapply(ctp, r1 = res_sol$r1, n1 = res_sol$n1, r= res_sol$r, s = res_sol$s, n = res_sol$n,
                    pc0 = res_sol$pc0, pt0 = res_sol$pt0, MoreArgs = list(alpha = sub1$getAlpha()))
  result$Solutions <- cbind(result$Solutions, ClosedTestProcedure = ctp_res)
  result
}

#' Calculates the conditional error.
#'
#' Calculates the conditional error at the interim analysis for a given Simon's design with "k" responses.
#'
#' @param design a dataframe containing all critical values for a Simon's two-stage design defined by the colums "r1", "n1", "r", "n" and "p0".
#' \itemize{
#'    \item r1 = critical value for the first stage (more than "r1" responses needed to proceed to the second stage).
#'    \item n1 = number of patients enrolled in the first stage.
#'    \item r = critical value for the whole trial (more than "r" responses needed at the end of the study to reject the null hypothesis).
#'    \item n = number of patients enrolled in the whole trial.
#'    \item p0 = response probability under the null hypothesis
#' }
#'
#' @param k number of responses observed at the interim analysis.
#'
#' @examples
#' design <- getSolutions()$Solutions[1,]
#' conditional_error <- getCE(design, 4)
getCE <-function(design, k)
{
  stopifnot(class(design) == "data.frame")
  dgn_names = names(design)
  stopifnot("r1" %in% dgn_names, "r" %in% dgn_names, "p0" %in% dgn_names, "n1" %in% dgn_names, "n" %in% dgn_names)
  stopifnot(class(k) == "numeric" | class(k) == "integer", k >= 0)

  l1 = design$r1
  u1 = design$r
  l2 = design$r
  pi0 = design$p0
  n2 = design$n - design$n1
  if(k <= l1)
  {
    return(0)
  }else if (k > u1)
  {
    return(1)
  }else
  {
    res = 1 - pbinom(l2 - k, n2, pi0)
    return(res)
  }
}

#liefert den p-wert (fuer eine stufe)

#' Calculates the p-value (binomial test).
#'
#' Helper-function for the function \code{\link{getCP}}
#'
#' @param l number of responses
#' @param pi0 response probability under the null hypothesis
#' @param n2 number of enrolled patients
#'
#' @seealso \code{\link{getCP}}
getP <-function(l,pi0,n2)
{
  stopifnot(class(l) == "numeric"  | class(l) == "integer", l >=0)
  stopifnot(class(pi0) == "numeric", pi0 > 0, pi0 <= 1)
  stopifnot(class(n2) == "numeric", n2 >=0)

  res = 1- pbinom(l-1, n2, pi0)
  res
}

#' Calculates the number of responses needed for the second stage.
#'
#' Calculates the number of responses needed for the second stage of a Simon's two-stage design if the flexible extension is chosen in the planning phase.
#'
#' @param ce conditional error for the second stage.
#' @param p0 probability for a response under the null hypothesis.
#' @param n2 sample size for the second stage.
#'
#' @seealso \code{\link{getD_proportionally}}, \code{\link{getD_equally}}, \code{\link{getD_distributeToOne}}, \code{\link{getD_none}}
#'
#' @examples
#' #Calculate a Simon's two-stage design
#' design <- getSolutions()$Solutions[3,] #minimax-design for the default values.
#' #Get the conditional error values using proportionally "rest"-alpha spending.
#' ce_df <- getD_proportionally(design, 0.05)
#' #Assume 5 responses were observed in the interim analysis.
#' ce <- ce_df[5+1,]$ce # conditional error for 5 responses is listed in the 6th row of "ce_df"
#' #Calculate the number of patients needed in the second stage.
#' n2 <- design$n - design$n1
#' r2 <- get_r2_flex(ce, design$p0, n2)
#' r2
#' #Assume 10 patients more should be recruited in the second stage.
#' #(This changes the number of needed responses.)
#' n2 <- n2 + 10
#' r2 <- get_r2_flex(ce, design$p0, n2)
#' r2
get_r2_flex <-function(ce, p0, n2)
{
  stopifnot(class(ce) == "numeric", ce >= 0, ce <= 1)
  stopifnot(class(p0) == "numeric", p0 >= 0, p0 <= 1)
  stopifnot(class(n2) == "numeric", n2 >=0)

  p_tmp <- 1
  r2 <- -1
  while(p_tmp > ce)
  {
    r2 <- r2 + 1
    p_tmp <- getP(r2, p0, n2)
  }
  r2
}


#' Calculates the conditional power.
#'
#' Calculates the conditional power for a given Simon's two-stage design in the interim analysis if the number of patients which should be enrolled in the second stage is altert to "n2".
#'
#' @param n2 number of patients to be enrolled in the second stage of the study.
#' @param p1 response probability under the alternative hypothesis
#' @param design a dataframe containing all critical values for a Simon's two-stage design defined by the colums "r1", "n1", "r", "n" and "p0".
#' \itemize{
#'    \item r1 = critical value for the first stage (more than r1 responses needed to proceed to the second stage).
#'    \item n1 = number of patients enrolled in the first stage.
#'    \item r = critical value for the whole trial (more than r responses needed at the end of the study to reject the null hypothesis).
#'    \item n = number of patients enrolled in the whole trial.
#'    \item p0 = response probability under the null hypothesis.
#' }
#'
#' @param k number of responses observed at the interim analysis.
#' @param mode a value out of \{0,1,2,3\} dedicating the methode spending the "rest alpha" (difference between nominal alpha level and actual alpha level for the given design).
#' \itemize{
#'    \item 0 = "rest alpha" is not used.
#'    \item 1 = "rest alpha" is spent proportionally.
#'    \item 2 = "rest alpha" is spent equally.
#'    \item 3 = "rest alpha" is spent only to the worst case scenario (minimal number of responses at the interim analysis so that the study can proceed to the second stage).
#' }
#'
#' @param alpha overall significance level the trial was planned for.
#'
#' @references Englert S., Kieser M. (2012): Adaptive designs for single-arm phase II trials in oncology. Pharmaceutical Statistics 11,241-249.
#'
#' @seealso \code{\link{getN2}}
#'
#' @examples
#' #Calculate a Simon's two-stage design
#' design <- getSolutions()$Solutions[3,] #minimax-design for the default values.
#'
#' #Assume 3 responses were observed in the interim analysis.
#' #Therefore the conditional power is only about 0.55.
#' #In order to raise the conditional power to 0.8 "n2" has to be increased.
#'
#' #get the current "n2"
#' n2 <- design$n - design$n1
#'
#' #set k to 3 (only 3 responses observed so far)
#' k = 3
#'
#' #get the current conditional power
#' cp <- getCP(n2, design$p1, design, k, mode = 1, alpha = 0.05)
#' cp
#'
#' \donttest{
#' #increase n2 until the conditional power is larger than 0.8
#' while(cp < 0.8){
#' n2 <- n2 + 1
#' # Assume we spent the "rest alpha" proportionally (in the planning phase)
#' # therefore we set "mode = 1".
#' cp <- getCP(n2, design$p1, design, k, mode = 1, alpha = 0.05)
#' }
#'
#' n2
#' }
getCP <-function(n2,p1, design, k, mode = 0, alpha = 0.05)
{
  stopifnot(class(n2) == "numeric", n2 >= 0)
  stopifnot(class(p1) == "numeric", p1 > 0, p1 <= 1)
  stopifnot(class(design) == "data.frame")
  dgn_names = names(design)
  stopifnot("r1" %in% dgn_names, "r" %in% dgn_names, "p0" %in% dgn_names, "n1" %in% dgn_names, "n" %in% dgn_names)
  stopifnot(class(k) == "numeric" | class(k) == "integer", k >=0)
  stopifnot(class(as.numeric(mode)) == "numeric", mode %in% 0:3)
  stopifnot(class(alpha) == "numeric", alpha > 0, alpha <= 1)

  p0 = design$p0
  cp = 0
  ce = 0;
  if(mode == 0)
  {
    ce = getCE(design, k)
  }else if (mode == 1) {
    ce = getD_proportionally(design, alpha)
    ce = ce$ce[k + 1]
  }else if (mode == 2) {
    ce = getD_equally(design, alpha)
    ce = ce$ce[k + 1]
  }else if (mode == 3) {
    ce = getD_distributeToOne(design, alpha)
    ce = ce$ce[k + 1]
  }

  for(l in 0:n2)
  {
    if(getP(l,p0,n2) <= ce)
    {
      cp = cp + dbinom(l, n2, p1)
    }
  }
  cp
}

#' Calculates the number of patients which should be enrolled in the second stage.
#'
#' Calculates the number of patients which should be enrolled in the second stage if the conditional power should be altert to "cp".
#'
#' @param cp conditional power to which the number of patients for the second stage should be adjusted.
#' @param p1 response probability under the alternative hypothesis.
#' @param design a dataframe containing all critical values for a Simon's two-stage design defined by the colums r1, n1, r, n and p0.
#' \itemize{
#'    \item r1 = critical value for the first stage (more than r1 responses needed to proceed to the second stage).
#'    \item n1 = number of patients enrolled in the first stage.
#'    \item r = critical value for the whole trial (more than r responses needed at the end of the study to reject the null hypothesis).
#'    \item n = number of patients enrolled in the whole trial.
#'    \item p0 = response probability under the null hypothesis.
#' }
#' @param k number of responses observed at the interim analysis.
#' @param mode a value out of \{0,1,2,3\} dedicating the methode spending the "rest alpha" (difference between nominal alpha level and actual alpha level for the given design).
#' \itemize{
#'    \item 0 = "rest alpha" is not used.
#'    \item 1 = "rest alpha" is spent proportionally.
#'    \item 2 = "rest alpha" is spent equally.
#'    \item 3 = "rest alpha" is spent only to the worst case scenario (minimal number of responses at the interim analysis so that the study can proceed to the second stage).
#' }
#' @param alpha overall significance level the trial was planned for.
#'
#' @references Englert S., Kieser M. (2012): Adaptive designs for single-arm phase II trials in oncology. Pharmaceutical Statistics 11,241-249.
#'
#' @examples
#' #Calculate a Simon's two-stage design
#' design <- getSolutions()$Solutions[3,] #minimax-design for the default values.
#'
#' #Assume we only observed 3 responses in the interim analysis.
#' #Therefore the conditional power is only about 0.55.
#' #In order to raise the conditional power to 0.8 "n2" has to be increased.
#'
#' #set k to 3 (only 3 responses observed so far)
#' k = 3
#'
#' # Assume we spent the "rest alpha" proportionally in the planning phase
#' # there for we set "mode = 1".
#' n2 <- getN2(cp = 0.8, design$p1, design, k, mode = 1, alpha = 0.05)
#' n2
getN2 <-function(cp, p1, design, k, mode = 0, alpha = 0.05)
{
  stopifnot(class(cp) == "numeric", cp > 0 , cp <= 1)
  stopifnot(class(p1) == "numeric", p1 > 0, p1 <= 1)
  stopifnot(class(design) == "data.frame")
  dgn_names = names(design)
  stopifnot("r1" %in% dgn_names, "r" %in% dgn_names, "p0" %in% dgn_names, "n1" %in% dgn_names, "n" %in% dgn_names)
  stopifnot(class(k) == "numeric" | class(k) == "integer", k >=0)
  stopifnot(class(mode) == "numeric", mode %in% 0:3)
  stopifnot(class(alpha) == "numeric", alpha > 0, alpha <= 1)

  if(cp == 1)
  {
    return(0)
  }
  cpInt = 0
  n2 = 0
  while(cpInt < cp)
  {
    n2 = n2 +1
    cpInt = getCP(n2, p1, design , k, mode, alpha)
  }
  n2
}


#' Get the conditional errors proportionally.
#'
#' Calculates the conditional error for all possible outcomes at the interim analysis (different number of responses)
#' spending "rest alpha" (difference between nominal alpha level and actual alpha level) proportionally.
#'
#' @param design a dataframe containing all critical values for a Simon's two-stage design defined by the colums "r1", "n1", "r", "n" and "p0".
#' \itemize{
#'    \item r1 = critical value for the first stage (more than "r1" responses needed to proceed to the second stage).
#'    \item n1 = number of patients enrolled in the first stage.
#'    \item r = critical value for the whole trial (more than "r" responses needed at the end of the study to reject the null hypothesis).
#'    \item n = number of patients enrolled in the whole trial.
#'    \item p0 = response probability under the null hypothesis.
#' }
#' @param alpha overall significance level the trial was planned for.
#'
#' @references Englert S., Kieser M. (2012): Adaptive designs for single-arm phase II trials in oncology. Pharmaceutical Statistics 11,241-249.
#'
#' @seealso \code{\link{getD_equally}}, \code{\link{getD_distributeToOne}}, \code{\link{getD_none}}
#'
#' @examples
#' #Calculate a Simon's two-stage design
#' design <- getSolutions()$Solutions[3,] #minimax-design for the default values.
#'
#' ce_prop <- getD_proportionally(design, 0.05)
#' ce_prop
getD_proportionally <-function(design, alpha)
{
  stopifnot(class(design) == "data.frame")
  dgn_names = names(design)
  stopifnot("r1" %in% dgn_names, "r" %in% dgn_names, "p0" %in% dgn_names, "n1" %in% dgn_names, "n" %in% dgn_names)
  stopifnot(class(alpha) == "numeric", alpha > 0, alpha <= 1)

  ces = data.frame()
  for(k in 0: design$n1)
  {
    ces = rbind(ces, data.frame(k = k, ce = getCE(design,k)))
  }
  alphaR= alpha - design$Alpha
  sum = 0
  cesTmp = ces[ces$ce > 0 & ces$ce <1,]
  n2 = design$n - design$n1

  for(i in 1:nrow(cesTmp))
  {
    sum = sum + dbinom(cesTmp[i,]$k, design$n1, design$p0)
  }
  DCE = ces[ces$ce == 0,]
  DCE = rbind(DCE, cesTmp + (alphaR / sum))
  DCE = rbind(DCE, ces[ces$ce == 1,])
  DCE
}


#' Get the conditional errors equally.
#'
#' Calculates the conditional error for all possible outcomes at the interim analysis (different number of responses)
#' spending the "rest alpha" (difference between nominal alpha level and actual alpha level) equally.
#'
#' @param design a dataframe containing all critical values for a Simon's two-stage design defined by the colums r1, n1, r, n and p0.
#' \itemize{
#'    \item r1 = critical value for the first stage (more than r1 responses needed to proceed to the second stage).
#'    \item n1 = number of patients enrolled in the first stage.
#'    \item r = critical value for the whole trial (more than r responses needed at the end of the study to reject the null hypothesis).
#'    \item n = number of patients enrolled in the whole trial.
#'    \item p0 = response probability under the null hypothesis.
#' }
#' @param alpha overall significance level the trial was planned for.
#'
#' @references Englert S., Kieser M. (2012): Adaptive designs for single-arm phase II trials in oncology. Pharmaceutical Statistics 11,241-249.
#'
#' @seealso \code{\link{getD_proportionally}}, \code{\link{getD_distributeToOne}}, \code{\link{getD_none}}
#'
#' @examples
#' #Calculate a Simon's two-stage design
#' design <- getSolutions()$Solutions[3,] #minimax-design for the default values.
#'
#' ce_equally <- getD_equally(design, 0.05)
#' ce_equally
getD_equally <-function(design, alpha)
{
  stopifnot(class(design) == "data.frame")
  dgn_names = names(design)
  stopifnot("r1" %in% dgn_names, "r" %in% dgn_names, "p0" %in% dgn_names, "n1" %in% dgn_names, "n" %in% dgn_names)
  stopifnot(class(alpha) == "numeric", alpha > 0, alpha <= 1)

  ces = data.frame()
  for(k in 0: design$n1)
  {
    ces = rbind(ces, data.frame(k = k, ce = getCE(design,k)))
  }
  alphaR= alpha - design$Alpha
  sum = 0
  cesTmp = ces[ces$ce > 0 & ces$ce <1,]
  numRows = nrow(cesTmp)
  for(i in 1:numRows)
  {
    cesTmp[i,]$ce = cesTmp[i,]$ce + (alphaR / numRows) / dbinom(cesTmp[i,]$k, design$n1, design$p0)
    if(cesTmp[i,]$ce > 1){
      cesTmp[i,]$ce = 1
    }
  }
  DCE = ces[ces$ce == 0,]
  DCE = rbind(DCE, cesTmp)
  DCE = rbind(DCE, ces[ces$ce == 1,])
  DCE
}


#' Get the conditional errors.
#'
#' Calculates the conditional error for all possible outcomes at the interim analysis (different number of responses)
#' spending the "rest alpha" (difference between nominal alpha level and actual alpha level) only to increase the worst case (smallest conditional error value that is not equal to 0).
#'
#' @param design a dataframe containing all critical values for a Simon's two-stage design defined by the colums "r1", "n1", "r", "n" and "p0".
#' \itemize{
#'    \item r1 = critical value for the first stage (more than r1 responses needed to proceed to the second stage).
#'    \item n1 = number of patients enrolled in the first stage.
#'    \item r = critical value for the whole trial (more than r responses needed at the end of the study to reject the null hypothesis).
#'    \item n = number of patients enrolled in the whole trial.
#'    \item p0 = response probability under the null hypothesis.
#' }
#' @param alpha overall significance level the trial was planned for.
#'
#' @references Englert S., Kieser M. (2012): Adaptive designs for single-arm phase II trials in oncology. Pharmaceutical Statistics 11,241-249.
#'
#' @seealso \code{\link{getD_proportionally}}, \code{\link{getD_equally}}, \code{\link{getD_none}}
#'
#' @examples
#' #Calculate a Simon's two-stage design
#' design <- getSolutions()$Solutions[3,] #minimax-design for the default values.
#'
#' ce_toOne <- getD_distributeToOne(design, 0.05)
#' ce_toOne
getD_distributeToOne <-function(design, alpha)
{
  stopifnot(class(design) == "data.frame")
  dgn_names = names(design)
  stopifnot("r1" %in% dgn_names, "r" %in% dgn_names, "p0" %in% dgn_names, "n1" %in% dgn_names, "n" %in% dgn_names)
  stopifnot(class(alpha) == "numeric", alpha > 0, alpha <= 1)

  ces = data.frame()
  for(k in 0: design$n1)
  {
    ces = rbind(ces, data.frame(k = k, ce = getCE(design,k)))
  }
  alphaR= alpha - design$Alpha
  sum = 0
  cesTmp = ces[ces$ce > 0 & ces$ce <1,]

  cesTmp[1,]$ce = cesTmp[1,]$ce + alphaR / dbinom(cesTmp[1,]$k, design$n1, design$p0)

  DCE = ces[ces$ce == 0,]
  DCE = rbind(DCE, cesTmp)
  DCE = rbind(DCE, ces[ces$ce == 1,])
  DCE
}

#' Get the conditional errors.
#'
#' Calculates the conditional error for all possible outcomes at the interim analysis (different number of responses)
#' using no "rest alpha" spending (difference between nominal alpha level and actual alpha level).
#'
#'
#' @param design a dataframe containing all critical values for a Simon's two-stage design defined by the colums "r1", "n1", "r", "n" and "p0".
#' \itemize{
#'    \item r1 = critical value for the first stage (more than r1 responses needed to proceed to the second stage).
#'    \item n1 = number of patients enrolled in the first stage.
#'    \item r = critical value for the whole trial (more than r responses needed at the end of the study to reject the null hypothesis).
#'    \item n = number of patients enrolled in the whole trial.
#'    \item p0 = response probability under the null hypothesis.
#' }
#'
#' @references Englert S., Kieser M. (2012): Adaptive designs for single-arm phase II trials in oncology. Pharmaceutical Statistics 11,241-249.
#'
#' @seealso \code{\link{getD_proportionally}}, \code{\link{getD_equally}}, \code{\link{getD_distributeToOne}}
#'
#' @examples
#' #Calculate a Simon's two-stage design
#' design <- getSolutions()$Solutions[3,] #minimax-design for the default values.
#'
#' ce_toOne <- getD_none(design)
#' ce_toOne
getD_none <-function(design)
{
  stopifnot(class(design) == "data.frame")
  dgn_names = names(design)
  stopifnot("r1" %in% dgn_names, "r" %in% dgn_names, "p0" %in% dgn_names, "n1" %in% dgn_names, "n" %in% dgn_names)

  ces = data.frame()
  for(k in 0: design$n1)
  {
    ces = rbind(ces, data.frame(k = k, ce = getCE(design,k)))
  }
  ces
}

#' Calculates the "uniformly minimal variance unbiased estimator".
#'
#' Calculates the "uniformly minimal variance unbiased estimator" (UMVUE) for the true response rate
#' based on the approach of Grishick, Mosteller and Savage.
#'
#' @param k overall observed responses.
#' @param r1 critical value for the first stage.
#' @param n1 sample size for the first stage.
#' @param n overall sample size.
#'
#' @references Girshick MA, Mosteller F, and Savage LJ (1946): Unbiased estimates for certain binomial sampling problems with applications. Annals of Mathematical Statistics, 17(1):13-23.
#'
#' @examples
#' #Calculate a Simon's two-stage design
#' design <- getSolutions()$Solutions[3,] #minimax-design for the default values.
#'
#' #Assume 9 responses were observed in the whole trial.
#' k = 9
#'
#' umvue <- get_UMVUE_GMS(k, design$r1, design$n1, design$n)
get_UMVUE_GMS <- function(k, r1, n1, n)
{
  stopifnot(class(k) == "numeric" | class(k) == "integer", k >= 0)
  stopifnot(class(r1) == "numeric", r1 >= 0)
  stopifnot(class(n1) == "numeric", n1 >= 0, n1 > r1)
  stopifnot(class(n) == "numeric", n >= 0, n > n1)

  umvue <- numeric(0)
  if(k <= r1)
  {
    umvue <- k/n1
  }else
  {
    upper_boundry_dividend <- min(k-1, n1-1)
    lower_boundry_dividend <- max(r1, k - 1 - n + n1)

    upper_boundry_divisor <- min(k, n1)
    lower_boundry_divisor <- max(r1 + 1, k - n + n1)

    i_upper <- lower_boundry_dividend : upper_boundry_dividend
    i_lower <- lower_boundry_divisor : upper_boundry_divisor

    umvue <- sum(choose(n1-1,i_upper) * choose(n-n1, k-1-i_upper)) / sum(choose(n1, i_lower) * choose(n-n1, k-i_lower))
  }
  umvue
}

#' Calculates the "uniformly minimal variance unbiased estimator".
#'
#' Calculates the "uniformly minimal variance unbiased estimator" (UMVUE) for the true response rate only for the superset endpoint
#' (response rate superset endpoint minus response rate subset endpoint) in a subset design.
#'
#' @param t observed responses in the subset endpoint.
#' @param u observed responses in the superset endpoint.
#' @param r1 critical value for the first stage.
#' @param n1 sample size for the first stage.
#' @param n overall sample size.
#'
#' @examples
#' #Setup "sub1"-object
#' sub1 <- setupSub1Design(pc0 = 0.5, pt0 = 0.6)
#'
#' #Calculate a subset design
#' design <- getSolutionsSub1(sub1, skipN1 = FALSE)$Solutions[4,]
#'
#' #Assume 9 responses in the subset endpoint and 13 responses in the superset endpoint were observed.
#' t = 9
#' u = 13
#' umvue_second <- get_UMVUE_GMS_subset_second_only(t, u, design$r1, design$n1, design$n)
get_UMVUE_GMS_subset_second_only <- function(t, u, r1, n1, n)
{
  stopifnot(class(t) == "numeric", t >= 0)
  stopifnot(class(u) == "numeric", u >= 0, u >= t)
  stopifnot(class(r1) == "numeric", r1 >= 0)
  stopifnot(class(n1) == "numeric", n1 >= 0, n1 > r1)
  stopifnot(class(n) == "numeric", n >= 0, n > n1)

  umvue <- numeric(0)
  if(t <= r1)
  {
    sum_index <- 0:r1

    umvue <- sum(choose(n1-1, sum_index) * choose(n1 - 1 - sum_index, u - 1 -1)) / sum(choose(n1, sum_index) * choose(n1 - sum_index, u - t))
  }else
  {
    upper_boundry_dividend <- min(t, n1 - 1)
    lower_boundry_dividend <- max(r1 + 1, t - n + n1)

    upper_boundry_divisor <- min(t, n1)
    lower_boundry_divisor <- max(r1 + 1, t - n + n1)

    i_upper <- lower_boundry_dividend : upper_boundry_dividend
    i_lower <- lower_boundry_divisor : upper_boundry_divisor

    umvue <- sum(choose(n1 - 1, i_upper) * choose(n - n1, t - i_upper) * choose(n - t - 1, u - t -1)) /
             sum(choose(n1, i_lower) * choose(n - n1, t - i_lower) * choose(n-t, u-t))
  }
  umvue
}

#' Calculates the "uniformly minimal variance unbiased estimator".
#'
#' Calculates the "uniformly minimal variance unbiased estimator" (UMVUE) for the true response rate for the superset endpoint.
#'
#' @param t observed responses in the subset endpoint.
#' @param u observed responses in the superset endpoint.
#' @param r1 critical value for the first stage.
#' @param n1 sample size for the first stage.
#' @param n overall sample size.
#'
#'
#' @examples
#' #Setup "sub1"-object
#' sub1 <- setupSub1Design(pc0 = 0.5, pt0 = 0.6)
#'
#' #Calculate a subset design
#' design <- getSolutionsSub1(sub1, skipN1 = FALSE)$Solutions[4,]
#'
#' #Assume 9 responses in the subset endpoint and 13 responses in the superset endpoint were observed.
#' t = 9
#' u = 13
#' umvue_second <- get_UMVUE_GMS_subset_second_total(t, u, design$r1, design$n1, design$n)
get_UMVUE_GMS_subset_second_total <- function(t, u, r1, n1, n)
{
  stopifnot(class(t) == "numeric", t >= 0)
  stopifnot(class(u) == "numeric", u >= 0, u >= t)
  stopifnot(class(r1) == "numeric", r1 >= 0)
  stopifnot(class(n1) == "numeric", n1 >= 0, n1 > r1)
  stopifnot(class(n) == "numeric", n >= 0, n > n1)

  umvue <- get_UMVUE_GMS(t,r1,n1,n) + get_UMVUE_GMS_subset_second_only(t, u, r1, n1, n)
  umvue
}

#' Calculates the p-value.
#'
#' Calculates the p-value for a Simon's two-stage design based on the work from Koyama and Chen.
#'
#' @param k overall observed responses.
#' @param r1 critical value for the first stage.
#' @param n1 sample size for the first stage.
#' @param n overall sample size.
#' @param p0 response probability under the null hypothesis.
#'
#' @references Koyama T and Chen H (2008): Proper inference from Simon's two-stage designs. Statistics in Medicine, 27(16):3145-3154.
#'
#' @examples
#' #Calculate a Simon's two-stage design
#' design <- getSolutions()$Solutions[3,] #minimax-design for the default values.
#'
#' #Assume 9 responses were observed in the whole trial.
#' k = 9
#'
#' p_val <- get_p_KC(k, design$r1, design$n1, design$n, design$p0)
get_p_KC <-function(k, r1, n1, n, p0)
{
  stopifnot(class(k) == "numeric" | class(k) == "integer", k >= 0)
  stopifnot(class(r1) == "numeric", r1 >= 0)
  stopifnot(class(n1) == "numeric", n1 >= 0, n1 > r1)
  stopifnot(class(n) == "numeric", n >= 0, n > n1)
  stopifnot(class(p0) == "numeric", p0 >= 0, p0 <= 1)

  p <- numeric(0)
  if(k > r1)
  {
    k1 <- (r1 + 1) :  n1
    k2 <- k - k1
    n2 <- n - n1
    p <- sum(dbinom(k1, n1, p0)* (1 - pbinom(k2 -1, n2, p0)))
  }else
  {
    p <- 1 - pbinom(k-1, n1, p0)
  }

  p
}

#' Calculates the confidence interval.
#'
#' Calculates the two sided 1-2*alpha confidence interval based on the work from Koyama and Chen.
#'
#' @param k overall observed responses (must be larger than r1).
#' @param r1 critical value for the first stage.
#' @param n1 sample size for the first stage.
#' @param n overall sample size.
#' @param alpha determining the two sided 1-2*alpha confidence interval.
#' @param precision gives the precision (in decimal numbers) to which the confidence interval should be calculated (should be less than 10).
#'
#' @references Koyama T and Chen H (2008): Proper inference from Simon's two-stage designs. Statistics in Medicine, 27(16):3145-3154.
#'
#' @examples
#' #Calculate a Simon's two-stage design
#' design <- getSolutions()$Solutions[3,] #minimax-design for the default values.
#'
#' #Assume 9 responses were observed in the whole trial.
#' k = 9
#'
#' ci <- get_CI(k, design$r1, design$n1, design$n)
get_CI <-function(k, r1, n1, n, alpha = 0.05, precision = 4)
{
  stopifnot(class(k) == "numeric" | class(k) == "integer", k >= 0)
  stopifnot(class(r1) == "numeric", r1 >= 0)
  stopifnot(class(n1) == "numeric", n1 >= 0, n1 > r1)
  stopifnot(class(n) == "numeric", n >= 0, n > n1)
  stopifnot(class(alpha) == "numeric", alpha > 0, alpha <= 1)
  stopifnot(class(precision) == "numeric", precision >= 0, precision < 10)

  #could the study procceed to the second stage?
  stopifnot(k>r1, precision > 1)


  eff <- seq(0,1,0.01)
  tmp <- sapply(eff, function(bla) get_p_KC(k=k, r1 = r1, n1 = n1, n= n, p0 = bla))
  index_low <- which(tmp >= alpha)[1]
  index_high <- which(tmp >= (1-alpha))[1] - 1

  if(precision >2)
  {
    eff_low <- seq(eff[index_low - 1], eff[index_low], (1 / (10^(3))))
    eff_high <- seq(eff[index_high], eff[index_high + 1], (1 / (10^(3))))
    for(i in 2:(precision)) #each iteration provides +1 digits precision
    {
      if(i > 2)
      {
        eff_low <- seq(eff_low[index_low - 1], eff_low[index_low], (1 / (10^(i + 1))))
        eff_high <- seq(eff_high[index_high], eff_high[index_high + 1], (1 / (10^(i + 1))))
      }
      tmp_low <- sapply(eff_low, function(bla) get_p_KC(k=k, r1 = r1, n1 = n1, n= n, p0 = bla))
      tmp_high <- sapply(eff_high, function(bla) get_p_KC(k=k, r1 = r1, n1 = n1, n= n, p0 = bla))

      index_low <- which(tmp_low >= alpha)[1]
      index_high <- which(tmp_high >= (1- alpha))[1] - 1
    }

    result <- data.frame(CI_low = eff_low[index_low], CI_high = eff_high[index_high])

  }else
  {
    result <- data.frame(CI_low = eff[index_low], CI_high = eff[index_high])
  }
  result
}

#' Calculates the exact p value.
#'
#' Calculates the exact p value for a given subset design.
#'
#' @param t observed responses in the subset endpoint.
#' @param u observed responses in the superset endpoint.
#' @param r1 critical value for the first stage.
#' @param n1 sample size for the first stage.
#' @param n overall sample size.
#' @param pc0 the response probability for the subset endpoint under the null hypothesis.
#' @param pt0 the response probability for the superset endpoint under the null hypothesis.
#' @param sub1 "sub1"-object used to calculate the p value in c++ .
#'
#' @seealso \code{\link{setupSub1Design}}
#'
#' @examples
#' #Setup "sub1"-object
#' sub1 <- setupSub1Design(pc0 = 0.5, pt0 = 0.6)
#'
#' #Calculate a  subset design
#' design <- getSolutionsSub1(sub1, skipN1 = FALSE)$Solutions[4,]
#'
#' #Assuming 9 responses in the subset endpoint and 13 responses
#' #in the superset endpoint were observed.
#' t = 9
#' u = 13
#'
#' p_val <- get_p_exact_subset(t, u, design$r1, design$n1, design$n, design$pc0, design$pt0, sub1)
#' p_val
get_p_exact_subset <- function(t, u, r1, n1, n, pc0, pt0, sub1 = setupSub1Design())
{
  stopifnot(class(t) == "numeric", t >= 0)
  stopifnot(class(u) == "numeric", u >= 0, u >= t)
  stopifnot(class(r1) == "numeric", r1 >= 0)
  stopifnot(class(n1) == "numeric", n1 >= 0, n1 > r1)
  stopifnot(class(n) == "numeric", n >= 0, n > n1)
  stopifnot(class(pc0) == "numeric", pc0 > 0, pc0 <= 1)
  stopifnot(class(pt0) == "numeric", pt0 > 0, pt0 <= 1, pt0 >= pc0)
  stopifnot(class(sub1)[1] == "Rcpp_Sub1Design")

  p <- sub1$get_p_exact(t, u, r1, n1, n, pc0, pt0)
  p
}

#' Calculates the conditional power.
#'
#' Calculates the conditional power of a given subset design.
#'
#' @param t observed responses in the subset endpoint.
#' @param u observed responses in the superset endpoint.
#' @param enrolled number of patients enrolled so far.
#' @param r1 critical value for the first stage.
#' @param n1 sample size for the first stage.
#' @param r critical value for the subset endpoint.
#' @param s critical value for the superset endpoint.
#' @param n overall sample size.
#' @param pc1 the response probability under the alternative hypothesis for the subset endpoint.
#' @param pt1 the response probability under the alternative hypothesis for the superset endpoint.
#' @param sub1 "sub1"-object used to calculate the p value in c++ .
#'
#' @seealso \code{\link{setupSub1Design}}
#'
#' @examples
#' #Setup "sub1"-object
#' sub1 <- setupSub1Design(pc0 = 0.5, pt0 = 0.6)
#'
#' #Calculate a  subset design
#' design <- getSolutionsSub1(sub1, skipN1 = FALSE)$Solutions[4,]
#'
#' t <- 5
#' u <- 7
#' enrolled <- 10
#'
#' con_p <- get_conditionalPower(t, u, enrolled, design$r1,
#' design$n1, design$r, design$s, design$n, design$pc1, design$pt1, sub1)
get_conditionalPower <- function(t, u, enrolled, r1, n1, r, s, n, pc1, pt1, sub1 = setupSub1Design())
{
  stopifnot(class(t) == "numeric", t >= 0)
  stopifnot(class(u) == "numeric", u >= 0)
  stopifnot(class(enrolled) == "numeric", enrolled >= 0)
  stopifnot(class(r1) == "numeric", r1 >= 0)
  stopifnot(class(n1) == "numeric", n1 >= 0)
  stopifnot(class(r) == "numeric", r >= 0)
  stopifnot(class(s) == "numeric", s >= 0)
  stopifnot(class(n) == "numeric", n >= 0)
  stopifnot(class(pc1) == "numeric", pc1 > 0, pc1 <= 1)
  stopifnot(class(pt1) == "numeric", pt1 > 0, pt1 <= 1)
  stopifnot(class(sub1)[1] == "Rcpp_Sub1Design")
  stopifnot(t <= u, n1 < n , pc1 <= pt1, r1 < n1, r <= s, s < n)

  cp = sub1$get_conditionalPower(t, u, enrolled, r1, n1, r, s, n, pc1, pt1)
  cp
}

#' Calculates the confidence set.
#'
#' The p-value of  Subset Designs depends on two endpoints e.g. the superset endpoint and the subset endpoint.
#' Therefore the confidence interval for the response rate of the subset endpoint depends on the response rate of the superset endpoint and vice versa.
#' This results in a conficence "area" which is called the confidence set. "get_conficence_set" returns a set of points which outline the border of the confidence set.
#'
#' @param t observed responses in subset endpoint.
#' @param u observed responses in the superset endpoint.
#' @param r1 critical value for the first stage.
#' @param n1 sample size for the first stage.
#' @param n overall sample size.
#' @param pc0 the response probability under the null hypothesis for the subset endpoint.
#' @param pt0 the response probability under the null hypothesis for the superset endpoint.
#' @param alpha significance level the study was planned for.
#'
#' @references Kunz, C. U.  (2011), Two-stage designs for phase II trials with one or two endpoints. http://d-nb.info/1024218457
#'
#' @seealso \code{\link{setupSub1Design}}, \code{\link{plot_confidence_set}}
#'
#' @examples
#' #Setup "sub1"-object
#' sub1 <- setupSub1Design(pc0 = 0.5, pt0 = 0.6)
#'
#' #Calculate a  subset design
#' design <- getSolutionsSub1(sub1, skipN1 = FALSE)$Solutions[4,]
#'
#' t <- 12
#' u <- 13
#' alpha = 0.1
#'
#' conf_set <- get_confidence_set(t, u, design$r1, design$n1, design$n, design$pc0, design$pt0, alpha)
get_confidence_set <- function(t, u, r1, n1, n, pc0, pt0, alpha)
{
  stopifnot(class(t) == "numeric", t >= 0)
  stopifnot(class(u) == "numeric", u >= 0, u >= t)
  stopifnot(class(r1) == "numeric", r1 >= 0)
  stopifnot(class(n1) == "numeric", n1 >= 0, n1 > r1)
  stopifnot(class(n) == "numeric", n >= 0, n > n1)
  stopifnot(class(pc0) == "numeric", pc0 > 0, pc0 <= 1)
  stopifnot(class(pt0) == "numeric", pt0 > 0, pt0 <= 1, pt0 >= pc0)
  stopifnot(class(alpha) == "numeric", alpha > 0, alpha <= 1)

  set_boaders <- matrix(0, 20000, 2)
  sub1 <- setupSub1Design()

  for(pi1_lower in seq(pc0*0.8,1,0.01))
  {
    for(pi2_lower in seq(pi1_lower,1,0.01))
    {
      if(get_p_exact_subset(t,u,r1, n1, n, pi1_lower, pi2_lower, sub1) < alpha)
      {
        point_Index <- (pi1_lower * 100) * 100  + pi2_lower*100
        set_boaders[point_Index, ] <- c(pi1_lower, pi2_lower)
      }
    }
  }
  points <- set_boaders[(set_boaders[,1] > 0) | (set_boaders[,2] > 0), ]

  xValues <- as.numeric(names(table(points[,1])))
  numXValues <- length(xValues)
  x <- xValues[numXValues]
  y1 <- max(points[points[,1] == x, 2])


  area_boader <- data.frame(responserate_ep1 = x, responserate_ep2 = y1)

  for(l in seq(1, x, -0.01))
  {
    area_boader <- rbind(area_boader, data.frame(responserate_ep1 = l, responserate_ep2 = l))
  }

  numXValues <- numXValues -1
  while((y1 < 1) & (numXValues > 0))
  {
    x <- xValues[numXValues]
    index <- sum(points[,1] <= x)
    y1 <- points[index, 2]
    area_boader <- rbind(area_boader, data.frame(responserate_ep1 = x, responserate_ep2 = y1))
    numXValues <- numXValues -1
  }
  area_boader <- area_boader[-1,]
  area_boader
}

#' Plots the "confidence set" according to the observed responses.
#'
#' Plots the "confidence set" which can be received by invoking "get_confidence_set". Also the "uniformly minimal variance unbiased estimator" and the acceptance area are included in the plot.
#'
#' @param t observed responses in the subset endpoint.
#' @param u observed responses in the superset endpoint.
#' @param r1 critical value for the first stage.
#' @param n1 sample size for the first stage.
#' @param n overall sample size.
#' @param pc0 the response probability for the subset endpoint under the null hypothesis.
#' @param pt0 the response probability for the superset endpoint under the null hypothesis.
#' @param alpha overall significance level the trial was planned for.
#'
#' @references Kunz, C. U.  (2011), Two-stage designs for phase II trials with one or two endpoints. http://d-nb.info/1024218457
#'
#' @seealso \code{\link{get_confidence_set}}, \code{\link{get_UMVUE_GMS_subset_second_total}}, \code{\link{get_UMVUE_GMS}}
#'
#' @examples
#' #Setup "sub1"-object
#' sub1 <- setupSub1Design(pc0 = 0.5, pt0 = 0.6)
#'
#' #Calculate a  subset design
#' design <- getSolutionsSub1(sub1, skipN1 = FALSE)$Solutions[4,]
#'
#' #Assume 11 responses in the subset endpoint and 12 responses in the superset endpoint were observed.
#' t = 10
#' u = 12
#'
#' plot_confidence_set(t, u, design$r1, design$n1, design$n, design$pc0, design$pt0, 0.1)
plot_confidence_set <- function(t, u, r1, n1, n, pc0, pt0, alpha)
{
  stopifnot(class(t) == "numeric", t >= 0)
  stopifnot(class(u) == "numeric", u >= 0, u >= t)
  stopifnot(class(r1) == "numeric", r1 >= 0)
  stopifnot(class(n1) == "numeric", n1 >= 0, n1 > r1)
  stopifnot(class(n) == "numeric", n >= 0, n > n1)
  stopifnot(class(pc0) == "numeric", pc0 > 0, pc0 <= 1)
  stopifnot(class(pt0) == "numeric", pt0 > 0, pt0 <= 1, pt0 >= pc0)
  stopifnot(class(alpha) == "numeric", alpha > 0, alpha <= 1)

  plot(x=NA,y=NA, xlim=c(0,1), ylim=c(0,1), xlab = "Response Rate Subset Endpoint", ylab ="Response Rate Superset Endpoint")
  confSet <- get_confidence_set(t, u, r1, n1, n, pc0, pt0, alpha)
  polygon(confSet$responserate_ep1, confSet$responserate_ep2, col="darkgreen")
  polygon(c(0,0,pc0,pc0), c(0,pt0,pt0,0), col = "darkred")
  abline(h=pt0)
  abline(v=pc0)
  umvue_ep1 <- get_UMVUE_GMS(t, r1, n1, n)
  umvue_ep2 <- get_UMVUE_GMS_subset_second_total(t, u, r1, n1, n)
  points(umvue_ep1, umvue_ep2, pch = 19)
  legend("topleft", c("Confidence Set", "Acceptance Area", "UMVUE"), fill = c("darkgreen", "darkred", "black"))
}

#' Plots the study state of a given subset design.
#'
#' Plots the study state of a given subset design displaying the already enrolled patients and the stopping rules for the given study.
#'
#' @param sr dataframe containing the stopping rules for the given subset design defined by 3 columns named "Enrolled_patients", "Needed_responses_ep1" and "Needed_responses_ep2".
#'           This way each row defines when the study has to be stopped for futility.
#' @param enrolledPat dataframe defined by two boolean vectors named "ep1" and "ep2" indicating which patient had a response in the subset and superset endpoint.
#' @param r1 critical value for the first stage.
#' @param n1 sample size for the first stage.
#' @param r critical value for the subset endpoint.
#' @param s critical value for the superset endpoint.
#' @param n overall sample size.
#'
#' @seealso \code{\link{getSolutionsSub1}}
#'
#' @examples
#' #Calculate a subset design.
#' sub1 <- setupSub1Design(alpha = 0.1, beta = 0.2, pc0 = 0.3, pt0 = 0.4)
#' design <- getSolutionsSub1(sub1)$Solutions[10,]
#' #Define the stopping rules according to the chosen design.
#' sr <- data.frame(Enrolled_patients = c(design$n1, design$n),
#' Needed_responses_ep1 = c(design$r1, design$r), Needed_responses_ep2 = c(0,design$s))
#' #Simulate 14 random generated outcomes.
#' tmp_ep1 <- rbinom(14,1, design$pc1)
#' tmp_ep2 <- tmp_ep1 | rbinom(14,1, design$pt1)
#' enrolledPat <- data.frame(ep1 = tmp_ep1, ep2 = tmp_ep2)
#' #Plot study state.
#' plot_sub1_study_state(sr, enrolledPat, design$r1, design$n1, design$r, design$s, design$n)
plot_sub1_study_state <-function(sr, enrolledPat = data.frame(ep1 = logical(), ep2 = logical()), r1, n1, r, s, n)
{
  stopifnot(class(sr) == "data.frame")
  sr_names = names(sr)
  stopifnot("Enrolled_patients" %in% sr_names, "Needed_responses_ep1" %in% sr_names, "Needed_responses_ep2" %in% sr_names)
  stopifnot(class(enrolledPat) == "data.frame")
  ep_names = names(enrolledPat)
  stopifnot("ep1" %in% ep_names, "ep2" %in% ep_names)
  stopifnot(class(r1) == "numeric", r1 >= 0)
  stopifnot(class(n1) == "numeric", n1 >= 0)
  stopifnot(class(r) == "numeric", r >= 0)
  stopifnot(class(s) == "numeric", s >= 0)
  stopifnot(class(n) == "numeric", n >= 0)

  pats_present <- nrow(enrolledPat) > 0

  if(pats_present)
  {
    plot(x=NA, y=NA, xlim =c(0,n+1), ylim = c(0, s+ sum(enrolledPat$ep2)), xlab="Patients Enrolled", ylab="Observed Responses")
  }else{
    plot(x=NA, y=NA, xlim =c(0,n+1), ylim = c(0, s), xlab="Patients Enrolled", ylab="Observed Responses")
  }

  en <- sr$Enrolled_patients
  ep1 <- sr$Needed_responses_ep1
  ep2 <- sr$Needed_responses_ep2

  offset <- 0.2

  needed_ep2 <- as.numeric(names(table(ep2)))

  ep2_short <- numeric()
  rules_ep2 <- data.frame(en=en,ep2=ep2)

  for(i in needed_ep2){ep2_short <- c(ep2_short, rules_ep2$en[rules_ep2$ep2 == i][1])}
  rules_short <- cbind(enrolled = ep2_short, ep1_response = needed_ep2)
  numRows <- nrow(rules_short)

  plot_rules_ep2 <- data.frame(enrolled=numeric(2*numRows + 2), ep2_response=numeric(2*numRows + 2))
  for(i in 1:(numRows))
  {
    plot_rules_ep2[i*2 -1, ] <- c(rules_short[i,1], rules_short[i,2] +offset)
    if(i != numRows)
      plot_rules_ep2[i*2,] <- cbind(rules_short[i+1,1], rules_short[i,2] +offset)
    else
      plot_rules_ep2[i*2,] <- cbind(rules_short[i,1] + 1, rules_short[i,2] +offset)
  }

  plot_rules_ep2 <- rbind(data.frame(enrolled=plot_rules_ep2$enrolled[1], ep2_response=0), plot_rules_ep2)
  plot_rules_ep2[2*numRows +1,] <- c(n + 1, plot_rules_ep2$ep2_response[2*numRows])
  plot_rules_ep2[2*numRows +2,] <- c(n + 1,0)
  polygon(plot_rules_ep2$enrolled, plot_rules_ep2$ep2_response, col="lightblue")

  needed_ep1 <- as.numeric(names(table(ep1)))

  ep1_short <- numeric()
  rules <- data.frame(en=en,ep1=ep1)

  for(i in needed_ep1){ep1_short <- c(ep1_short, rules$en[rules$ep1 == i][1])}
  rules_short <- cbind(enrolled = ep1_short, ep1_response = needed_ep1)
  numRows <- nrow(rules_short)

  plot_rules <- data.frame(enrolled=numeric(2*numRows + 2), ep1_response=numeric(2*numRows + 2))
  for(i in 1:(numRows))
  {
    plot_rules[i*2 -1, ] <- c(rules_short[i,1], rules_short[i,2] +offset)
    if(i != numRows)
      plot_rules[i*2,] <- cbind(rules_short[i+1,1], rules_short[i,2] +offset)
    else
      plot_rules[i*2,] <- cbind(rules_short[i,1] + 1, rules_short[i,2] +offset)
  }

  plot_rules <- rbind(data.frame(enrolled=plot_rules$enrolled[1], ep1_response=0), plot_rules)
  plot_rules[2*numRows +1,] <- c(n + 1, plot_rules$ep1_response[2*numRows])
  plot_rules[2*numRows +2,] <- c(n + 1,0)
  polygon(plot_rules$enrolled, plot_rules$ep1_response, col="red")
  abline(v = c(n1, n), lty = 3, col = "darkred")
  abline(h = c(r1,r), lty = 3, col = "green")
  abline(h = s, lty=3, col = "darkblue")

  if(pats_present)
  {
    legend_items <- c("r1, r", "s", "n1, n", "Stopping Rules Ep1", "Stopping Rules Ep2", "Response Ep1", "Response Ep2")
    legend_col <- c("green", "darkblue", "darkred", "red", "blue", NA, NA)
    legend_pch <- c(NA,NA,NA,NA,NA,"o ","x ")
  }else{
    legend_items <- c("r1, r", "s", "n1, n", "Stopping Rules Ep1", "Stopping Rules Ep2")
    legend_col <- c("green", "darkblue", "darkred", "red", "lightblue")
    legend_pch <- c(NA,NA,NA,NA,NA)
  }

  legend("topleft", legend=legend_items , fill = legend_col, pch = legend_pch, border = "white")

  #plot enrolled patients

  if(pats_present)
  {
    numRes <- nrow(enrolledPat)
    responses <- data.frame(ep1 = numeric(numRes), ep2 = numeric(numRes), enrolled = (1:numRes))
    ep1_count <- 0
    ep2_count <- 0
    for(i in 1:numRes)
    {
      ep1_count <- ep1_count + enrolledPat$ep1[i]
      ep2_count <- ep2_count + enrolledPat$ep2[i]
      responses[i,1:2] <- c(ep1_count, ep2_count)
    }
    points(responses$enrolled, responses$ep1, type="b")
    points(responses$enrolled, responses$ep2, pch="x", type="b")
  }
}

#' Plots the study state of a given Simon's two-stage design.
#'
#' Plots the study state of a given Simon's two-stage design displaying the already enrolled patients and the stopping rules.
#'
#' @param sr dataframe containing the stopping rules for the given Simon's two-stage design defined by two columns named "Enrolled_patients" and "Needed_responses_ep1".
#'           This way each row defines when the study has to be stopped for futility.
#' @param enrolledPat dataframe defined by a boolean vector in one column named "ep1" indicating which patient had a response.
#' @param r1 critical value for the first stage.
#' @param n1 sample size for the first stage.
#' @param r critical value for the subset endpoint.
#' @param n overall sample size.
#'
#' @examples
#' #Calculate a Simon's two-stage design
#' design <- getSolutions()$Solutions[3,] #minimax-design for the default values.
#' #Define the stopping rules according to the chosen design
#' sr <- data.frame(Enrolled_patients = c(design$n1, design$n),
#' Needed_responses_ep1 = c(design$r1, design$r))
#' #Simulate 18 random generated outcomes.
#' enrolledPat <- data.frame(ep1 = rbinom(18,1, design$p1))
#' #Plot study state
#' plot_simon_study_state(sr, enrolledPat, design$r1, design$n1, design$r, design$n)
plot_simon_study_state <-function(sr, enrolledPat = data.frame(ep1 = logical()), r1, n1, r, n)
{
  stopifnot(class(sr) == "data.frame")
  sr_names = names(sr)
  stopifnot("Enrolled_patients" %in% sr_names, "Needed_responses_ep1" %in% sr_names)
  stopifnot(class(enrolledPat) == "data.frame")
  ep_names = names(enrolledPat)
  stopifnot("ep1" %in% ep_names)
  stopifnot(class(r1) == "numeric", r1 >= 0)
  stopifnot(class(n1) == "numeric", n1 >= 0)
  stopifnot(class(r) == "numeric", r >= 0, r > r1)
  stopifnot(class(n) == "numeric", n >= 0, n > n1)

  pats_present <- nrow(enrolledPat) > 0

  if(pats_present)
  {
    plot(x=NA, y=NA, xlim =c(0,n+1), ylim = c(0, r + sum(enrolledPat$ep1)), xlab="Patients Enrolled", ylab="Observed Responses")
  }else{
    plot(x=NA, y=NA, xlim =c(0,n+1), ylim = c(0, r), xlab="Patients Enrolled", ylab="Observed Responses")
  }

  en <- sr$Enrolled_patients
  ep1 <- sr$Needed_responses_ep1

  offset <- 0.2

  needed_ep1 <- as.numeric(names(table(ep1)))

  ep1_short <- numeric()
  rules <- data.frame(en=en,ep1=ep1)

  for(i in needed_ep1){ep1_short <- c(ep1_short, rules$en[rules$ep1 == i][1])}
  rules_short <- cbind(enrolled = ep1_short, ep1_response = needed_ep1)
  numRows <- nrow(rules_short)

  plot_rules <- data.frame(enrolled=numeric(2*numRows + 2), ep1_response=numeric(2*numRows + 2))
  for(i in 1:(numRows))
  {
    plot_rules[i*2 -1, ] <- c(rules_short[i,1], rules_short[i,2] +offset)
    if(i != numRows)
      plot_rules[i*2,] <- cbind(rules_short[i+1,1], rules_short[i,2] +offset)
    else
      plot_rules[i*2,] <- cbind(rules_short[i,1] + 1, rules_short[i,2] +offset)
  }

  plot_rules <- rbind(data.frame(enrolled=plot_rules$enrolled[1], ep1_response=0), plot_rules)
  plot_rules[2*numRows +1,] <- c(n + 1, plot_rules$ep1_response[2*numRows])
  plot_rules[2*numRows +2,] <- c(n + 1,0)
  polygon(plot_rules$enrolled, plot_rules$ep1_response, col="red")
  abline(v = c(n1, n), lty = 3, col = "blue")
  abline(h = c(r1,r), lty = 3, col = "green")

  if(pats_present)
  {
    legend_items <- c("r1, r", "n1, n", "Stopping Rules", "Observed Response")
    legend_col <- c("green", "blue", "red", "black")
  }else{
    legend_items <- c("r1, r", "n1, n", "Stopping Rules")
    legend_col <- c("green", "blue", "red")
  }

  legend("topleft", legend=legend_items , fill = legend_col, border = "white")

  #plot enrolled patients

  if(pats_present)
  {
    numRes <- nrow(enrolledPat)
    responses <- data.frame(ep1 = numeric(numRes), ep2 = numeric(numRes), enrolled = (1:numRes))
    ep1_count <- 0
    ep2_count <- 0
    for(i in 1:numRes)
    {
      ep1_count <- ep1_count + enrolledPat$ep1[i]
      responses[i,1:2] <- c(ep1_count, ep2_count)
    }
    points(responses$enrolled, responses$ep1, type="b")
  }
}

#'Returns the conditional power.
#'
#'Returns the conditional power when "k" responses where observed out of "numPat" Patients for the given Simon's tow stage design.
#'
#' @param k number of observed responses
#' @param numPat number of enrolled patients.
#' @param r1 critical value for the first stage.
#' @param n1 sample size for the first stage.
#' @param r critical value for the subset endpoint.
#' @param n overall sample size.
#' @param p1 response rate under the alternative hypothesis.
#'
#' @examples
#' #Calculate a Simon's two-stage design
#' design <- getSolutions()$Solutions[3,] #minimax-design for the default values.
#' #Assume 3 out of 20 patients had a response.
#' getCP_simon(3,20,design$r1, design$n1, design$r, design$n, design$p1)
getCP_simon <-function(k, numPat, r1, n1, r, n, p1)
{
  stopifnot(class(k) == "numeric" | class(k) == "integer", k >= 0)
  stopifnot(class(numPat) == "numeric", numPat >= 0, numPat >= k)
  stopifnot(class(r1) == "numeric", r1 >= 0)
  stopifnot(class(n1) == "numeric", n1 >= 0)
  stopifnot(class(r) == "numeric", r >= 0)
  stopifnot(class(n) == "numeric", n >= 0)
  stopifnot(class(p1) == "numeric", p1 > 0, p1 <= 1)

  simon <- setupSimon()
  cp <- simon$getConditionalPower(k, numPat, r1, n1, r, n, p1)
  cp
}