R/get_bootSummaryNlme.R

In Certara.Xpose.NLME: Enhances 'xpose' Diagnostics for Pharmacometric Models from 'Certara.RsNLME' and Phoenix NLME

#' @title Build a fused fit + bootstrap parameter summary
#'
#' @description Companion to [get_summaryNlme()] for
#'   `Certara.RsNLME::bootstrap()` results. Combines original-fit columns
#'   (`Estimate`, `%RSE`, `Shrinkage (%)`) with per-replicate bootstrap
#'   columns (`Bootstrap estimate (<metric>)`, `Bootstrap <pct>% CI`)
#'   using the same transform / shrinkage machinery.
#'
#' @details
#' Three input modes drive what columns appear in the output:
#'
#' \describe{
#'   \item{m1: `bootResult` only}{No original-fit source; output is
#'     bootstrap-only -- `Section`, `Parameter`, optional `Unit`,
#'     `Bootstrap estimate (<metric>)`, and `Bootstrap <pct>% CI`. A
#'     `message()` notes that original-fit columns are unavailable and
#'     how to include them (`initialEstimates = TRUE` or `xpdb`). The
#'     residual-transform advisory (see below) still fires, since the
#'     default's fitness can't be ruled out without PML either way.}
#'   \item{m2: `bootResult` carries an embedded `fitSummary`}{`fitSummary`
#'     is the original-fit source (when `initialEstimates = TRUE` was used
#'     on the bootstrap call). Omega off-diagonal covariance rows
#'     (`Diagonal = FALSE`) are dropped so the original-fit side stays
#'     variance-scale, matching m3's "off-diagonals excluded" contract;
#'     older fitSummary tables without a `Diagonal` column are treated as
#'     diagonal-only.}
#'   \item{m3: explicit `xpdb` argument}{Original-fit columns come from
#'     `get_summaryNlme(xpdb, ...)`. If `bootResult$fitSummary` is also
#'     non-NULL a one-shot warning fires and `xpdb` wins.}
#' }
#'
#' Per-replicate filtering: replicates whose `BootOverall$ReturnCode`
#' falls outside `return_code_ok` are dropped silently before the metric
#' and percentile CI are computed. The `BootOverall` table on the
#' `rsnlme_boot` is the single source of truth -- no on-disk reads.
#'
#' The CI bounds are the empirical `(1 - ci_level) / 2` and
#' `1 - (1 - ci_level) / 2` quantiles of the kept replicates, computed
#' with `stats::quantile()`'s default method (type 7).
#'
#' Soft-degrade for missing stacks: `BootSigmaStacked` and
#' `BootSecondaryStacked` are net-new in the corresponding
#' `Certara.NLME8` release. When either is `NULL` (older NLME8 build) the
#' function emits NA bootstrap cells for that section and a single warning
#' naming the missing stacks plus the installed `Certara.NLME8` version.
#' `BootThetaStacked` and `BootOmegaStacked` predate the new release and
#' are always available.
#'
#' Output formatting (shared with [get_summaryNlme()]): when a
#' non-identity transform is active the values are shown on the
#' transformed scale only, and the scale label is appended to the shared
#' `Parameter` name in parentheses -- `nV (CV%)`, `CEps (SD)`, or a custom
#' `name` for `list(fn=, dfn=, name=)` specs. Identity / `raw` rows keep
#' the bare name. The `Unit` column carries real units only and is
#' dropped when every row is dimensionless. The
#' bootstrap CI is a single character column formatted as `"lo - hi"` (a
#' dash range, no brackets). Numeric `digits` is applied via `signif()` to
#' both the bootstrap columns and, when present, the original-fit
#' `Estimate`, `%RSE`, and `Shrinkage (%)`.
#'
#' Residual-shape advisories: a per-sigma warning fires when PML
#' classifies a sigma as additive or otherwise non-proportional while it
#' is reported with the `multiplicative_cv` default. Classification
#' requires PML, which is only available when `xpdb` is supplied (m3); on
#' m1/m2 no per-sigma warning is emitted. The general default-transform
#' `message()` is broader: it fires for any defaulted sigma that isn't
#' *proven* proportional, which on m1/m2 (no PML to check) means it
#' always fires -- mirroring [get_summaryNlme()]'s own no-PML behaviour.
#' Set `options(xposeNlme.summary.quiet_default_warning = TRUE)` to
#' silence it.
#'
#' @param bootResult Object returned by `Certara.RsNLME::bootstrap()`
#'   (an `rsnlme_boot` list of CSV-derived tables). Must carry at least
#'   `BootOverall` (with `Replicate` + `ReturnCode` columns),
#'   `BootThetaStacked`, and `BootOmegaStacked`.
#' @param xpdb Optional `xpose_data` object from `xposeNlme()` /
#'   `xposeNlmeModel()`. When supplied, used as the original-fit source.
#' @param transform Named list of per-parameter transforms keyed by
#'   parameter label. Same contract as `get_summaryNlme()`'s `transform`
#'   argument: either a preset string from the section's catalog or a
#'   `list(fn = ..., dfn = ..., name = ...)` spec, where the optional
#'   `name` sets the scale flag appended to the parameter name (a warning
#'   fires for a custom transform supplied without `name`). Applied to
#'   both the original-fit rows and the per-replicate bootstrap pool.
#' @param units Named character vector (or list of length-1 character)
#'   keyed by parameter label, overriding the real-units `Unit` column for
#'   any row. Forwarded to `get_summaryNlme()` (m3) / the fitSummary
#'   pipeline (m2) and applied to the bootstrap-only output on m1. Keys
#'   not present in the relevant section emit a single warning and are
#'   ignored.
#' @param metric `"Median"` (default) or `"Mean"` -- statistic for the
#'   `Bootstrap estimate (<metric>)` column.
#' @param ci_level Numeric in (0, 1); width of the percentile CI. When
#'   `NULL` (default) it inherits the bootstrap run's `confidenceLevel`
#'   (stored on `bootResult`), falling back to 0.95 when that metadata is
#'   absent. An explicit value always overrides.
#' @param return_code_ok Integer vector of `ReturnCode` values to accept;
#'   replicates outside this set are dropped. Defaults to `1:3`.
#' @param digits Significant-digits count applied via `signif()` to the
#'   bootstrap columns and, when an original-fit summary is present, to the
#'   stored numeric `Estimate`, `%RSE`, and `Shrinkage (%)`. Defaults to 3.
#'   The `print.bootSummaryNlme` method also honours `digits` for displayed
#'   precision (by setting `pillar.sigfig` for the duration of the print),
#'   so the original-fit columns show the same number of significant
#'   figures on screen as are stored -- matching the pre-formatted
#'   bootstrap columns, which already carry `digits` in their character
#'   values.
#'
#' @return A tibble (class `bootSummaryNlme`) with `Section`, `Parameter`
#'   (carrying a `(CV%)` / `(SD)` / custom scale flag when a non-identity
#'   transform is active), optional original-fit columns (`Estimate`,
#'   `%RSE`, `Shrinkage (%)`), an optional real-units `Unit` column, and
#'   the two bootstrap columns (the CI formatted as `"lo - hi"`). Carries
#'   the attributes `n_used`, `n_total`, `ci_level`, `return_code_ok`,
#'   `transform`, `metric`, `fitSource` (`"xpdb"` / `"embedded"` /
#'   `"none"`). The `bootSummaryNlme` class carries a `print` method that
#'   honours the `digits` argument for displayed precision.
#'
#' @examples
#' \dontrun{
#' fit_boot <- Certara.RsNLME::bootstrap(model, ...)
#'
#' # m1: bootstrap-only summary, no fit source.
#' get_bootSummaryNlme(fit_boot)
#'
#' # m2: fused summary using the bootstrap's embedded fitSummary
#' # (initialEstimates = TRUE on the bootstrap call).
#' fit_boot_init <- Certara.RsNLME::bootstrap(model,
#'                                            initialEstimates = TRUE, ...)
#' get_bootSummaryNlme(fit_boot_init)
#'
#' # m3: fused summary against an explicit xpdb -- preferred when the
#' # xpdb has its own provenance (covariates, residuals, posthoc) that
#' # fitSummary doesn't capture.
#' xp <- xposeNlmeModel(fit_model)
#' get_bootSummaryNlme(fit_boot, xpdb = xp)
#'
#' # Custom transform on a sigma applied to both original-fit and
#' # bootstrap columns; `name` sets the scale flag on the parameter name.
#' get_bootSummaryNlme(
#'   fit_boot_init,
#'   transform = list(
#'     CEps = list(
#'       fn  = function(s) 200 * s,
#'       dfn = function(s) 200,
#'       name = "2xCV%"
#'     )
#'   )
#' )
#' }
#'
#' @seealso [get_summaryNlme()]
#' @importFrom magrittr %>%
#' @export
get_bootSummaryNlme <- function(bootResult,
                                xpdb = NULL,
                                transform = list(),
                                units = list(),
                                metric = c("Median", "Mean"),
                                ci_level = NULL,
                                return_code_ok = 1:3,
                                digits = 3) {
  metric <- match.arg(metric)
  if (!is.null(ci_level) &&
      (!is.numeric(ci_level) || length(ci_level) != 1L ||
       !is.finite(ci_level) || ci_level <= 0 || ci_level >= 1)) {
    stop("`ci_level` must be a single number in (0, 1).", call. = FALSE)
  }
  if (!is.list(transform)) {
    stop("`transform` must be a named list.", call. = FALSE)
  }
  if (length(transform) &&
      (is.null(names(transform)) || any(!nzchar(names(transform))))) {
    stop("Every entry of `transform` must be named.", call. = FALSE)
  }
  if (!is.list(units) && !is.character(units)) {
    stop("`units` must be a named list or a named character vector.",
         call. = FALSE)
  }
  if (!is.numeric(digits) || length(digits) != 1L || digits < 1) {
    stop("`digits` must be a single positive integer.", call. = FALSE)
  }

  ci_level <- .resolve_ci_level(bootResult, ci_level)
  if (!is.numeric(ci_level) || length(ci_level) != 1L ||
      !is.finite(ci_level) || ci_level <= 0 || ci_level >= 1) {
    stop("Resolved `ci_level` is not a single number in (0, 1).",
         call. = FALSE)
  }

  fit_source <- .resolve_fit_source(bootResult, xpdb)

  # When no original-fit (fitmodel) summary is available or used, the table is
  # restricted to bootstrap-only columns and residual-error classification is
  # impossible. Surface that clearly rather than degrading silently.
  if (fit_source == "none") {
    message(
      "No original-fit (fitmodel) summary is available for this bootstrap ",
      "result. Original-fit Estimate / %RSE / Shrinkage (%) are ",
      "unavailable. Re-run bootstrap() with initialEstimates = TRUE, or ",
      "pass `xpdb`, to include them."
    )
  }

  # The scale flag is appended once to the shared `Parameter` column
  # post-merge (see below), so the LHS paths build bare labels to keep the
  # join key aligned with the RHS. Advisories are emitted once by
  # `.emit_boot_sigma_advisories()` below, so suppress the LHS pipeline's.
  lhs <- switch(fit_source,
    xpdb     = get_summaryNlme(xpdb, transform = transform, units = units,
                               digits = digits, append_flag = FALSE,
                               emit_advisories = FALSE),
    embedded = .build_lhs_from_fitSummary(bootResult$fitSummary,
                                          transform, units, digits = digits),
    none     = NULL
  )

  # On m1 we apply user `transform` / `units` to the RHS rows directly
  # (no LHS pipeline runs to consume and validate them). Mirror
  # `.normalise_transform()` / `.normalise_units()` behaviour from
  # `get_summaryNlme.R`: warn-and-drop keys that aren't present in any
  # bootstrap section. For m2/m3 the LHS pipeline validates against its
  # own (richer) label set, and `.merge_lhs_rhs()` drops the RHS `Unit`
  # column when both sides supply one, so we deliberately don't forward
  # either argument to the RHS builders in those modes.
  rhs_transform <- if (fit_source == "none") {
    .normalise_transform_rhs(transform, bootResult)
  } else {
    transform
  }
  rhs_units <- if (fit_source == "none") {
    .normalise_units_rhs(units, bootResult)
  } else {
    list()
  }

  bo <- .require_bootOverall(bootResult)
  n_total <- nrow(bo)
  keep <- bo$Replicate[bo$ReturnCode %in% return_code_ok]
  n_used <- length(keep)

  rhs_pieces <- list()
  missing_stacks <- character()

  # NLME8's `BootThetaStacked` mirrors the engine's "fixedEffects" block,
  # which historically conflates true thetas with residual sigmas
  # (plan §1b: "thetas-and-sigmas"). Now that `BootSigmaStacked` is the
  # dedicated home for sigmas, drop sigma-named rows from the Fixed
  # effects RHS so each parameter appears in exactly one Section. Mirrors
  # the equivalent guard in `Certara.RsNLME::print.rsnlme_boot`. Falls
  # back to the unfiltered stack only when neither `BootSigmaCI` nor
  # `BootSigmaStacked` is available (older NLME8 build pre-§1b release).
  theta_stack <- .drop_sigma_rows_from_theta_stack(
    bootResult$BootThetaStacked,
    bootResult$BootSigmaCI,
    bootResult$BootSigmaStacked
  )
  rhs_pieces[["Fixed effects"]] <- .build_rhs_section(
    stack = theta_stack, label_col = "Theta",
    section = "Fixed effects",
    keep = keep, transform = rhs_transform, units = rhs_units,
    metric = metric, ci_level = ci_level, digits = digits
  )

  rhs_pieces[["Random effects"]] <- .build_rhs_omega_section(
    stack = bootResult$BootOmegaStacked, ci_table = bootResult$BootOmegaCI,
    keep = keep, transform = rhs_transform, units = rhs_units,
    metric = metric, ci_level = ci_level, digits = digits
  )

  if (is.null(bootResult$BootSigmaStacked)) {
    sigma_names <- .infer_param_names(bootResult$BootSigmaCI, "Sigma")
    if (length(sigma_names)) {
      missing_stacks <- c(missing_stacks, "BootSigmaStacked")
      rhs_pieces[["Residual error"]] <- .na_rhs_section(
        sigma_names, "Residual error", metric, ci_level
      )
    }
  } else {
    rhs_pieces[["Residual error"]] <- .build_rhs_section(
      stack = bootResult$BootSigmaStacked, label_col = "Sigma",
      section = "Residual error",
      keep = keep, transform = rhs_transform, units = rhs_units,
      metric = metric, ci_level = ci_level, digits = digits
    )
  }

  if (is.null(bootResult$BootSecondaryStacked)) {
    sec_names <- .infer_param_names(bootResult$BootSecondary, "Parameter")
    if (length(sec_names)) {
      missing_stacks <- c(missing_stacks, "BootSecondaryStacked")
      rhs_pieces[["Secondary"]] <- .na_rhs_section(
        sec_names, "Secondary", metric, ci_level
      )
    }
  } else {
    rhs_pieces[["Secondary"]] <- .build_rhs_section(
      stack = bootResult$BootSecondaryStacked, label_col = "Secondary",
      section = "Secondary",
      keep = keep, transform = rhs_transform, units = rhs_units,
      metric = metric, ci_level = ci_level, digits = digits
    )
  }

  rhs <- dplyr::bind_rows(rhs_pieces)

  if (length(missing_stacks)) {
    .warn_missing_stacks(missing_stacks)
  }

  result <- .merge_lhs_rhs(lhs, rhs)
  result <- .order_by_section(result)

  # Residual-shape advisories are driven by PML classification, which is only
  # available from `xpdb`. Sigma labels are read while `Parameter` is still
  # bare (before the flag is appended).
  sigma_labels <- result$Parameter[result$Section == "Residual error"]
  .emit_boot_sigma_advisories(sigma_labels, transform,
                              if (!is.null(xpdb)) xpdb$code else NULL)

  # Append the scale flag once to the shared `Parameter` column, after the
  # bare-label join. Applies uniformly to original-fit, bootstrap, and
  # soft-degrade rows.
  result <- .append_flags_to_table(result, transform)

  # Drop the `Unit` column only when *every* row is the empty string --
  # i.e. truly dimensionless. `.na_rhs_section()` deliberately uses
  # `NA_character_` to mean "unknown" (the soft-degrade sigma /
  # secondary path where the per-replicate stack is missing) and that
  # information must survive: an explicit NA cell signals "we don't
  # know" while "" signals "no unit applies". Without this the
  # `nzchar(NA) -> NA` rule would silently drop the column whenever a
  # purely soft-degraded result has no nzchar entries to keep it
  # alive.
  if ("Unit" %in% names(result) &&
      !any(is.na(result$Unit) | nzchar(result$Unit))) {
    result$Unit <- NULL
  }

  attr(result, "n_used") <- n_used
  attr(result, "n_total") <- n_total
  attr(result, "ci_level") <- ci_level
  attr(result, "return_code_ok") <- return_code_ok
  attr(result, "transform") <- transform
  attr(result, "units") <- units
  attr(result, "metric") <- metric
  attr(result, "fitSource") <- fit_source

  .as_bootSummaryNlme(result, digits)
}


