R/get_summaryNlme.R

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

#' @title Build a parameter summary table for an NLME `xpdb`
#'
#' @description Produces a single tibble combining fixed effects, random
#'   effects, residual errors, and secondary parameters, along with their
#'   estimates and `%RSE` on the chosen scale. Shrinkage for random effects
#'   and residual errors is also included. The transform applied to each row
#'   controls both the displayed `Estimate` and the corresponding `%RSE`,
#'   which is computed on the transformed scale via the delta method.
#'   Built-in presets carry analytic derivatives; for transforms outside the
#'   catalog, supply a custom `fn` with its derivative `dfn`.
#'
#' @details Default transforms per section:
#'
#' \itemize{
#'   \item Fixed effects: `raw`.
#'   \item Random effects (omegas): `lognormal_cv`, defined as
#'     \eqn{100 \sqrt{\exp(\omega^2) - 1}}, where \eqn{\omega^2} is the
#'     variance stored in `prmTable$value`.
#'   \item Residual error (sigmas): `multiplicative_cv` (\eqn{100\sigma}).
#'     This default is applied **uniformly to every sigma**, regardless of
#'     the error-model shape declared in PML -- and choosing the wrong
#'     scale doesn't correct itself, it just mislabels the number. For
#'     purely additive or combined error models the %CV label is
#'     misleading -- override with `transform = list(<sigma> = "raw")`
#'     (the engine reports the residual error as a standard deviation) or
#'     a custom `fn`.
#'
#'   When the PML source is available, each sigma's shape is inferred
#'   directly from its `observe()` expression by differentiating it with
#'   respect to the error variable (base R's `stats::D()`) and checking
#'   whether that derivative is a constant (additive error) or
#'   proportional to the noise-free prediction with no curvature
#'   (proportional error -- the one shape `multiplicative_cv` is exact
#'   for). This inference is syntax-only and best-effort: it works
#'   directly on whatever PML the `xpdb` happens to embed, without
#'   assuming it came from any particular model-building tool. The
#'   advisory message is skipped only when every defaulted sigma is
#'   proven proportional; it fires for additive error, for any other
#'   non-proportional shape (combined, power, ...), and whenever the
#'   shape can't be determined at all (a function outside `D()`'s
#'   derivative table, unusual syntax, or no PML source) -- in that last
#'   case the safer default is to warn rather than assume. An extra
#'   per-sigma warning fires for any sigma whose inferred role is
#'   additive or otherwise non-proportional. Set `emit_advisories =
#'   FALSE` to silence both the message and the per-sigma warnings, or
#'   `options(xposeNlme.summary.quiet_default_warning = TRUE)` to
#'   silence only the message. [get_bootSummaryNlme()]'s `bootResult`-only
#'   and embedded-`fitSummary` input modes have no PML to check either, so
#'   they inherit this same "no PML source" behaviour: the default message
#'   always fires for a defaulted sigma, while the per-sigma warning (which
#'   needs a known shape to name) stays silent.
#'   \item Secondary: `raw`.
#' }
#'
#' Built-in preset catalog:
#'
#' \itemize{
#'   \item Fixed effects: `raw` (custom `fn` always available).
#'   \item Random effects: `raw`, `lognormal_cv`, `normal_sd`. `normal_sd`
#'     returns the standard deviation \eqn{\sqrt{\Omega}} (i.e.
#'     `fn = sqrt(value)`), not the variance.
#'   \item Residual error: `raw`, `log_additive_cv`, `multiplicative_cv`.
#'     `log_additive_cv` is \eqn{100 \sqrt{\exp(\sigma^2) - 1}} and
#'     `multiplicative_cv` is \eqn{100\sigma}.
#' }
#'
#' Custom transform spec:
#' `transform = list(<label> = list(fn = function(x) ..., dfn = function(x) ..., name = ...))`.
#' `dfn` is optional; when omitted, %RSE falls back to the raw scale and a
#' single warning per call lists the affected parameters. `name` is the
#' optional scale flag appended to the parameter name; a custom transform
#' supplied without `name` triggers a warning and leaves the name unflagged.
#'
#' Scale flag and units: when a non-identity transform is active the scale
#' label is appended to `Parameter` in parentheses -- `nV (CV%)`,
#' `CEps (SD)`, or the custom `name`. Identity / `raw` rows keep the bare
#' name. The `Unit` column carries physical units only (user `units` or the
#' model's structural-parameter units) and is dropped when every row is
#' dimensionless.
#'
#' Off-diagonal omega/sigma rows are excluded entirely. `transform` and
#' `units` are keyed by `prmTable$label` (the PML-source name -- `tvCl`,
#' `nV`, `CEps`); names not present among the labels emit a single warning
#' per call and are ignored.
#'
#' @param xpdb An `xpose_data` object created by `xposeNlme()` or
#'   `xposeNlmeModel()`.
#' @param .problem Problem number (default 1). Mirrors `get_prmNlme()`.
#' @param .subprob Subproblem number (default 0). Mirrors `get_prmNlme()`.
#' @param .method Estimation method filter (default `NULL`). Mirrors
#'   `get_prmNlme()`.
#' @param transform Named list of per-parameter transforms keyed by
#'   `prmTable$label`. Each value is 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.
#' @param units Named character vector keyed by `prmTable$label` that
#'   populates or overrides the `Unit` column with physical units for
#'   matching rows.
#' @param shrinkage Shrinkage calculation method, one of `"engine"`
#'   (default), `"sd"`, or `"var"`.
#'   \itemize{
#'     \item `"engine"`: uses the standard-deviation-based shrinkage values
#'       reported directly by the engine (read from `xpdb$summary`) --
#'       eta shrinkage \eqn{1 - SD(\eta)/\omega} (with \eqn{\omega =
#'       \sqrt{\Omega}}, the model standard deviation) and eps shrinkage
#'       \eqn{1 - SD(IWRES)}. Note the engine computes the eta SD with
#'       denominator \eqn{n} (population) but the eps SD with denominator
#'       \eqn{n - 1} (sample).
#'     \item `"sd"`: recomputes the standard-deviation-based shrinkage using
#'       R's `sd()` (denominator \eqn{n - 1}). This differs from `"engine"`
#'       only for eta shrinkage, since the engine's eps path already uses
#'       \eqn{n - 1}.
#'     \item `"var"`: recomputes a variance-based shrinkage using R's
#'       `var()` (denominator \eqn{n - 1}) -- eta shrinkage
#'       \eqn{1 - Var(\eta)/\Omega} and eps shrinkage \eqn{1 - Var(IWRES)}.
#'   }
#'   The recompute paths (`"sd"` / `"var"`) use subject-level etas for eta
#'   shrinkage when the model has random effects (skipped for naive-pooled
#'   or any other no-`ranef()` fit) and the IWRES column in `xpdb$data` for
#'   eps shrinkage; multi-residual models additionally need the embedded
#'   PML source to map each `ObsName` row to its driving sigma.
#' @param digits Optional significant-digits count. When non-`NULL`,
#'   `signif()` is applied to the stored numeric `Estimate`, `%RSE`, and
#'   `Shrinkage (%)`, and the `print.summaryNlme` method shows that many
#'   significant figures (by setting `pillar.sigfig` for the duration of
#'   the print). `NULL` (default) keeps full precision in both the stored
#'   values and the printed display.
#' @param append_flag When `TRUE` (default), the scale flag (`CV%` / `SD` /
#'   custom `name`) is appended to `Parameter` for non-identity transforms.
#' @param emit_advisories When `TRUE` (default), the residual-transform
#'   advisory message and the per-sigma type warnings are emitted. Set to
#'   `FALSE` to silence both.
#'
#' @return A tibble (class `summaryNlme`) with columns `Section`,
#'   `Parameter` (carrying a `(CV%)` / `(SD)` / custom scale flag when a
#'   non-identity transform is active), `Estimate`, `%RSE`, `Shrinkage (%)`,
#'   and -- only when at least one row resolves to a non-empty physical
#'   unit -- `Unit`. The `summaryNlme` class carries a `print` method that honours
#'   the `digits` argument for displayed precision.
#'
#' @examples
#' \dontrun{
#' # 1) Default output on a log-normal IIV + proportional error model.
#' xp <- xposeNlmeModel(fit)
#' get_summaryNlme(xp)
#'
#' # 2) Per-parameter override on an additive error model. The engine
#' #    reports residual error as a standard deviation, so `raw` shows the
#' #    SD directly (no misleading %CV flag).
#' get_summaryNlme(
#'   xp,
#'   transform = list(EEps = "raw"),
#'   units     = list(EEps = "ng/mL")
#' )
#'
#' # 3) Custom transform for combined add-mult error.
#' get_summaryNlme(
#'   xp,
#'   transform = list(
#'     CEps = list(
#'       fn  = function(s) 100 * s,
#'       dfn = function(s) 100
#'     )
#'   )
#' )
#' }
#'
#' @seealso [get_prmNlme()], [get_overallNlme()], [get_etaSubjectNlme()]
#' @importFrom magrittr %>%
#' @export
get_summaryNlme <- function(xpdb,
                            .problem = 1,
                            .subprob = 0,
                            .method = NULL,
                            transform = list(),
                            units = list(),
                            shrinkage = c("engine", "sd", "var"),
                            digits = NULL,
                            append_flag = TRUE,
                            emit_advisories = TRUE) {
  xpdb <- .ensure_xpose_data(xpdb)
  shrinkage <- match.arg(shrinkage)

  prm <- .pull_prmTable(xpdb, .problem, .subprob, .method)
  prm <- .drop_offDiagonal(prm)
  prm <- .annotate_section(prm)

  shrink_map <- .resolve_shrinkage_map(xpdb, prm, shrinkage, .problem, .subprob)
  sigma_roles <- .resolve_sigma_roles(xpdb, prm)

  result <- .apply_summary_pipeline(
    prm_like = prm,
    transform = transform,
    units = units,
    shrink_map = shrink_map,
    param_units = xpdb$nlme_param_units,
    sigma_roles = sigma_roles,
    append_flag = append_flag,
    digits = digits,
    emit_advisories = emit_advisories
  )

  .as_summaryNlme(result, digits)
}


