R/getPPKinits.R

In nlmixr2autoinit: Automatic Generation of Initial Estimates for Population Pharmacokinetic Modeling

#' Automated pipeline for generating initial estimates in population PK models
#'
#' Provides a unified and fully automated workflow to generate initial
#' pharmacokinetic and residual variability parameters for population PK models
#' using concentration–time data from bolus, infusion, or oral administration.
#'
#' @param dat A data frame containing pharmacokinetic records in standard nlmixr2
#'   format, including ID, TIME, EVID, and DV.
#' @param control A list created by `initsControl()` specifying configuration for
#'   pooling, non-compartmental analysis, steady-state detection, fallback rules,
#'   statistical model components, and parameter selection metrics.
#' @param verbose Logical (default = TRUE); when TRUE, displays key progress
#'   messages and stepwise updates during the initialization process. When FALSE,
#'   the function runs quietly without printing intermediate information.
#'
#' @details
#' The pipeline integrates four model-informed analytical components applied to
#' raw or pooled concentration–time profiles:
#' \enumerate{
#'   \item Adaptive single-point methods
#'   \item Naive pooled graphic methods
#'   \item Naive pooled non-compartmental analysis (NCA) with optional Wagner–Nelson Ka calculation for oral dosing
#'   \item Parameter sweeping across one-, two-, three-compartment and Michaelis–Menten models
#' }
#'
#' In addition to structural PK parameters, the framework also initializes
#' statistical model components:
#' \itemize{
#'   \item Inter-individual variability (IIV): pragmatic fixed \eqn{\omega^2}
#'     values are assigned to random effects.
#'   \item Residual unexplained variability (RUV): estimated either using a
#'     data-driven method based on trimmed residual summaries or a
#'     fixed-fraction approach consistent with NONMEM User Guide recommendations.
#'   \item Model applicability: the automated and model-informed strategy
#'     generates robust initial values suitable for both linear and nonlinear
#'     mixed-effects pharmacokinetic models.
#' }
#'
#' @return An object of class `getPPKinits` containing recommended initial
#' parameter estimates, intermediate results, and computation diagnostics.
#'
#' @seealso \code{\link{initsControl}}, \code{\link{run_single_point}},
#'   \code{\link{run_graphcal}}, \code{\link{run_pooled_nca}},
#'   \code{\link{sim_sens_1cmpt_mm}}, \code{\link{sim_sens_2cmpt}},
#'   \code{\link{sim_sens_3cmpt}}, \code{\link{metrics.}}
#'
#' @author Zhonghui Huang
#'
#' @examples
#' \donttest{
#' ## Bolus example
#' getPPKinits(Bolus_1CPT,verbose = TRUE)
#' ## Oral example (run quietly)
#' getPPKinits(Oral_1CPT,verbose = FALSE)
#' }
#'
#' @export