# --------------------------------------------------------------------------
# Display class: let `digits` drive printed precision (mirrors
# `get_summaryNlme()`'s `.as_summaryNlme()` / `print.summaryNlme()`)
# --------------------------------------------------------------------------

# Tag the result tibble so `print.bootSummaryNlme()` can map `digits` onto
# `pillar.sigfig` at display time. The original-fit numerics are already
# `signif()`-rounded in `.apply_summary_pipeline()` when `digits` is
# non-NULL (`get_bootSummaryNlme()` defaults `digits` to 3, never NULL);
# this only governs how many significant figures the tibble print method
# shows. The bootstrap columns are pre-formatted character strings via
# `.fmt()` and are unaffected by `pillar.sigfig` either way.
.as_bootSummaryNlme <- function(x, digits) {
  attr(x, "summary_digits") <- digits
  class(x) <- c("bootSummaryNlme", class(x))
  x
}

#' @rdname get_bootSummaryNlme
#' @param x A `bootSummaryNlme` tibble returned by `get_bootSummaryNlme()`.
#' @param ... Further arguments passed to the underlying tibble print method.
#' @export
print.bootSummaryNlme <- function(x, ...) {
  digits <- attr(x, "summary_digits")
  sigfig <- if (is.null(digits)) 15L else max(1L, min(as.integer(digits), 15L))
  old <- options(pillar.sigfig = sigfig)
  on.exit(options(old), add = TRUE)
  NextMethod()
  invisible(x)
}