# --------------------------------------------------------------------------
# Display class: let `digits` drive printed precision
# --------------------------------------------------------------------------

# Tag the summary tibble so `print.summaryNlme()` can map `digits` onto
# `pillar.sigfig` at display time. The numeric columns are already
# `signif()`-rounded in the pipeline when `digits` is non-NULL; this only
# governs how many significant figures the tibble print method shows.
.as_summaryNlme <- function(x, digits) {
  attr(x, "summary_digits") <- digits
  class(x) <- c("summaryNlme", class(x))
  x
}

#' @rdname get_summaryNlme
#' @param x A `summaryNlme` tibble returned by `get_summaryNlme()`.
#' @param ... Further arguments passed to the underlying tibble print method.
#' @export
print.summaryNlme <- 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)
}


# --------------------------------------------------------------------------
# Shared pipeline: row-build + advisory emissions + Unit-column collapse
# --------------------------------------------------------------------------

# `prm_like` must carry columns: section, type, label, value, se, diagonal.
# Both `get_summaryNlme()` and `get_bootSummaryNlme()`'s fitSummary path feed
# this helper; their inputs differ only in how `prm_like` is constructed
# (xpose's `prmTable` vs. RsNLME's `fitSummary` reshaped into prm_like
# shape). The three advisory emissions live here so any future LHS path
# inherits them automatically -- bypassing them would require deliberately
# not calling the pipeline.
.apply_summary_pipeline <- function(prm_like, transform, units,
                                    shrink_map, param_units,
                                    sigma_roles = NULL,
                                    append_flag = TRUE,
                                    digits = NULL,
                                    emit_advisories = TRUE) {
  transform <- .normalise_transform(transform, prm_like$label)
  units <- .normalise_units(units, prm_like$label)

  rows <- vector("list", nrow(prm_like))
  custom_no_dfn <- character()
  custom_no_name <- character()
  for (i in seq_len(nrow(prm_like))) {
    row <- prm_like[i, ]
    spec <- .pick_transform(row, transform[[row$label]])

    out <- .apply_transform(spec, row$value, row$se, row$label)
    custom_no_dfn <- c(custom_no_dfn, out$missing_dfn_for)
    if (!isTRUE(spec$preset) &&
        (is.null(spec$flag) || is.na(spec$flag))) {
      custom_no_name <- c(custom_no_name, row$label)
    }

    unit <- .resolve_unit(row, spec, units, param_units)

    shrink <- if (length(shrink_map) && row$label %in% names(shrink_map)) {
      unname(shrink_map[[row$label]])
    } else {
      NA_real_
    }

    estimate <- out$estimate
    rse <- out$rse
    if (!is.null(digits)) {
      estimate <- signif(estimate, digits)
      rse <- signif(rse, digits)
      shrink <- if (is.na(shrink)) shrink else signif(shrink, digits)
    }

    rows[[i]] <- tibble::tibble(
      Section = row$section,
      Parameter = if (isTRUE(append_flag)) {
        .append_flag(row$label, spec)
      } else {
        row$label
      },
      Estimate = estimate,
      `%RSE` = rse,
      `Shrinkage (%)` = shrink,
      Unit = unit
    )
  }

  result <- dplyr::bind_rows(rows)
  if (isTRUE(emit_advisories)) {
    .emit_residual_default_message(prm_like, transform, sigma_roles)
    .emit_sigma_role_warnings(prm_like, transform, sigma_roles)
  }
  if (length(custom_no_dfn)) {
    warning(
      "Custom transform without `dfn` for: ",
      paste(unique(custom_no_dfn), collapse = ", "),
      ". %RSE reported on the raw scale.",
      call. = FALSE
    )
  }
  if (isTRUE(append_flag) && 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
    )
  }

  # Drop only when every row is the empty string -- a `NA_character_`
  # entry (e.g. from a user-supplied `units = c(EEps = NA)`) means
  # "unknown" and is preserved, distinct from the dimensionless `""`
  # case. Same predicate used by `get_bootSummaryNlme()` for the
  # merged LHS+RHS table.
  if (!any(is.na(result$Unit) | nzchar(result$Unit))) {
    result$Unit <- NULL
  }
  result
}


