R/Solve.beta.given.power.R

In gsMeanFreq: Group Sequential Clinical Trial Designs for Composite Endpoints

#' @title Function to solve for the treatment effect parameter to achieve target power by Monte Carlo simulation.
#' @description
#' For each design scenario in \code{params}, this function solves for the treatment effect coefficient \code{beta.trt} that achieves the
#' desired power using an iterative Monte Carlo calibration procedure. For scenarios labeled \code{TypeI}, the function sets \code{beta.trt = 0}.
#' For scenarios labeled \code{Power}, it repeatedly simulates two-sample composites data, estimates calibration quantities (a power constant and
#' variance) using either the generalized log-rank (\code{"LR"}) or generalized-t (\code{"GT"}) approach, updtaes \code{beta.trt} using
#' \code{find.beta.trt()}, and iterates until convergence within \code{tol}.
#'
#'
#' @param nsim Integer giving the number of Monte Carlo replicates used for each iteration. Default is \code{1000}.
#' @param params A data frame where each rows defines a simulation/design scenario. Must include a column \code{Type} with
#' values \code{"TypeI"} or \code{"Power"}. For \code{"Power"} scenarios, the function expects the following columns:
#' \code{placebo.lambda_0}, \code{grp.size}, \code{recruitment},
#' and \code{random.censor.rate}. Additional columns are carried through to the output.
#' @param estimator Character string specifying which calibration method to use:
#'   \code{"LR"} for the generalized log-rank statistic or \code{"GT"} for the
#'   generalized-t statistic.
#' @param tol Positive numeric value giving the convergence tolerance for the
#'   fixed-point iteration in \code{beta.trt}. The algorithm stops when
#'   \code{abs(beta.trt.old - beta.trt.new) < tol}. Default is \code{0.001}.
#' @param seed Seed for reproducibility.
#'
#' @returns  A data.frame with the same rows as \code{params} and an additional
#' column \code{beta.trt} containing the solved treatment effect coefficient.
#' For \code{Type == "TypeI"}, \code{beta.trt} is set to \code{0}. For
#' \code{Type == "Power"}, \code{beta.trt} is the converged solution from the
#' Monte Carlo calibration procedure.
#' @export
#' @importFrom dplyr bind_cols bind_rows
#' @importFrom foreach foreach %dopar%
#'
Solve.beta.given.power <- function(nsim = 1000, params, estimator, tol = 0.001, seed = NULL){
  # Handle the seed dynamically
  if (!is.null(seed)) set.seed(seed)
  base_seed <- sample.int(1e6, 1)


  results <- data.frame()
  for (j in 1:nrow(params)){
    param <- params[j,]
    if (param$Type == "TypeI"){
      temp <- bind_cols(param, beta.trt = 0)
      results <- bind_rows(results, temp)
    } else if (param$Type == "Power") {

      # for each combination, estimates the constants and variance from Monte Carlo simulations

      beta.trt.old = 0
      # beta.iterations <- c(beta.trt.old)  # <- Track beta values

      # Initial estimates: solve the first beta.trt under the null
      estimates <- foreach(i = 1*(1:nsim), .combine = rbind, .packages = c("pracma", "bdsmatrix"),
                           .errorhandling = "remove") %dopar%{
                             df <- TwoSample.generate.sequential(lambda_0vec = rep(param$placebo.lambda_0, 2),
                                                                 sizevec = rep(param$grp.size, 2), beta.trt = beta.trt.old,
                                                                 recruitment = param$recruitment,
                                                                 random.censor.rate = param$random.censor.rate,
                                                                 seed = base_seed+i)
                             if (estimator == 'LR'){
                               est = TwoSample.Constant.LR(data = df)
                             } else if (estimator == 'GT'){
                               est = TwoSample.Constant.GT(data = df)
                             }
                             # est <- estimator.fn(data = df)
                             est
                           }

      beta.trt.new <- find.beta.trt(power = 0.8, power.c = mean(estimates[,1]), var = mean(estimates[,2]))
      # beta.iterations <- c(beta.iterations, beta.trt.new)

      # use beta.trt to get the updated constants and variance until converge

      while (abs(beta.trt.old - beta.trt.new) >= tol) {
        beta.trt.old <- beta.trt.new
        estimates <- foreach(i = 1*(1:nsim), .combine = rbind, .packages = c("pracma", "bdsmatrix"),
                             .errorhandling = "remove") %dopar%{
                               df <- TwoSample.generate.sequential(lambda_0vec = rep(param$placebo.lambda_0, 2),
                                                                   sizevec = rep(param$grp.size, 2), beta.trt = beta.trt.old,
                                                                   recruitment = param$recruitment,
                                                                   random.censor.rate = param$random.censor.rate,
                                                                   seed = base_seed+i)
                               if (estimator == 'LR'){
                                 est = TwoSample.Constant.LR(data = df)
                               } else if (estimator == 'GT'){
                                 est = TwoSample.Constant.GT(data = df)
                               }
                               # est <- estimator.fn(data = df)
                               est
                             }

        beta.trt.new <- find.beta.trt(power = 0.8, power.c = mean(estimates[,1]), var = mean(estimates[,2]))
        # beta.iterations <- c(beta.iterations, beta.trt.new)  # <- Append
      }
      temp <- bind_cols(param, beta.trt = beta.trt.new)
      # temp$beta.path <- list(beta.iterations)  # <- Save full path
      results <- bind_rows(results, temp)
    }
  } # end of the "j" loop
  return(results)
}