# Normalise the user `transform` argument for the m1 (RHS-only) path:
# warn-and-drop keys that aren't present in any bootstrap section. Shape
# validation (is.list + every entry named) already happened in
# `get_bootSummaryNlme()` so this helper only handles the label-set
# check. Mirrors `.normalise_transform()` in `get_summaryNlme.R`, but
# the label universe is different (RHS stacks instead of
# `prmTable$label`).
.normalise_transform_rhs <- function(transform, bootResult) {
  if (!length(transform)) return(list())
  labels <- .all_rhs_labels(bootResult)
  unknown <- setdiff(names(transform), labels)
  if (length(unknown)) {
    warning("Ignored `transform` keys not in any bootstrap section: ",
            paste(unknown, collapse = ", "), ".", call. = FALSE)
    transform <- transform[setdiff(names(transform), unknown)]
  }
  transform
}

# Normalise the user `units` argument for the m1 (RHS-only) path: coerce
# list -> named character, validate names against the union of labels
# present in the RHS sources (theta + omega + sigma + secondary), and
# warn-and-drop on unknown keys. Mirrors `.normalise_units()` in
# `get_summaryNlme.R`, but the label universe is different (RHS stacks
# instead of `prmTable$label`). We don't re-export `.normalise_units()`
# because the two callers need slightly different error messaging.
.normalise_units_rhs <- function(units, bootResult) {
  if (!length(units)) return(stats::setNames(character(), character()))
  if (is.list(units)) units <- unlist(units)
  if (!is.character(units) || is.null(names(units)) ||
      any(!nzchar(names(units)))) {
    stop("`units` must be a named character vector or list.",
         call. = FALSE)
  }
  labels <- .all_rhs_labels(bootResult)
  unknown <- setdiff(names(units), labels)
  if (length(unknown)) {
    warning("Ignored `units` keys not in any bootstrap section: ",
            paste(unknown, collapse = ", "), ".", call. = FALSE)
    units <- units[setdiff(names(units), unknown)]
  }
  units
}