# --------------------------------------------------------------------------
# Defensive xpdb coercion
# --------------------------------------------------------------------------

# `xposeNlme()` returns an `xpdb` with class `c("xpose_data", "uneval")`,
# but older ggplot2 (<= 3.x, e.g. R 4.0.x build hosts) ships an
# `[[<-.uneval` method that calls `new_aes(NextMethod())` and forces
# the class to `c("uneval")` -- silently stripping `xpose_data` on any
# `xpdb$x <- y` mutation in user or test code. Subsequent
# `xpose::is.xpdb(xpdb)` then returns FALSE.
#
# Restore the class when the object still looks like an xpdb (canonical
# slots present). Hard-error only when the object is something else
# entirely.
.ensure_xpose_data <- function(xpdb) {
  if (xpose::is.xpdb(xpdb)) return(xpdb)
  if (is.list(xpdb) &&
      all(c("code", "files", "summary", "data") %in% names(xpdb))) {
    class(xpdb) <- c("xpose_data", "uneval")
    return(xpdb)
  }
  stop("`xpdb` is not an xpose_data object.", call. = FALSE)
}


# --------------------------------------------------------------------------
# Section / classification helpers
# --------------------------------------------------------------------------

.section_for <- c(the = "Fixed effects", ome = "Random effects",
                  sig = "Residual error", sec = "Secondary")