getPPKinits <- function(dat,
                        control = initsControl(),
                        verbose = TRUE) {

  # load control
  nlmixr2est::nlsControl()
  nlmixr2est::nlmControl()
  nlmixr2est::nlminbControl()
  nlmixr2est::saemControl()
  nlmixr2est::foceiControl()

  `%>%` <- magrittr::`%>%`

  # Unpack control components
  .pooledctrl    <- control$pooled.control
  .ncactrl       <- control$nca.control
  .ssctrl        <-  control$ss.control
  .fbctrl        <- control$fallback.control
   selmetrics    <- control$selmetrics
   hybrid.base   <- control$hybrid.base
   preferNCA     <- control$preferNCA

  ################# 1. Data preprocessing #################
  # Record start time
  start.time <- Sys.time()
  process_result <- processData(dat,verbose=verbose)
  dat <- process_result$dat
  dose_type <-
    process_result$Datainfo$Value[process_result$Datainfo$Infometrics == "Dose Type"]
  route <-
  process_result$Datainfo$Value[process_result$Datainfo$Infometrics == "Dose Route"]
  Datainfo<-process_result$Datainfo

  # Reset ID
  dat <- dat %>%
    dplyr::mutate(ID = ID + (resetflag - 1) * max(dat$ID, na.rm = TRUE))

  # Adjusting CMT column name for internal model calculation
  if (route == "oral") {
    # Define expected CMT based on EVID
    expected_cmt <- ifelse(dat$EVID == 1, "depot",
                           ifelse(dat$EVID == 0, "centre", dat$CMT))

    # Check if any value is different
    if (any(dat$CMT != expected_cmt)) {
      dat$CMT <- expected_cmt
    }
  } else {
    # For bolus or infusion: all CMT should be "centre"
    if (any(dat$CMT != "centre")) {
      dat$CMT <- "centre"
    }
  }

  # obtain the pooled data
  pooled_data <- get_pooled_data(dat,
                                 dose_type,
                                 .pooledctrl)

  # estimate the half-life
  half_life_out <- get_hf(dat = dat,pooled = pooled_data,verbose = verbose)
  half_life <-half_life_out$half_life_median
  ##################2. Single-point method  #################
  sp_result <- run_single_point(
    dat = dat,
    route = route,
    half_life = half_life,
    dose_type = dose_type,
    pooled_ctrl = .pooledctrl,
    ssctrl = .ssctrl
  )

  dat <- sp_result$dat
  sp_out <- sp_result$singlepoint.results

  used_sp_ka_fallback  <- FALSE

  if (.fbctrl$enable_ka_fallback &&
      route == "oral" &&
      (is.na(sp_out$ka) || sp_out$ka < 0)){
      sp_out$ka <- 1

   used_sp_ka_fallback<- TRUE
  }

  ####################3. Graphical residuals  ##################
  graph_out <- run_graphcal(
    dat,
    route,
    data_type,
    pooled =  pooled_data,
    pooled_ctrl = .pooledctrl,
    nlastpoints = 3
  )

  used_graph_ka_fallback  <- FALSE
  if (.fbctrl$enable_ka_fallback &&
      route == "oral" &&
      (is.na(graph_out$ka) || graph_out$ka < 0)) {
       graph_out$ka <- 1
       used_graph_ka_fallback  <- TRUE
  }


  ###################4. Naive pooled NCA   ##################
  nca_out <- run_pooled_nca(dat,
                            dose_type=dose_type,
                            pooled=pooled_data,
                            route=route,
                            pooled_ctrl=pooled_control(),
                            nca_ctrl=nca_control()
                            )
  ka_nca_fd <- NA
  ka_nca_efd <- NA
  ka_nca_all <- NA
  ka_wn <- NA
  used_nca_ka_fallback  <- FALSE

  if (route == "oral" &&
      !is.na(nca_out$nca.fd.results$clobs) &&
      !is.na(nca_out$nca.fd.results$vzobs)) {
    ka_wn <- ka_wanger_nelson(pooled_data$datpooled_fd$binned.df,
                              nca_out$nca.fd.results)
    ka_nca_fd <- signif(ka_wn$ka, 3)

    if (.fbctrl$enable_ka_fallback) {
      if ((is.na(ka_nca_fd) || ka_nca_fd < 0)) {
        ka_nca_fd <- 1
        used_nca_ka_fallback <- TRUE
      }
      ka_nca_efd <- 1
      ka_nca_all <- 1
    }
  }

  ###############5. Base parameter predictive performance   ##################
  # This section performs model evaluation using parameters from each method separately,
  # without mixing ka, CL, and Vd across methods (i.e., not hybrid).
  if (!hybrid.base){
    if (verbose){
    message(crayon::black(
      paste0(
        "Evaluating the predictive performance of calculated one-compartment model parameters",
        strrep(".", 20)
      )
    ))
    }
  base_perf.out <- list(
    simpcal = eval_perf_1cmpt(dat, "rxSolve", sp_out$ka, sp_out$cl, sp_out$vd, route),
    graph   = eval_perf_1cmpt(dat, "rxSolve", graph_out$ka, graph_out$cl, graph_out$vd, route),
    nca_fd  = eval_perf_1cmpt(
      dat,
      "rxSolve",
      ka_nca_fd,
      nca_out$nca.fd.results$clobs,
      nca_out$nca.fd.results$vzobs,
      route
    ),
    nca_efd = eval_perf_1cmpt(
      dat,
      "rxSolve",
      ka_nca_efd,
      nca_out$nca.efd.results$clobs,
      nca_out$nca.efd.results$vzobs,
      route
    ),
    nca_all = eval_perf_1cmpt(
      dat,
      "rxSolve",
      ka_nca_all,
      nca_out$nca.all.results$clobs,
      nca_out$nca.all.results$vzobs,
      route
    )
  )
  ka <- c(NA,NA,NA,NA,NA)

  if (route=="oral"){
    ka <- c(
      sp_out$ka,
      graph_out$ka,
      ka_nca_fd,     # by wanger nelson
      ka_nca_efd,         # no nca for md
      ka_nca_all         # no nca for md
    )
  }

  base.out <- data.frame(
    method = c(
      "Adaptive single-point method",
      "Graphic methods",
      "Naive pooled NCA (FD)",
      "Naive pooled NCA (MD)",
      "Naive pooled NCA (all)"
    ),

    ka = ka,

    cl = c(
      sp_out$cl,
      graph_out$cl,
      nca_out$nca.fd.results$clobs,
      nca_out$nca.efd.results$clobs,
      nca_out$nca.all.results$clobs
    ),

    vd = c(
      sp_out$vd,
      graph_out$vd,
      nca_out$nca.fd.results$vzobs,
      nca_out$nca.efd.results$vzobs,
      nca_out$nca.all.results$vzobs
    ),

    simAPE = c(
      unname(base_perf.out$simpcal["APE"]),
      unname(base_perf.out$graph["APE"]),
      unname(base_perf.out$nca_fd["APE"]),
      unname(base_perf.out$nca_efd["APE"]),
      unname(base_perf.out$nca_all["APE"])
    ),

    simMAE = c(
      unname(base_perf.out$simpcal["MAE"]),
      unname(base_perf.out$graph["MAE"]),
      unname(base_perf.out$nca_fd["MAE"]),
      unname(base_perf.out$nca_efd["MAE"]),
      unname(base_perf.out$nca_all["MAE"])
    ),

    simMAPE = c(
      unname(base_perf.out$simpcal["MAPE"]),
      unname(base_perf.out$graph["MAPE"]),
      unname(base_perf.out$nca_fd["MAPE"]),
      unname(base_perf.out$nca_efd["MAPE"]),
      unname(base_perf.out$nca_all["MAPE"])
    ),

    simRMSE = c(
      unname(base_perf.out$simpcal["RMSE"]),
      unname(base_perf.out$graph["RMSE"]),
      unname(base_perf.out$nca_fd["RMSE"]),
      unname(base_perf.out$nca_efd["RMSE"]),
      unname(base_perf.out$nca_all["RMSE"])
    ),

    simrRMSE1 = c(
      unname(base_perf.out$simpcal["rRMSE1"]),
      unname(base_perf.out$graph["rRMSE1"]),
      unname(base_perf.out$nca_fd["rRMSE1"]),
      unname(base_perf.out$nca_efd["rRMSE1"]),
      unname(base_perf.out$nca_all["rRMSE1"])
    ),

    simrRMSE2 = c(
      unname(base_perf.out$simpcal["rRMSE2"]),
      unname(base_perf.out$graph["rRMSE2"]),
      unname(base_perf.out$nca_fd["rRMSE2"]),
      unname(base_perf.out$nca_efd["rRMSE2"]),
      unname(base_perf.out$nca_all["rRMSE2"])
    ),

    time.spent = c(
      sp_out$time.spent,
      graph_out$time.spent,
      nca_out$nca.fd.results$time.spent,
      nca_out$nca.efd.results$time.spent,
      nca_out$nca.all.results$time.spent
    )
  )

  colnames(base.out) <- c(
    "Method",
    "Calculated Ka",
    "Calculated CL",
    "Calculated Vd",
    "Absolute Predicted Error (APE)",
    "Mean Absolute Error (MAE)",
    "Mean Absolute Percentage Error (MAPE)",
    "Root Mean Squared Error (RMSE)",
    "Relative Root Mean Squared Error (rRMSE1)",
    "Relative Root Mean Squared Error (rRMSE2)",
    "Time spent"
  )
  }

  # Hybrid mode: evaluates predictive performance using mixed-source parameters.
  if (hybrid.base){
    if (verbose){
    message(crayon::black(
      paste0(
        "Evaluating the predictive performance of calculated one-compartment model parameters",
        strrep(".", 20)
      )
    ))
    }
    base.out <- hybrid_eval_perf_1cmpt(
      route = route,
      dat = dat,
      sp_out_ka     = sp_out$ka,
      sp_out_cl     = sp_out$cl,
      sp_out_vd     = sp_out$vd,
      graph_out_ka  = graph_out$ka,
      graph_out_cl  = graph_out$cl,
      graph_out_vd  = graph_out$vd,
      nca_fd_ka     = ka_nca_fd,
      nca_fd_cl     = nca_out$nca.fd.results$clobs,
      nca_fd_vd     = nca_out$nca.fd.results$vzobs,
      nca_efd_ka    = ka_nca_efd,
      nca_efd_cl    = nca_out$nca.efd.results$clobs,
      nca_efd_vd    = nca_out$nca.efd.results$vzobs,
      nca_all_ka    = ka_nca_all,
      nca_all_cl    = nca_out$nca.all.results$clobs,
      nca_all_vd    = nca_out$nca.all.results$vzobs,
      verbose = verbose
    )

    colnames(base.out) <- c(
      "Ka Method",
      "CL Method",
      "Vd Method",
      "Calculated Ka",
      "Calculated CL",
      "Calculated Vd",
      "Absolute Predicted Error (APE)",
      "Mean Absolute Error (MAE)",
      "Mean Absolute Percentage Error (MAPE)",
      "Root Mean Squared Error (RMSE)",
      "Relative Root Mean Squared Error (rRMSE1)",
      "Relative Root Mean Squared Error (rRMSE2)"
    )
  }

  cols_to_format <-
    c("Calculated Ka", "Calculated CL", "Calculated Vd")
  base.out[cols_to_format] <-
    lapply(base.out[cols_to_format], function(x)
      signif(as.numeric(x), digits = 3))

#################7. Select the base parameters by metrics####################

  # Full list of available performance metric column names
  metrics.all <- c(
    "Absolute Predicted Error (APE)",
    "Mean Absolute Error (MAE)",
    "Mean Absolute Percentage Error (MAPE)",
    "Root Mean Squared Error (RMSE)",
    "Relative Root Mean Squared Error (rRMSE1)",
    "Relative Root Mean Squared Error (rRMSE2)"
  )

  # Extract abbreviated metric keys from column names (e.g., "APE", "MAE", etc.)
  metric_keys <- sub(".*\\((.*)\\)", "\\1", metrics.all)

  # Identify selected metric columns based on user-specified selmetrics
  stat_cols  <- metrics.all[metric_keys %in% selmetrics]

  # Extract the column-wise minimum values for each performance metric
  mins <- sapply(base.out[stat_cols], min, na.rm = TRUE)
  lmat <- sweep(base.out[stat_cols], 2, mins, FUN = "==")
  base.out$min_count <- rowSums(lmat, na.rm = TRUE)

  if (length(stat_cols)==1){
    base.out$metrics.rank = rank(base.out[[stat_cols]],)
  }

  base.best <- base.out %>%
    dplyr::arrange(
      `Relative Root Mean Squared Error (rRMSE2)`,
      `Mean Absolute Percentage Error (MAPE)`
    ) %>%
    dplyr::filter(min_count == max(min_count)) %>%
    dplyr::slice(1)

  # Defensive check: if selected model has rRMSE1 > 10x the global minimum, fallback to rRMSE1 + MAPE
  rRMSE1_selected <-
    base.best$`Relative Root Mean Squared Error (rRMSE1)`
  rRMSE1_min <-
    min(base.out$`Relative Root Mean Squared Error (rRMSE1)`, na.rm = TRUE)

  # Only apply defensive fallback when rRMSE2 is the only selected metric
  if (identical(selmetrics, "rRMSE2") &&
      rRMSE1_selected > 10 * rRMSE1_min) {
    # Define fallback metric set
    fallback_metrics <- c("Relative Root Mean Squared Error (rRMSE1)")

    # Recalculate min_count using fallback metrics
    fallback_mins <-
      sapply(base.out[fallback_metrics], min, na.rm = TRUE)
    fallback_lmat <-
      sweep(base.out[fallback_metrics], 2, fallback_mins, FUN = "==")
    base.out$min_count <- rowSums(fallback_lmat, na.rm = TRUE)

    # Re-select best model using fallback metrics
    base.best <- base.out %>%
      dplyr::arrange(`Relative Root Mean Squared Error (rRMSE1)`,
                     `Mean Absolute Percentage Error (MAPE)`) %>%
      dplyr::filter(min_count == max(min_count)) %>%
      dplyr::slice(1)
  }

  # Prefer NCA
  if (preferNCA) {
    if (identical(selmetrics, "rRMSE2")) {
    rrmse2_best <- base.out %>%
      dplyr::arrange(`Relative Root Mean Squared Error (rRMSE2)`) %>%
      dplyr::slice(1)

    if (!(grepl("nca", rrmse2_best$`CL Method`, ignore.case = TRUE) &&
          grepl("nca", rrmse2_best$`Vd Method`, ignore.case = TRUE))) {

      nca_candidates <- base.out %>%
        dplyr::filter(
          grepl("nca", `CL Method`, ignore.case = TRUE) &
            grepl("nca", `Vd Method`, ignore.case = TRUE)
        )

      nca_better <- nca_candidates %>%
        dplyr::filter(
          `Relative Root Mean Squared Error (rRMSE1)` < rrmse2_best$`Relative Root Mean Squared Error (rRMSE1)`
        )

      if (nrow(nca_better) > 0) {
        base.best <- nca_better %>%
          dplyr::arrange(`Relative Root Mean Squared Error (rRMSE1)`) %>%
          dplyr::slice(1)
      }
    }
  }
  }

  base.ka.best <- base.best$`Calculated Ka`
  base.cl.best <- base.best$`Calculated CL`
  base.vd.best <- base.best$`Calculated Vd`

  message_text <- paste0(
    "Base PK parameter analysis finished. Estimated ka: ", base.ka.best,
    ", estimated CL: ", base.cl.best,
    ", estimated Vd: ", base.vd.best
  )
  if (verbose){
  cat(message_text, "\n")
################# 8. Parameter Sweeping on Vmax and Km #######################
  message(crayon::black(
    paste0("Run parameter sweeping on nonlinear elimination kinetics PK parameters",strrep(".", 20))))
  }

  sim.vmax.km.results.all <- NULL
  approx.vc.value <- sp_result$approx.vc.out$approx.vc.value

  # defensive fallback
  if (is.na(approx.vc.value)) {
    approx.vc.value <-suppressWarnings(c(signif(base.out[base.out$metrics.rank == 1, ]$`Calculated Vd` / 5, 3),
        signif(base.out[base.out$metrics.rank == 1, ]$`Calculated Vd` / 10, 3)))

  }

  if (route %in% c("bolus", "infusion")) {
    sim.vmax.km.results.all <- sim_sens_1cmpt_mm(
      dat = dat,
      # sim_vmax = list(mode = "auto", est.cl = base.cl.best),
      sim_vmax = list(mode = "auto", est.cl = base.out[base.out$metrics.rank<=2,]$`Calculated CL`),
      sim_km   = list(mode = "auto"),
      sim_vd   = list(mode = "manual", values = base.vd.best),
      sim_ka   = list(mode = "manual", values = NA),
      route    = "iv",
      verbose = verbose
    )
  }

  if (route %in% c("oral")) {
    sim.vmax.km.results.all <- sim_sens_1cmpt_mm(
      dat = dat,
      # sim_vmax = list(mode = "auto", est.cl = base.cl.best),
      sim_vmax = list(mode = "auto", est.cl = base.out[base.out$metrics.rank<=2,]$`Calculated CL`),
      sim_km   = list(mode = "auto"),
      sim_vd   = list(mode = "manual", values = c(approx.vc.value, base.vd.best)),
      sim_ka   = list(mode = "manual", values = base.ka.best),
      route    = "oral",
      verbose = verbose
    )
  }

  ########### 9. Parameter Sweeping on Multi-Compartmental Model Parameters#####
  if (verbose){
  message(crayon::black(
    paste0("Run parameter sweeping on multi-compartmental PK parameters",strrep(".", 20))))
  }
  # Collect identified vc from single-point extra and base.best.vd
  # Two-compartment model simulation
  sim.2cmpt.results.all <- NULL

  if (route %in% c("bolus", "infusion")) {
    sim.2cmpt.results.all <- sim_sens_2cmpt(
      dat = dat,
      sim_cl  = list(mode = "manual", values = base.cl.best),
      sim_vc  = list(mode = "manual", values = c( approx.vc.value, base.out[base.out$metrics.rank<=2,]$`Calculated Vd`)),
      sim_vp  = list(mode = "auto"),
      sim_q   = list(mode = "auto", auto.strategy = "scaled"),
      sim_ka  = list(mode = "manual", values = NA),
      route   = "iv",
      verbose = verbose
    )
  }

  if (route %in% c("oral")) {
    sim.2cmpt.results.all <- sim_sens_2cmpt(
      dat = dat,
      sim_cl  = list(mode = "manual", values = base.cl.best),
      sim_vc  = list(mode = "manual", values = c(approx.vc.value, base.out[base.out$metrics.rank<=2,]$`Calculated Vd`)),
      sim_vp  = list(mode = "auto"),
      sim_q   = list(mode = "auto", auto.strategy = "scaled"),
      sim_ka  = list(mode = "manual", values = base.ka.best),
      route   = "oral",
      verbose = verbose
     )
    }

  # Three-compartment model simulation
  sim.3cmpt.results.all <- NULL

  if (route %in% c("bolus", "infusion")) {
    sim.3cmpt.results.all <- sim_sens_3cmpt(
      dat = dat,
      sim_cl   = list(mode = "manual", values = base.cl.best),
      sim_vc   = list(mode = "manual", values = c(approx.vc.value, base.out[base.out$metrics.rank<=2,]$`Calculated Vd`)),
      sim_vp   = list(mode = "auto"),
      sim_vp2  = list(mode = "auto"),
      sim_q    = list(mode = "auto", auto.strategy = "scaled"),
      sim_q2   = list(mode = "auto", auto.strategy = "scaled"),
      sim_ka   = list(mode = "manual", values = NA),
      route    = "iv",
      verbose = verbose
    )
  }

  if (route %in% c("oral")) {
    sim.3cmpt.results.all <- sim_sens_3cmpt(
      dat = dat,
      sim_cl   = list(mode = "manual", values = base.cl.best),
      sim_vc   = list(mode = "manual", values = c(approx.vc.value, base.out[base.out$metrics.rank<=2,]$`Calculated Vd`)),
      sim_vp   = list(mode = "auto"),
      sim_vp2  = list(mode = "auto"),
      sim_q    = list(mode = "auto", auto.strategy = "scaled"),
      sim_q2   = list(mode = "auto", auto.strategy = "scaled"),
      sim_ka   = list(mode = "manual", values = base.ka.best),
      route    = "oral",
      verbose = verbose
    )
  }


  colnames(sim.vmax.km.results.all) <-
    c(
      "Simulated Vmax",
      "Simulated Km",
      "Simulated Vd",
      "Simulated Ka",
      "Absolute Predicted Error (APE)",
      "Mean Absolute Error (MAE)",
      "Mean Absolute Percentage Error (MAPE)",
      "Root Mean Squared Error (RMSE)",
      "Relative Root Mean Squared Error (rRMSE1)",
      "Relative Root Mean Squared Error (rRMSE2)",
      "Time spent"
    )
  colnames(sim.2cmpt.results.all) <-
    c(
      "Simulated Vc",
      "Simulated Vp",
      "Simulated Q",
      "Simulated CL",
      "Simulated Ka",
      "Absolute Predicted Error (APE)",
      "Mean Absolute Error (MAE)",
      "Mean Absolute Percentage Error (MAPE)",
      "Root Mean Squared Error (RMSE)",
      "Relative Root Mean Squared Error (rRMSE1)",
      "Relative Root Mean Squared Error (rRMSE2)",
      "Time spent"
    )
  colnames(sim.3cmpt.results.all) <-
    c(
      "Simulated Vc",
      "Simulated Vp",
      "Simulated Vp2",
      "Simulated Q",
      "Simulated Q2",
      "Simulated CL",
      "Simulated Ka",
      "Absolute Predicted Error (APE)",
      "Mean Absolute Error (MAE)",
      "Mean Absolute Percentage Error (MAPE)",
      "Root Mean Squared Error (RMSE)",
      "Relative Root Mean Squared Error (rRMSE1)",
      "Relative Root Mean Squared Error (rRMSE2)",
      "Time spent"
    )
############### 10. Residual error sigma estimation#####################

  method_additive <- .fbctrl$sigma_method_additive
  method_proportional <-.fbctrl$sigma_method_proportional
  sigma_fallback_fraction <- .fbctrl$sigma_fallback_fraction

  # Always compute full model-based result (reused if needed)
  sigma.out <- try(getsigma(dat),silent = T)
  sigma_add <- NA_real_
  sigma_prop <- NA_real_

  # Initialize flags for fallback detection
  used_add_fallback <- FALSE
  used_prop_fallback <- FALSE

  # --- Additive error model: decide method used ---
  if (method_additive == "model") {
    # Try model-based estimate
    sigma_add <- sigma.out$summary$sigma_additive
    if (is.na(sigma_add)) {
      used_add_fallback <- TRUE
      dv_obs <- dat[dat$EVID == 0, "DV", drop = TRUE]
      dv_mean <- mean(dv_obs, na.rm = TRUE, trim = 0.05)
      sigma_add <-  sigma_fallback_fraction *  dv_mean
    }
  } else if (method_additive == "fixed_fraction") {
    used_add_fallback <- TRUE
    dv_obs <- dat[dat$EVID == 0, "DV", drop = TRUE]
    dv_mean <- mean(dv_obs, na.rm = TRUE, trim = trim)
    sigma_add <-  sigma_fallback_fraction * dv_mean
  }

  # --- Proportional error model: decide method used ---
  if (method_proportional == "model") {
    sigma_prop <- sigma.out$summary$sigma_proportional
    if (is.na(sigma_prop)) {
      used_prop_fallback <- TRUE
      sigma_prop <- sigma_fallback_fraction
    }
  } else if (method_proportional == "fixed_fraction") {
    used_prop_fallback <- TRUE
    sigma_prop <-sigma_fallback_fraction
  }

  # --- Determine method labels used ---
  method_add_used <- if (used_add_fallback) {
    "Fallback (fixed fraction)"
  } else {
    "Model-based"
  }

  method_prop_used <- if (used_prop_fallback) {
    "Fallback (fixed fraction)"
  } else {
    "Model-based"
  }

  recommended_sigma_add_init <- sigma_add
  recommended_sigma_prop_init <-sigma_prop

################## 11. Parameter Selection Selection############################

  # Extract the column-wise minimum values for each performance metric
  mins <- sapply(sim.vmax.km.results.all[stat_cols], min, na.rm = TRUE)
  lmat <- sweep(sim.vmax.km.results.all[stat_cols], 2, mins, FUN = "==")
  sim.vmax.km.results.all$min_count <- rowSums(lmat, na.rm = TRUE)

  mins <- sapply(sim.2cmpt.results.all[stat_cols], min, na.rm = TRUE)
  lmat <- sweep(sim.2cmpt.results.all[stat_cols], 2, mins, FUN = "==")
  sim.2cmpt.results.all$min_count <- rowSums(lmat, na.rm = TRUE)

  mins <- sapply(sim.3cmpt.results.all[stat_cols], min, na.rm = TRUE)
  lmat <- sweep(sim.3cmpt.results.all[stat_cols], 2, mins, FUN = "==")
  sim.3cmpt.results.all$min_count <- rowSums(lmat, na.rm = TRUE)

  recommended_mm <- sim.vmax.km.results.all %>%
    dplyr::arrange(
      `Relative Root Mean Squared Error (rRMSE2)`,
      `Mean Absolute Percentage Error (MAPE)`
    ) %>%
    dplyr::filter(min_count == max(min_count)) %>%
    dplyr::slice(1)  # Take the first row in case of ties

  recommended.multi1 <- sim.2cmpt.results.all %>%
    dplyr::arrange(
      `Relative Root Mean Squared Error (rRMSE2)`,
      `Mean Absolute Percentage Error (MAPE)`
    ) %>%
    dplyr::filter(min_count == max(min_count)) %>%
    dplyr::slice(1)

  recommended.multi2 <- sim.3cmpt.results.all %>%
    dplyr::arrange(
      `Relative Root Mean Squared Error (rRMSE2)`,
      `Mean Absolute Percentage Error (MAPE)`
    ) %>%
    dplyr::filter(min_count == max(min_count)) %>%
    dplyr::slice(1)

  recommended_vmax_init <- recommended_mm$`Simulated Vmax`
  recommended_km_init <- recommended_mm$`Simulated Km`

  recommended_vc2cmpt_init <- recommended.multi1$`Simulated Vc`
  recommended_vp2cmpt_init <- recommended.multi1$`Simulated Vp`
  recommended_q2cmpt_init <- recommended.multi1$`Simulated Q`

  recommended_vc3cmpt_init <- recommended.multi2$`Simulated Vc`
  recommended_vp3cmpt_init <- recommended.multi2$`Simulated Vp`
  recommended_vp23cmpt_init <- recommended.multi2$`Simulated Vp2`
  recommended_q3cmpt_init <- recommended.multi2$`Simulated Q`
  recommended_q23cmpt_init <- recommended.multi2$`Simulated Q2`

  ############## 12. Finally selection####################
  #
  f_init_ka <- base.best$`Calculated Ka`[1]
  f_init_cl <- base.best$`Calculated CL`[1]
  f_init_vd <- base.best$`Calculated Vd`[1]
  f_init_vmax <- recommended_vmax_init
  f_init_km <- recommended_km_init
  # Multi-compartmental parameters
  f_init_vc2cmpt <-  recommended_vc2cmpt_init
  f_init_vp2cmpt <-  recommended_vp2cmpt_init
  f_init_q2cmpt <-   recommended_q2cmpt_init
  f_init_vc3cmpt <- recommended_vc3cmpt_init
  f_init_vp3cmpt <- recommended_vp3cmpt_init
  f_init_vp23cmpt <- recommended_vp23cmpt_init
  f_init_q3cmpt <- recommended_q3cmpt_init
  f_init_q23cmpt <- recommended_q23cmpt_init

  # Method identification
  if (!hybrid.base) {
    # Unified method — all parameters from one method
    method_name <- base.best$Method[1]

    # Use same method for all parameters
    sel.method.ka <- method_name
    sel.method.cl <- method_name
    sel.method.vd <- method_name
    sel.method.multi <- "Parameter sweeping"

    # Map method name to corresponding Ka estimation label
    ka_label_map <- c(
      "Adaptive single-point method" = "Adaptive single-point method",
      "Graphic methods" = "Graphic methods",
      "Naive pooled NCA (FD)" = "Wanger-nelson method",
      "Naive pooled NCA (MD)" = "Wanger-nelson method",
      "Naive pooled NCA (all)" = "Wanger-nelson method"
    )

    # Apply mapping
    sel.method.ka <- if (!is.na(ka_method)) ka_label_map [[ka_method]] else "IV"

    # Fallback handling for Ka (if applicable)
    if (method_name == "Adaptive single-point method" &&
        route == "oral" && isTRUE(used_sp_ka_fallback)) {
      sel.method.ka <- "Fallback (fixed Ka)"
    }

    if (method_name == "Graphic methods" &&
        isTRUE(used_graph_ka_fallback)) {
      sel.method.ka <- "Fallback (fixed Ka)"
    }

    if (startsWith(method_name, "Naive pooled NCA") &&
        isTRUE(used_nca_ka_fallback)) {
      sel.method.ka <- "Fallback (fixed Ka)"
    }
  }

  if (hybrid.base) {
    # Hybrid mode — parameter-level method assignment
    sel.method.multi <- "Parameter sweeping"

    # Extract method source names for each parameter
    ka_method <- base.best$`Ka Method`[1]
    cl_method <- base.best$`CL Method`[1]
    vd_method <- base.best$`Vd Method`[1]

    # Mapping of method source
    method_label_map <- c(
      simpcal = "Adaptive single-point method",
      graph   = "Graphic methods",
      nca_fd  = "Naive pooled NCA (FD)",
      nca_efd = "Naive pooled NCA (MD)",
      nca_all = "Naive pooled NCA (all)"
    )
    # Assign base method labels
    sel.method.ka <- if (!is.na(ka_method)) method_label_map[[ka_method]] else "IV"
    sel.method.cl <- method_label_map[[cl_method]]
    sel.method.vd <- method_label_map[[vd_method]]

    # Fallback handling for Ka (if needed)
    if (ka_method == "simpcal" &&
        route == "oral" && isTRUE(used_sp_ka_fallback)) {
      sel.method.ka <- "Fallback (fixed Ka)"
    }

    if (ka_method == "graph" && isTRUE(used_graph_ka_fallback)) {
      sel.method.ka <- "Fallback (fixed Ka)"
    }

    if (ka_method %in% c("nca_fd", "nca_efd", "nca_all") &&
        isTRUE(used_nca_ka_fallback)) {
      sel.method.ka <- "Fallback (fixed Ka)"
    }
  }

  init.params.out.ka <- data.frame(method = sel.method.ka,
                                   vd = f_init_ka)
  init.params.out.cl <- data.frame(method = sel.method.cl,
                                   vd = f_init_cl)
  init.params.out.vd <- data.frame(method = sel.method.vd,
                                   vd = f_init_vd)
  init.params.out.vmax.km <-
    data.frame(method = sel.method.multi,
               vmax = f_init_vmax,
               km =  f_init_km)

  init.params.out.vc.vp <- data.frame(
    method  =  sel.method.multi,
    vc2cmpt =  f_init_vc2cmpt,
    vp2cmpt =   f_init_vp2cmpt,
    q2cmpt =   f_init_q2cmpt,
    vc3cmpt =  f_init_vc3cmpt,
    vp3cmpt =  f_init_vp3cmpt,
    vp23cmpt =   f_init_vp23cmpt,
    q3cmpt =   f_init_q3cmpt,
    q23cmpt =   f_init_q23cmpt
  )

  colnames(init.params.out.ka) <- c("Method", "Ka")
  colnames(init.params.out.cl) <- c("Method", "CL")
  colnames(init.params.out.vd) <- c("Method", "Vd")
  colnames(init.params.out.vmax.km) <- c("Method", "Vmax", "Km")
  colnames(init.params.out.vc.vp) <-
    c(
      "Method",
      "Vc2cmpt",
      "Vp2cmpt",
      "Q2cmpt",
      "Vc3cmpt",
      "Vp3cmpt",
      "Vp23cmpt",
      "Q3cmpt",
      "Q23cmpt"
    )

  # Create sigma output data frame
  init.params.out.sigma <- data.frame(
    Method_add = method_add_used,
    Method_prop = method_prop_used,
    Sigma_additive = sigma_add,
    Sigma_proportional = sigma_prop
  )

  colnames(init.params.out.sigma) <-
    c("Method (sigma additive)", "Method (sigma proportional)" , "Sigma additive", "Sigma proportional")

  init.params.out.all <- list(
    init.params.ka = init.params.out.ka,
    init.params.cl = init.params.out.cl,
    init.params.vd = init.params.out.vd,
    init.params.vmax.km = init.params.out.vmax.km,
    init.params.multi = init.params.out.vc.vp,
    init.params.sigma = init.params.out.sigma
  )

  Recommended_inits_df = data.frame(
    Parameters = c(
      "Ka",
      "CL",
      "Vd",
      "Vmax",
      "Km",
      "Vc(2CMPT)",
      "Vp(2CMPT)",
      "Q(2CMPT)",
      "Vc(3CMPT)",
      "Vp(3CMPT)",
      "Vp2(3CMPT)",
      "Q(3CMPT)",
      "Q2(3CMPT)",
      "Sigma additive",
      "Sigma proportional"
    ),

    Methods = c(
      init.params.out.all$init.params.ka$Method,
      init.params.out.all$init.params.cl$Method,
      init.params.out.all$init.params.vd$Method,
      init.params.out.all$init.params.vmax.km$Method,
      init.params.out.all$init.params.vmax.km$Method,
      init.params.out.all$init.params.multi$Method,
      init.params.out.all$init.params.multi$Method,
      init.params.out.all$init.params.multi$Method,
      init.params.out.all$init.params.multi$Method,
      init.params.out.all$init.params.multi$Method,
      init.params.out.all$init.params.multi$Method,
      init.params.out.all$init.params.multi$Method,
      init.params.out.all$init.params.multi$Method,
      init.params.out.all$init.params.sigma$`Method (sigma additive)`,
      init.params.out.all$init.params.sigma$`Method (sigma proportional)`

    ),
    Values = c(
      init.params.out.all$init.params.ka$Ka,
      init.params.out.all$init.params.cl$CL,
      init.params.out.all$init.params.vd$Vd,
      init.params.out.all$init.params.vmax.km$Vmax,
      init.params.out.all$init.params.vmax.km$Km,
      init.params.out.all$init.params.multi$Vc2cmpt,
      init.params.out.all$init.params.multi$Vp2cmpt,
      init.params.out.all$init.params.multi$Q2cmpt,
      init.params.out.all$init.params.multi$Vc3cmpt,
      init.params.out.all$init.params.multi$Vp3cmpt,
      init.params.out.all$init.params.multi$Vp23cmpt,
      init.params.out.all$init.params.multi$Q3cmpt,
      init.params.out.all$init.params.multi$Q23cmpt,
      init.params.out.all$init.params.sigma$`Sigma additive`,
      init.params.out.all$init.params.sigma$`Sigma proportional`
    )
  )

  Recommended_inits_df$Values <- format(Recommended_inits_df$Values,
                                        scientific = FALSE,
                                        digits = 3)

  init.history <- list(
    base.out = base.out,
    single.point.out = sp_result,
    nca.out = nca_out,
    ka.wanger.nelson.out =  ka_wn,
    graph.out =   graph_out,
    sim.vmax.km = sim.vmax.km.results.all,
    sim.2cmpt = sim.2cmpt.results.all,
    sim.3cmpt = sim.3cmpt.results.all,
    sigma.out = sigma.out
  )

  params.descriptions <- c(
    "Ka: absorption constant rate",
    "CL: clearance",
    "Vd: volume of distribution",
    "Vmax: maximum metabolic rate",
    "Km: Michaelis constant",
    "Vc: volume of distribution of the central compartment",
    "Vp: volume of distribution of the peripheral compartment",
    "Vp2: volume of distribution of the second peripheral compartment",
    "Q: inter-compartmental clearance",
    "Q2: inter-compartmental clearance between central and second peripheral compartment",
    "Sigma additive: standard deviation of additive residual error",
    "Sigma proportional: standard deviation of proportional residual error"
  )

  end.time <- Sys.time()
  time.spent <-
    round(as.numeric(difftime(end.time, start.time, units = "secs")), 3)

  output_env <- new.env()
  output_env$Datainfo <- Datainfo
  output_env$Recommended_initial_estimates <- Recommended_inits_df
  output_env$Run.history <- init.history
  output_env$time.spent <- time.spent
  output_env$Parameter.descriptions <- params.descriptions
  output_env$Omegas<- getOmegas()

  class(output_env) <- "getPPKinits"

  return(output_env)

} # end of function