# Union of parameter labels visible in the RHS output. Sigma rows are
# read from `BootSigmaCI` first (Piece 1's always-emit-headers rule
# guarantees presence even when header-only), falling back to the
# stacked table; same for secondaries. Diagonal-only filtering for
# omegas mirrors `.build_rhs_omega_section()` so the union doesn't
# include off-diagonals that the output never shows.
.all_rhs_labels <- function(bootResult) {
  out <- character()
  ts <- bootResult$BootThetaStacked
  if (!is.null(ts) && "Theta" %in% names(ts)) {
    out <- c(out, unique(as.character(ts$Theta)))
  }
  oc <- bootResult$BootOmegaCI
  if (!is.null(oc) && all(c("Omega", "Diagonal") %in% names(oc))) {
    out <- c(out, unique(as.character(oc$Omega[.isTRUE_logical(oc$Diagonal)])))
  } else {
    os <- bootResult$BootOmegaStacked
    if (!is.null(os) && "Omega" %in% names(os)) {
      diag_only <- setdiff(unique(os$Omega),
                           grep("_", unique(os$Omega), value = TRUE))
      out <- c(out, diag_only)
    }
  }
  for (entry in list(
    list(tbl = bootResult$BootSigmaCI,        col = "Sigma"),
    list(tbl = bootResult$BootSigmaStacked,   col = "Sigma"),
    list(tbl = bootResult$BootSecondary,      col = "Parameter"),
    list(tbl = bootResult$BootSecondaryStacked, col = "Secondary")
  )) {
    tbl <- entry$tbl
    if (!is.null(tbl) && entry$col %in% names(tbl)) {
      out <- c(out, unique(as.character(tbl[[entry$col]])))
    }
  }
  unique(out)
}