.section_default <- c("Fixed effects" = "raw",
                      "Random effects" = "lognormal_cv",
                      "Residual error" = "multiplicative_cv",
                      "Secondary" = "raw")

.section_catalog <- list(
  "Fixed effects" = c("raw"),
  "Random effects" = c("raw", "lognormal_cv", "normal_sd"),
  "Residual error" = c("raw", "log_additive_cv", "multiplicative_cv"),
  "Secondary" = c("raw")
)

.pull_prmTable <- function(xpdb, .problem, .subprob, .method) {
  rows <- xpdb$files %>%
    dplyr::filter(name == "prmTable" &
                    problem == get(".problem") &
                    subprob == get(".subprob"))
  if (!is.null(.method)) {
    rows <- dplyr::filter(rows, .data$method == get(".method"))
  }
  if (nrow(rows) == 0) {
    stop("No prmTable for problem ", .problem,
         " / subprob ", .subprob, ".", call. = FALSE)
  }
  if (nrow(rows) > 1) {
    warning("More than one prmTable matches; using the last.",
            call. = FALSE)
    rows <- rows[nrow(rows), ]
  }
  rows$data[[1]]
}

.drop_offDiagonal <- function(prm) {
  drop <- prm$type %in% c("ome", "sig") &
    !is.na(prm$diagonal) & !prm$diagonal
  prm[!drop, , drop = FALSE]
}

