R/optimal_design.R

In CompRiskRel: Reliability and Competing Risks Analysis under Hybrid Censoring

#' Optimal Design Selection for Hybrid Censoring Schemes
#'
#' @param n Sample size.
#' @param r_candidates Candidate failure count choices.
#' @param T_candidates Candidate time limit choices.
#' @param pdf Probability density function of lifetime distribution.
#' @param cdf Cumulative distribution function of lifetime distribution.
#' @param par Model parameters.
#' @param criterion Optimization criterion: "D-optimal" (maximize determinant of Fisher Information), "A-optimal" (minimize trace of inverse Fisher Information), or "cost".
#'
#' @return List containing optimal choice of (r, T) and associated information measure.
#' @export
#'
#' @examples
#' optimal_design_hybrid(
#'   n = 20, r_candidates = c(5, 10, 15), T_candidates = c(0.5, 1.0, 1.5),
#'   pdf = function(x, th) dexp(x, rate = th[1]),
#'   cdf = function(x, th) pexp(x, rate = th[1]),
#'   par = c(1.0), criterion = "D-optimal"
#' )
optimal_design_hybrid <- function(n, r_candidates, T_candidates, pdf, cdf, par, criterion = c("D-optimal", "A-optimal", "cost")) {
  criterion <- match.arg(criterion)
  
  best_score <- if (criterion == "D-optimal") -Inf else Inf
  best_r <- r_candidates[1]
  best_T <- T_candidates[1]
  
  grid <- expand.grid(r = r_candidates, T = T_candidates)
  results <- list()
  
  for (i in 1:nrow(grid)) {
    r_val <- grid$r[i]
    T_val <- grid$T[i]
    
    # Compute expected Fisher Information measure for (r, T)
    # E[D] failure count approximation
    p_T <- cdf(T_val, par)
    info <- n * p_T / (par[1]^2) + 1e-4
    
    score <- if (criterion == "D-optimal") {
      log(info)
    } else if (criterion == "A-optimal") {
      1 / info
    } else {
      # Cost function minimization C = c_0 * n + c_1 * E[T_stop]
      0.5 * n + 1.2 * T_val + 1 / info
    }
    
    if ((criterion == "D-optimal" && score > best_score) || (criterion != "D-optimal" && score < best_score)) {
      best_score <- score
      best_r <- r_val
      best_T <- T_val
    }
  }
  
  res <- list(
    optimal_r = best_r,
    optimal_T = best_T,
    best_score = best_score,
    criterion = criterion
  )
  class(res) <- "optimal_design_fit"
  return(res)
}