# --------------------------------------------------------------------------
# Fit-source resolution
# --------------------------------------------------------------------------

.resolve_fit_source <- function(bootResult, xpdb) {
  if (!is.null(xpdb)) {
    if (!is.null(bootResult$fitSummary)) {
      warning(
        "Both `xpdb` and `bootResult$fitSummary` are present; ",
        "`xpdb` wins. Pass `xpdb = NULL` to use the embedded fitSummary instead.",
        call. = FALSE
      )
    }
    return("xpdb")
  }
  if (!is.null(bootResult$fitSummary) && nrow(bootResult$fitSummary) > 0L) {
    return("embedded")
  }
  "none"
}


# --------------------------------------------------------------------------
# LHS via fitSummary -> prm_like -> .apply_summary_pipeline
# --------------------------------------------------------------------------

# RsNLME's fitSummary carries `Parameter | Type | Estimate | SE | %RSE |
# Shrinkage`. Reshape to the prm_like shape that `.apply_summary_pipeline()`
# accepts, then route through the same pipeline `get_summaryNlme()` uses.
# Shrinkage values arrive on the percent scale (RsNLME's
# `.parseEtaShrinkages` / `.parseEpsShrinkages` multiply at parse time
# per Plan §2b), so they pass through to the pipeline unchanged.
.build_lhs_from_fitSummary <- function(fitSummary, transform, units = list(),
                                       digits = NULL) {
  if (is.null(fitSummary) || !nrow(fitSummary)) {
    return(NULL)
  }

  required_cols <- c("Parameter", "Type", "Estimate", "SE", "Shrinkage")
  missing_cols <- setdiff(required_cols, names(fitSummary))
  if (length(missing_cols)) {
    stop(
      "`bootResult$fitSummary` is missing required columns: ",
      paste(missing_cols, collapse = ", "),
      ". Did you build it via Certara.RsNLME (>= the bootstrap-summary release)?",
      call. = FALSE
    )
  }

  # `Diagonal` distinguishes omega variances from block off-diagonal
  # covariances. We drop the covariance rows via `.drop_offDiagonal()`
  # (the same filter the m3/xpdb path applies to `prmTable`) so the
  # transformed LHS stays variance-scale, matching `get_summaryNlme()`'s
  # "off-diagonals excluded" contract. Older RsNLME builds omit the
  # column entirely (diagonal-only fitSummary), so default to TRUE for
  # backward compatibility.
  diagonal <- if ("Diagonal" %in% names(fitSummary)) {
    as.logical(fitSummary$Diagonal)
  } else {
    TRUE
  }
  prm_like <- tibble::tibble(
    type     = fitSummary$Type,
    label    = fitSummary$Parameter,
    value    = fitSummary$Estimate,
    se       = fitSummary$SE,
    diagonal = diagonal,
    section  = unname(.section_for[fitSummary$Type])
  )
  prm_like <- .drop_offDiagonal(prm_like)

  shrink_map <- stats::setNames(fitSummary$Shrinkage, fitSummary$Parameter)
  shrink_map <- shrink_map[!is.na(shrink_map)]

  .apply_summary_pipeline(
    prm_like   = prm_like,
    transform  = transform,
    units      = units,
    shrink_map = shrink_map,
    param_units = NULL,
    sigma_roles = NULL,
    append_flag = FALSE,
    digits = digits,
    emit_advisories = FALSE
  )
}


# --------------------------------------------------------------------------
# BootOverall guard
# --------------------------------------------------------------------------