.annotate_section <- function(prm) {
  prm$section <- unname(.section_for[prm$type])
  prm
}


# --------------------------------------------------------------------------
# Transform argument validation and dispatch
# --------------------------------------------------------------------------

.normalise_transform <- function(transform, labels) {
  if (!is.list(transform)) {
    stop("`transform` must be a named list.", call. = FALSE)
  }
  if (!length(transform)) return(list())
  if (is.null(names(transform)) || any(!nzchar(names(transform)))) {
    stop("Every entry of `transform` must be named.", call. = FALSE)
  }
  unknown <- setdiff(names(transform), labels)
  if (length(unknown)) {
    warning("Ignored `transform` keys not in prmTable$label: ",
            paste(unknown, collapse = ", "), ".", call. = FALSE)
    transform <- transform[setdiff(names(transform), unknown)]
  }
  transform
}

.normalise_units <- function(units, labels) {
  if (!length(units)) return(stats::setNames(character(), character()))
  if (is.list(units)) units <- unlist(units)
  if (!is.character(units) || is.null(names(units))) {
    stop("`units` must be a named character vector or list.", call. = FALSE)
  }
  unknown <- setdiff(names(units), labels)
  if (length(unknown)) {
    warning("Ignored `units` keys not in prmTable$label: ",
            paste(unknown, collapse = ", "), ".", call. = FALSE)
    units <- units[setdiff(names(units), unknown)]
  }
  units
}

# Resolve the transform spec for a single row: user override, then section
# default. Validate preset against the section's catalog; pass custom specs
# through.
.pick_transform <- function(row, user_spec) {
  if (is.null(user_spec)) {
    return(list(name = unname(.section_default[row$section]),
                preset = TRUE,
                section = row$section))
  }

  if (is.character(user_spec) && length(user_spec) == 1L) {
    catalog <- .section_catalog[[row$section]]
    if (!user_spec %in% catalog) {
      stop("Transform '", user_spec, "' for '", row$label,
           "' is not in the ", row$section, " catalog (",
           paste(catalog, collapse = ", "), ").", call. = FALSE)
    }
    return(list(name = user_spec, preset = TRUE, section = row$section))
  }

  if (is.list(user_spec) && is.function(user_spec$fn)) {
    return(list(name = "custom", preset = FALSE,
                section = row$section,
                fn = user_spec$fn,
                dfn = if (is.function(user_spec$dfn)) user_spec$dfn else NULL,
                flag = if (is.character(user_spec$name) &&
                           length(user_spec$name) == 1L &&
                           nzchar(user_spec$name)) user_spec$name
                       else NA_character_))
  }

  stop("Invalid transform spec for '", row$label,
       "': must be a preset string or list(fn = ...) ",
       "(optional dfn, name).",
       call. = FALSE)
}


# --------------------------------------------------------------------------
# Estimate / %RSE computation (delta method on the transformed scale)
# --------------------------------------------------------------------------

# Each preset declares fn(value) and dfn(value). The delta method gives
# rse_t = |dfn(value) / fn(value)| * se * 100. raw / multiplicative_cv
# collapse to the same numeric rse as the untransformed value (linear or
# identity scale change).
# `flag` is the scale label appended to the parameter name (e.g. `nV (CV%)`),
# not a physical unit. The `Unit` column carries real units only. `raw` is the
# only identity transform and carries no flag.
.preset_specs <- list(
  raw = list(
    fn  = function(x) x,
    dfn = function(x) 1,
    flag = ""
  ),
  lognormal_cv = list(
    fn  = function(x) 100 * sqrt(exp(x) - 1),
    dfn = function(x) 50 * exp(x) / sqrt(exp(x) - 1),
    flag = "CV%"
  ),
  normal_sd = list(
    fn  = function(x) sqrt(x),
    dfn = function(x) 1 / (2 * sqrt(x)),
    flag = "SD"
  ),
  log_additive_cv = list(
    fn  = function(x) 100 * sqrt(exp(x^2) - 1),
    dfn = function(x) 100 * x * exp(x^2) / sqrt(exp(x^2) - 1),
    flag = "CV%"
  ),
  multiplicative_cv = list(
    fn  = function(x) 100 * x,
    dfn = function(x) 100,
    flag = "CV%"
  )
)