#' Print method for `getPPKinits` objects
#'
#' Prints a summary of the results from the initial parameter estimation pipeline,
#' including recommended initial estimates, ETA variance estimates, and parameter descriptions.
#' It is the default S3 `print` method for objects of class `getPPKinits`.
#'
#' @param x An object of class `getPPKinits` containing the initial parameter estimation results.
#' Expected components include:
#' \itemize{
#'   \item \code{Recommended_initial_estimates}: A data frame with estimated values and selection methods.
#'   \item \code{Parameter.descriptions}: A character vector explaining the meaning of each parameter.
#'   \item \code{time.spent}: Time taken to compute the estimates.
#' }
#' @param ... Additional arguments (for compatibility with the generic \code{print()}).
#'
#' @return Prints a formatted summary to the console.
#'
#' @examples
#' \donttest{
#' ## Oral example
#' inits.out <- getPPKinits(Bolus_1CPT)
#' print(inits.out)
#' }
#'
#' @export

print.getPPKinits <- function(x, ...) {

  cat("===============Initial Parameter Estimation Summary ===============\n")
  cat("\nRecommended initial estimates :\n")
  print(utils::head(x$Recommended_initial_estimates, 15))

  cat("\nTime spent :\n")
  print(paste0(x$time.spent,"s"))

  if (!is.null(x$Omegas)) {
    cat("\nETA variances and derived covariances:\n")
    print(x$Omegas)
    cat("Note: The ETA variances and covariances listed above are predefined default initialization values automatically assigned by the package.\n")
  }

  cat("\nParameter descriptions:\n")
  print(x$Parameter.descriptions)

  cat("\n=============== End of Summary ===============\n")
}