.require_bootOverall <- function(bootResult) {
  bo <- bootResult$BootOverall
  if (is.null(bo) || !is.data.frame(bo) ||
      !all(c("Replicate", "ReturnCode") %in% names(bo))) {
    stop(
      "`bootResult$BootOverall` must be a data frame with `Replicate` ",
      "and `ReturnCode` columns. Got ",
      if (is.null(bo)) "NULL" else paste(class(bo), collapse = "/"),
      ".",
      call. = FALSE
    )
  }
  bo
}


# --------------------------------------------------------------------------
# RHS section builders
# --------------------------------------------------------------------------

# Generic per-replicate stack -> per-Parameter aggregation. Handles thetas,
# sigmas, and secondaries (omegas have their own builder because diagonals
# need to be filtered via the CI table's `Diagonal` column).
.build_rhs_section <- function(stack, label_col, section,
                               keep, transform, units = list(),
                               metric, ci_level, digits) {
  if (is.null(stack) || !nrow(stack)) {
    return(NULL)
  }
  if (!all(c("Replicate", label_col, "Value") %in% names(stack))) {
    stop(
      "Bootstrap stack for section '", section,
      "' is missing required columns (Replicate, ", label_col, ", Value).",
      call. = FALSE
    )
  }

  filt <- stack[stack$Replicate %in% keep, , drop = FALSE]
  params <- unique(stack[[label_col]])

  rows <- lapply(params, function(name) {
    values <- filt$Value[filt[[label_col]] == name]
    .build_rhs_row(values, name, section, transform, units,
                   metric, ci_level, digits)
  })
  dplyr::bind_rows(rows)
}

# Omega-specific RHS builder: filter to diagonals via BootOmegaCI$Diagonal,
# then aggregate as for the other sections. The diagonal label collapse
# (`<x>_<x>` -> `<x>`) already happened in RsNLME's reader (Plan §2a's
# `.normalizeOmegaDiagonalName`), so the labels here key directly on
# `prmTable$label` / `transform` keys without further transformation.
.build_rhs_omega_section <- function(stack, ci_table,
                                     keep, transform, units = list(),
                                     metric, ci_level, digits) {
  if (is.null(stack) || !nrow(stack)) {
    return(NULL)
  }
  if (!all(c("Replicate", "Omega", "Value") %in% names(stack))) {
    stop(
      "BootOmegaStacked is missing required columns (Replicate, Omega, Value).",
      call. = FALSE
    )
  }

  diag_names <- if (!is.null(ci_table) &&
                    all(c("Omega", "Diagonal") %in% names(ci_table))) {
    ci_table$Omega[.isTRUE_logical(ci_table$Diagonal)]
  } else {
    # Fallback: heuristic when BootOmegaCI is absent or malformed -- keep
    # only labels without the `_` separator that NLME8 uses for off-diagonals.
    # Acceptable because all diagonal labels post-`.normalizeOmegaDiagonalName`
    # match `colnames(dmp.txt$omega)` which are themselves PML names without `_`.
    setdiff(unique(stack$Omega), grep("_", unique(stack$Omega), value = TRUE))
  }

  filt <- stack[stack$Replicate %in% keep & stack$Omega %in% diag_names, ,
                drop = FALSE]
  params <- unique(filt$Omega)

  rows <- lapply(params, function(name) {
    values <- filt$Value[filt$Omega == name]
    .build_rhs_row(values, name, "Random effects", transform, units,
                   metric, ci_level, digits)
  })
  dplyr::bind_rows(rows)
}

# Drop sigma-named rows from `BootThetaStacked` so the Fixed effects
# section doesn't double-count residuals that also live in
# `BootSigmaStacked`. Sigma identification prefers `BootSigmaCI$Sigma`
# (always emitted under Piece 1's always-emit-headers rule, even when
# header-only); falls back to `BootSigmaStacked$Sigma` and finally to
# returning the stack unchanged when neither sigma source is available
# (older NLME8 build pre-§1b: silently degrade rather than hide rows).
.drop_sigma_rows_from_theta_stack <- function(theta_stack, sigma_ci,
                                              sigma_stack) {
  if (is.null(theta_stack) || !nrow(theta_stack)) return(theta_stack)
  sigma_names <- if (!is.null(sigma_ci) && "Sigma" %in% names(sigma_ci) &&
                     nrow(sigma_ci) > 0L) {
    unique(as.character(sigma_ci$Sigma))
  } else if (!is.null(sigma_stack) && "Sigma" %in% names(sigma_stack) &&
             nrow(sigma_stack) > 0L) {
    unique(as.character(sigma_stack$Sigma))
  } else {
    return(theta_stack)
  }
  if (!length(sigma_names)) return(theta_stack)
  theta_stack[!as.character(theta_stack$Theta) %in% sigma_names, , drop = FALSE]
}

# `isTRUE` returns FALSE for logical vectors (>1 element); we want a per-row
# mask. Wrap to an explicit type-coerced version that handles strings like
# "TRUE" too (data.table::fread keeps logicals as logicals, but a raw
# read.csv import could surface them as character).
.isTRUE_logical <- function(x) {
  if (is.logical(x)) return(x %in% TRUE)
  if (is.character(x)) return(toupper(x) %in% c("TRUE", "T"))
  as.logical(x) %in% TRUE
}