# Scale flag for a resolved transform spec. Presets read `.preset_specs$flag`;
# custom (`list(fn=, dfn=, name=)`) specs use the user-supplied `name` and
# return "" when none was given.
.transform_flag <- function(spec) {
  if (isTRUE(spec$preset)) {
    f <- .preset_specs[[spec$name]]$flag
    if (is.null(f)) "" else f
  } else {
    if (is.null(spec$flag) || is.na(spec$flag)) "" else spec$flag
  }
}

# Append the scale flag to a parameter label in parentheses, e.g.
# `nV` + `CV%` -> `nV (CV%)`. Identity/`raw` rows (empty flag) stay bare.
.append_flag <- function(label, spec) {
  f <- .transform_flag(spec)
  if (nzchar(f)) paste0(label, " (", f, ")") else label
}

.apply_transform <- function(spec, value, se, label) {
  if (spec$preset) {
    p <- .preset_specs[[spec$name]]
    est <- p$fn(value)
    rse <- if (is.na(se)) NA_real_ else .delta_rse(p$dfn(value), se, est)
    return(list(estimate = est, rse = rse, missing_dfn_for = character()))
  }

  est <- spec$fn(value)
  if (is.null(spec$dfn)) {
    # Fallback: report %RSE on the raw scale -- this is exactly the
    # delta-method form with derivative 1, so reuse `.delta_rse()` to
    # pick up the same NA / zero / non-finite guards that the preset
    # path enforces. Avoids `100 * se / 0 = Inf` slipping through to
    # the output tibble for frozen-at-zero or pathological values.
    rse <- if (is.na(se)) NA_real_ else .delta_rse(1, se, value)
    return(list(estimate = est, rse = rse, missing_dfn_for = label))
  }
  rse <- if (is.na(se)) NA_real_ else .delta_rse(spec$dfn(value), se, est)
  list(estimate = est, rse = rse, missing_dfn_for = character())
}

.delta_rse <- function(deriv, se, est) {
  if (is.na(est) || est == 0 || !is.finite(est)) return(NA_real_)
  100 * abs(deriv * se / est)
}


# --------------------------------------------------------------------------
# Shrinkage source dispatch
# --------------------------------------------------------------------------

.resolve_shrinkage_map <- function(xpdb, prm, mode, .problem, .subprob) {
  if (mode == "engine") {
    return(.parse_engine_shrinkage(xpdb, .problem, .subprob))
  }

  eta_shrink <- .recompute_eta_shrinkage(xpdb, prm, mode, .problem)
  eps_shrink <- .recompute_eps_shrinkage(xpdb, prm, mode, .problem)

  c(eta_shrink, eps_shrink)
}

.parse_engine_shrinkage <- function(xpdb, .problem, .subprob) {
  pull <- function(target) {
    s <- xpdb$summary
    hit <- s$problem == .problem & s$subprob == .subprob & s$label == target
    if (!any(hit)) return(NULL)
    .parse_named_value_string(s$value[which(hit)[1]])
  }
  out <- c(pull("etashk"), pull("epsshk"))
  if (!length(out)) return(stats::setNames(numeric(), character()))
  100 * out
}

# Engine values land in xpdb$summary$value as "name1 = v1, name2 = v2, ...".
# Returns a named numeric. Silently drops malformed entries.
.parse_named_value_string <- function(s) {
  if (is.null(s) || is.na(s) || !nzchar(s)) {
    return(stats::setNames(numeric(), character()))
  }
  parts <- trimws(strsplit(s, ",", fixed = TRUE)[[1]])
  m <- regmatches(parts, regexec("^([^=]+?)\\s*=\\s*(.+)$", parts))
  ok <- vapply(m, function(x) length(x) == 3L, logical(1))
  if (!any(ok)) return(stats::setNames(numeric(), character()))
  m <- m[ok]
  nms <- vapply(m, function(x) trimws(x[[2L]]), character(1))
  vals <- suppressWarnings(as.numeric(vapply(m, function(x) trimws(x[[3L]]),
                                             character(1))))
  stats::setNames(vals, nms)
}

# Sample-SD / variance-ratio recompute on the eta_subject table using R's
# sd() / var() (denominator n-1). The engine's eta shrinkage uses
# denominator n (population variance; see nlme-engine .../writelog.f90
# lines ~528-540), so this recompute intentionally differs from the
# "engine" mode for etas. (The engine's eps path does use n-1, so the eps
# recompute matches there.)
.recompute_eta_shrinkage <- function(xpdb, prm, mode, .problem) {
  eta_sub <- tryCatch(get_etaSubjectNlme(xpdb, .problem),
                      error = function(e) NULL)
  if (is.null(eta_sub)) return(stats::setNames(numeric(), character()))

  ome <- prm[prm$type == "ome" & prm$diagonal == TRUE, ]
  ome_var <- stats::setNames(ome$value, ome$label)

  etas <- intersect(unique(eta_sub$Eta), names(ome_var))
  out <- vapply(etas, function(eta) {
    vals <- eta_sub$ETA_VAL[eta_sub$Eta == eta]
    if (length(vals) < 2L) return(NA_real_)
    if (mode == "sd") {
      1 - stats::sd(vals) / sqrt(ome_var[[eta]])
    } else {
      1 - stats::var(vals) / ome_var[[eta]]
    }
  }, numeric(1))
  100 * out
}

# Per-sigma IWRES pooling. Reproduces the FORTRAN engine algorithm in
# nlme-engine .../writelog.f90 lines 779-806 in R. For single-sigma
# models every IWRES row pools to the lone sigma and no PML map is
# needed; only multi-sigma models require the ObsName -> sigma map
# from `.map_obs_to_sigma()`.
.recompute_eps_shrinkage <- function(xpdb, prm, mode, .problem) {
  sigma_labels <- prm$label[prm$type == "sig" & prm$diagonal]
  if (!length(sigma_labels)) return(stats::setNames(numeric(), character()))

  data_row <- xpdb$data[xpdb$data$problem == .problem, ]
  if (!nrow(data_row)) {
    stop("No data found for problem ", .problem, ".", call. = FALSE)
  }
  d <- data_row$data[[1]]
  if (!"IWRES" %in% names(d)) {
    stop("IWRES not in residuals; use shrinkage = \"engine\".",
         call. = FALSE)
  }

  pool <- if (length(sigma_labels) == 1L) {
    # Single sigma: every IWRES row trivially maps to it.
    list(d$IWRES)
  } else {
    # Multi-sigma: need the ObsName -> sigma map from the embedded PML.
    if (is.null(xpdb$code) || !length(xpdb$code)) {
      stop("xpdb$code is empty; use shrinkage = \"engine\".",
           call. = FALSE)
    }
    if (!"ObsName" %in% names(d)) {
      stop("ObsName not in data; use shrinkage = \"engine\".",
           call. = FALSE)
    }

    m <- .map_obs_to_sigma(xpdb$code, sigma_labels)$map
    if (!length(m)) {
      stop("Could not recover ObsName -> sigma map; ",
           "use shrinkage = \"engine\".", call. = FALSE)
    }

    # Strip a trailing parenthetical unit suffix the engine sometimes
    # writes into ObsName (e.g. `CObs(ng/mL)` when the input dataset
    # carried a `#@` units row). The PML map is keyed by the bare
    # observe name (`CObs`); without this normalisation the lookup
    # silently returns NA on every row and the per-sigma shrinkage
    # collapses to NA without any visible failure.
    obs_key <- sub("\\([^)]*\\)\\s*$", "", as.character(d$ObsName))
    d_sigma <- unname(m[obs_key])
    lapply(sigma_labels, function(sig) d$IWRES[d_sigma == sig])
  }

  out <- vapply(seq_along(sigma_labels), function(i) {
    vals <- pool[[i]]
    vals <- vals[is.finite(vals)]
    if (length(vals) < 2L) return(NA_real_)
    if (mode == "sd") 1 - stats::sd(vals) else 1 - stats::var(vals)
  }, numeric(1))
  stats::setNames(100 * out, sigma_labels)
}