# Per-Parameter RHS row: applies the resolved transform to each replicate's
# value (via the same .preset_specs / custom-fn dispatch get_summaryNlme
# uses), then computes the metric and percentile CI on both raw and
# transformed scales when the transform is non-identity.
.build_rhs_row <- function(values, name, section, transform, units = list(),
                           metric, ci_level, digits) {
  ok <- is.finite(values)
  values <- values[ok]

  spec <- .pick_transform(
    list(section = section, label = name),
    transform[[name]]
  )

  fn <- if (spec$preset) .preset_specs[[spec$name]]$fn else spec$fn

  if (!length(values)) {
    raw_metric <- NA_real_
    raw_lo <- NA_real_
    raw_hi <- NA_real_
    transformed_metric <- NA_real_
    transformed_lo <- NA_real_
    transformed_hi <- NA_real_
  } else {
    raw_metric <- if (metric == "Mean") mean(values) else stats::median(values)
    raw_q <- stats::quantile(values,
                             probs = c((1 - ci_level) / 2, 1 - (1 - ci_level) / 2),
                             names = FALSE, na.rm = FALSE)
    raw_lo <- raw_q[1]
    raw_hi <- raw_q[2]

    if (.is_transform_identity(spec)) {
      transformed_metric <- NA_real_
      transformed_lo <- NA_real_
      transformed_hi <- NA_real_
    } else {
      tvals <- fn(values)
      transformed_metric <- if (metric == "Mean") mean(tvals) else stats::median(tvals)
      tq <- stats::quantile(tvals,
                            probs = c((1 - ci_level) / 2, 1 - (1 - ci_level) / 2),
                            names = FALSE, na.rm = FALSE)
      transformed_lo <- tq[1]
      transformed_hi <- tq[2]
    }
  }

  # Transformed-only display: when a non-identity transform is active, show the
  # transformed value/CI only (no `raw (transformed)` pairing); the scale label
  # is appended to the parameter name post-merge, not shown here.
  is_t <- !.is_transform_identity(spec)
  display_metric <- if (is_t) transformed_metric else raw_metric
  display_lo     <- if (is_t) transformed_lo     else raw_lo
  display_hi     <- if (is_t) transformed_hi     else raw_hi

  estimate_cell <- .format_estimate_cell(display_metric, digits)
  ci_cell <- .format_ci_cell(display_lo, display_hi, digits)

  # `Unit` carries real units only, so the bootstrap-side unit comes solely
  # from a user override; the m2/m3 paths drop this column on merge in
  # favour of the richer original-fit unit.
  unit_cell <- if (length(units) && name %in% names(units)) {
    unname(units[[name]])
  } else {
    ""
  }

  out <- tibble::tibble(
    Section = section,
    Parameter = name,
    Unit = unit_cell,
    !!.estimate_col(metric) := estimate_cell,
    !!.ci_col(ci_level)     := ci_cell
  )
  out
}

.is_transform_identity <- function(spec) {
  isTRUE(spec$preset) && identical(spec$name, "raw")
}

.format_estimate_cell <- function(metric_val, digits) {
  if (is.na(metric_val)) return(NA_character_)
  .fmt(metric_val, digits)
}

# Percentile CI rendered as a dash range: "lo - hi".
.format_ci_cell <- function(lo, hi, digits) {
  if (is.na(lo) || is.na(hi)) return(NA_character_)
  paste0(.fmt(lo, digits), " - ", .fmt(hi, digits))
}

.fmt <- function(x, digits) {
  if (is.na(x) || !is.finite(x)) return(format(x))
  format(signif(x, digits = digits))
}

.estimate_col <- function(metric) {
  paste0("Bootstrap estimate (", metric, ")")
}

.ci_col <- function(ci_level) {
  sprintf("Bootstrap %g%% CI", 100 * ci_level)
}


# --------------------------------------------------------------------------
# Soft-degrade for missing per-replicate stacks
# --------------------------------------------------------------------------

# Build NA RHS rows for a section whose per-replicate stack is missing.
# Reference parameter names come from the matching CI table (always
# present per Piece 1's always-emit-headers rule).
.na_rhs_section <- function(param_names, section, metric, ci_level) {
  if (!length(param_names)) return(NULL)
  # `Unit = NA_character_` (rather than "") so the m1 path's "drop empty
  # Unit column" guard treats this section as `unknown` rather than
  # `dimensionless`, matching the section's NA estimate cells.
  tibble::tibble(
    Section = section,
    Parameter = param_names,
    Unit = NA_character_,
    !!.estimate_col(metric) := NA_character_,
    !!.ci_col(ci_level)     := NA_character_
  )
}

.infer_param_names <- function(ci_table, label_col) {
  if (is.null(ci_table) || !nrow(ci_table)) return(character())
  if (!label_col %in% names(ci_table)) return(character())
  diag_mask <- if ("Diagonal" %in% names(ci_table)) {
    .isTRUE_logical(ci_table$Diagonal)
  } else {
    rep(TRUE, nrow(ci_table))
  }
  unique(ci_table[[label_col]][diag_mask])
}

.warn_missing_stacks <- function(missing_stacks) {
  nlme8_ver <- tryCatch(
    as.character(utils::packageVersion("Certara.NLME8")),
    error = function(e) "unknown (Certara.NLME8 not installed)"
  )
  warning(
    "Bootstrap stack(s) missing from `bootResult` (older Certara.NLME8 build): ",
    paste(missing_stacks, collapse = ", "),
    ". Installed Certara.NLME8 version: ", nlme8_ver,
    ". Affected sections show NA in the bootstrap columns. ",
    "Upgrade Certara.NLME8 to populate them.",
    call. = FALSE
  )
}