# --------------------------------------------------------------------------
# Unit resolution (priority chain) and conditional column emission
# --------------------------------------------------------------------------

# Resolve the real unit for a row. Priority: user `units` override, then the
# model's structural-parameter units, else dimensionless "".
.resolve_unit <- function(row, spec, user_units, param_units) {
  if (!is.null(user_units) && row$label %in% names(user_units)) {
    return(unname(user_units[[row$label]]))
  }

  if (!is.null(param_units) && row$label %in% names(param_units)) {
    u <- unname(param_units[[row$label]])
    if (!is.na(u) && nzchar(u)) return(u)
  }

  ""
}


# --------------------------------------------------------------------------
# Sigma-role classification and warnings
# --------------------------------------------------------------------------

.resolve_sigma_roles <- function(xpdb, prm) {
  sigma_labels <- prm$label[prm$type == "sig" & prm$diagonal]
  if (!length(sigma_labels) ||
      is.null(xpdb$code) || !length(xpdb$code)) {
    return(stats::setNames(character(), character()))
  }
  tryCatch(
    .map_obs_to_sigma(xpdb$code, sigma_labels)$roles,
    error = function(e) stats::setNames(character(), character())
  )
}

.emit_residual_default_message <- function(prm, transform, roles = NULL) {
  sigma_labels <- prm$label[prm$type == "sig" & prm$diagonal]
  defaulted <- setdiff(sigma_labels, names(transform))
  if (!length(defaulted)) return(invisible())
  if (isTRUE(getOption("xposeNlme.summary.quiet_default_warning"))) {
    return(invisible())
  }
  # Suppress the advisory only when every defaulted sigma's role has been
  # symbolically *proven* proportional (see `.classify_sigma_role()`) --
  # the one shape `multiplicative_cv` is exact for. Keep emitting for
  # additive, combined/power/other, mixed, unknown, or when `roles` is
  # empty (no PML source to check): in all of those cases assuming
  # proportional error would be a guess, and the safer default is to warn.
  if (length(roles)) {
    # Single-bracket lookup: `roles` is a plain named character vector, so
    # `roles[[s]]` would error ("subscript out of bounds") for a sigma with
    # no matching observe() block instead of the list-like `NULL` a reader
    # might expect. `roles[s]` returns `NA_character_` for a missing name,
    # which the `is.na()` check below correctly treats as "not proven
    # proportional".
    role_is_proportional <- vapply(defaulted, function(s) {
      r <- unname(roles[s])
      !is.na(r) && identical(r, "proportional")
    }, logical(1))
    if (all(role_is_proportional)) return(invisible())
  }
  message(
    "Default residual transform is `multiplicative_cv`, appropriate for a ",
    "genuinely proportional error model. If your model is additive, ",
    "combined, or otherwise non-proportional, override the affected ",
    "sigma(s) with `transform = list(<sigma> = ...)` (see ?get_summaryNlme)."
  )
  invisible()
}

.emit_sigma_role_warnings <- function(prm, transform, roles) {
  if (!length(roles)) return(invisible())
  sigma_labels <- prm$label[prm$type == "sig" & prm$diagonal]
  for (sig in sigma_labels) {
    if (sig %in% names(transform)) next
    # Single-bracket lookup (see `.emit_residual_default_message()`): `[[`
    # would error for a sigma absent from `roles` instead of yielding `NA`.
    role <- unname(roles[sig])
    if (is.na(role) || !nzchar(role) || identical(role, "unknown")) next
    if (identical(role, "additive")) {
      warning(
        "Sigma `", sig, "` looks additive in PML, but is reported with ",
        "the `multiplicative_cv` default. Consider ",
        "`transform = list(", sig, " = \"raw\")` ",
        "(the engine reports the residual error as a standard deviation).",
        call. = FALSE
      )
    } else if (role %in% c("other", "mixed")) {
      warning(
        "Sigma `", sig, "` does not look proportional in PML, but is ",
        "reported with the `multiplicative_cv` default (100 * sigma). ",
        "Review the error model and supply a matching ",
        "`transform = list(", sig, " = ...)` if this is misleading ",
        "(see ?get_summaryNlme).",
        call. = FALSE
      )
    }
  }
  invisible()
}