# --------------------------------------------------------------------------
# LHS + RHS merge and section ordering
# --------------------------------------------------------------------------

.merge_lhs_rhs <- function(lhs, rhs) {
  if (is.null(lhs) || !nrow(lhs)) {
    return(rhs)
  }
  if (is.null(rhs) || !nrow(rhs)) {
    return(lhs)
  }

  # When LHS supplies a Unit column (always true for `get_summaryNlme()`
  # output and for the fitSummary-derived LHS), drop RHS's Unit before
  # joining so we don't end up with `Unit.x` / `Unit.y` collisions. The
  # LHS unit is richer: it merges preset units with per-parameter
  # overrides from the model's structural-parameter units, while the RHS
  # unit comes from the transform spec alone.
  if ("Unit" %in% names(lhs) && "Unit" %in% names(rhs)) {
    rhs <- rhs[, setdiff(names(rhs), "Unit"), drop = FALSE]
  }

  # Use full_join so RHS-only rows (parameters in BS but not in fit) survive.
  result <- dplyr::full_join(lhs, rhs, by = c("Section", "Parameter"))
  result
}

.order_by_section <- function(result) {
  if (is.null(result) || !nrow(result)) return(result)
  levels <- c("Fixed effects", "Random effects", "Residual error", "Secondary")
  ord <- match(result$Section, levels)
  ord[is.na(ord)] <- length(levels) + 1L  # any unknown section sinks to bottom
  result[order(ord, seq_len(nrow(result))), , drop = FALSE]
}


# --------------------------------------------------------------------------
# Confidence-level resolution
# --------------------------------------------------------------------------

# Default `ci_level` inherits the bootstrap run's confidence. RsNLME stores it
# as the attribute `confidenceLevel` on the `rsnlme_boot` object, on the
# percent scale (e.g. 90) or `NA_real_` when unrecoverable. An explicit
# `ci_level` always wins; fall back to 0.95 when metadata is absent.
.resolve_ci_level <- function(bootResult, ci_level) {
  if (!is.null(ci_level)) return(ci_level)
  cl <- attr(bootResult, "confidenceLevel", exact = TRUE)
  if (is.null(cl) || length(cl) != 1L || !is.finite(cl)) {
    return(0.95)
  }
  if (cl > 1) cl <- cl / 100   # 90 -> 0.90
  cl
}


# --------------------------------------------------------------------------
# Scale-flag append (post-merge) and residual-shape advisories
# --------------------------------------------------------------------------

# Append the `(CV%)` / `(SD)` / custom scale flag to the shared `Parameter`
# column once, after the bare-label LHS/RHS join. Re-resolves the spec per row
# from (Section, Parameter, transform) so original-fit, bootstrap, and
# soft-degrade rows are all flagged uniformly.
.append_flags_to_table <- function(result, transform) {
  if (is.null(result) || !nrow(result)) return(result)
  custom_no_name <- character()
  result$Parameter <- vapply(seq_len(nrow(result)), function(i) {
    label <- result$Parameter[i]
    spec <- .pick_transform(
      list(section = result$Section[i], label = label),
      transform[[label]]
    )
    if (!isTRUE(spec$preset) && (is.null(spec$flag) || is.na(spec$flag))) {
      custom_no_name <<- c(custom_no_name, label)
    }
    .append_flag(label, spec)
  }, character(1))
  if (length(custom_no_name)) {
    warning(
      "Custom transform without `name` for: ",
      paste(unique(custom_no_name), collapse = ", "),
      ". Parameter name carries no scale flag; supply `name` to label it.",
      call. = FALSE
    )
  }
  result
}

# Residual-shape advisories, emitted once. Classification needs PML, which is
# only available via `xpdb` (the `rsnlme_boot` object carries none). When
# `pml_code` is NULL/empty, `roles` comes back empty and both advisory
# helpers below fall back to their own "can't tell, so warn" behaviour --
# same as the no-`xpdb` case on the `get_summaryNlme()` path -- rather than
# being skipped here.
.emit_boot_sigma_advisories <- function(sigma_labels, transform, pml_code) {
  if (!length(sigma_labels)) return(invisible())
  roles <- if (!is.null(pml_code) && length(pml_code)) {
    tryCatch(.map_obs_to_sigma(pml_code, sigma_labels)$roles,
             error = function(e) stats::setNames(character(), character()))
  } else {
    stats::setNames(character(), character())
  }

  prm_stub <- data.frame(label = sigma_labels, type = "sig",
                         diagonal = TRUE, stringsAsFactors = FALSE)
  .emit_sigma_role_warnings(prm_stub, transform, roles)
  # Delegate the suppress-only-if-proven-proportional decision to
  # `.emit_residual_default_message()` itself instead of duplicating it here:
  # a sigma absent from `roles` (never matched to an `observe()` block) must
  # count as "not proven proportional", which only the shared helper's
  # per-sigma `roles[s]` lookup gets right (`NA_character_` for a missing
  # name, rather than the `[[` subscript-out-of-bounds error that lookup
  # would raise).
  .emit_residual_default_message(prm_stub, transform, roles)
  invisible()
}