R/compare_prmNlme.R

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

#' @title Compare NLME parameter estimates across multiple runs
#'
#' @description Builds a single wide table comparing parameter estimates
#'   (and optional `%RSE`) across two or more NLME runs, alongside a block of
#'   run-level diagnostics (`-2LL`, `OFV diff`, `method`, `RetCode`,
#'   `condition`, `condition basis`, `nSub`, `nObs`, and total runtime). It is the multi-model
#'   sibling of [get_summaryNlme()]: where `get_summaryNlme()` summarises one
#'   `xpose_data` object, `compare_prmNlme()` lines several up side by side for
#'   run-record style model comparison.
#'
#' @details Runs can be supplied three ways:
#'
#' \itemize{
#'   \item As a pre-built named `list` of `xpose_data` objects (or a single
#'     `xpose_data`) via `x`. List names become the column labels.
#'   \item By explicit run name via `runs`, loaded from `dir` with
#'     [xposeNlme()].
#'   \item By auto-detection (`auto_detect = TRUE`, the default when `x` and
#'     `runs` are both `NULL`): `dir` is scanned for model files (`*.mdl` or
#'     `*.mmdl`, case-insensitive) whose same-named run output folder exists,
#'     in alphanumeric order. Name each model so its output folder matches the
#'     model file (e.g. `run001.mmdl` beside a `run001/` folder). A run is
#'     only included when its `nlme7engine.log` is present; runs that fail to
#'     load are recorded in `log_file` (written into `dir`) and skipped.
#' }
#'
#' The `transform` setting affects only diagonal OMEGA rows. Diagonal SIGMA
#' rows are kept on the reported [get_prmNlme()] scale (for `CEps` this is the
#' SD scale), so SIGMA values and `%RSE` do not change across transformation
#' settings. `%RSE` on transformed OMEGA is propagated with the delta method.
#'
#' The returned object carries `n_header` / `n_rse` attributes marking the
#' leading diagnostic block and the trailing RSE block, which
#' [as_flextable.prmComparisonNlme()] uses to draw separators. Rendering with
#' `flextable` and CSV export via `output_file` are both optional; the core
#' computation depends only on packages already imported by
#' `Certara.Xpose.NLME`.
#'
#' ## Timing
#'
#' The diagnostic row `total runtime (sec)` is engine-reported **CPU** time:
#' the sum of the `runtime` and `covtime` rows in `xpdb$summary`, which are
#' parsed from `nlme7engine.log` at import. That total can differ
#' substantially from the wall-clock elapsed time shown by
#' `print.rsnlme_fit` / a fit object's `runTime` (especially on multi-core
#' runs). See also [get_overallNlme()] for the same distinction, including
#' optional `runtime_wallclock` when an xpdb was built via
#' [xposeNlmeModel()].
#'
#' @param x Optional named `list` of `xpose_data` objects (or a single
#'   `xpose_data`). When supplied, `dir` / `runs` / `auto_detect` are ignored.
#'   List names become the column headers and must be unique.
#' @param dir Directory scanned for runs when `x` is `NULL` (default the
#'   working directory).
#' @param runs Optional character vector of run names (subfolders of `dir`)
#'   to load explicitly, in the given order. Overrides auto-detection.
#' @param auto_detect Logical; when `x` and `runs` are `NULL`, scan `dir` for
#'   completed runs (default `TRUE`).
#' @param max_runs Optional cap on the number of **successfully loaded** runs
#'   to include. `NULL` or `""` includes all runs. Failed imports do not count
#'   toward the cap.
#' @param transform One of `"untransformed"` (default), `"sqrt_om2"`, or
#'   `"sqrt_exp_om2_minus_1"`; applied to diagonal OMEGA only.
#' @param param_order One of `"original"` (default: first run's model order,
#'   with parameters unique to later runs appended in those runs' order) or
#'   `"alphabetical"`.
#' @param rse_separate Logical; when `TRUE`, `%RSE` is shown in its own set of
#'   rows (suffixed `" (RSE)"`) rather than in-line with each estimate.
#'   Accepts the strings `"YES"`/`"NO"` for backward compatibility.
#' @param drop_dOFV Logical; when `TRUE`, omit the `OFV diff` row. Accepts
#'   `"YES"`/`"NO"`.
#' @param output_file Optional path; when set, the full table is written as a
#'   CSV in the chosen `format`.
#' @param format One of `"column"` (default, models as columns) or `"row"`
#'   (models as rows, transposed). Controls the CSV layout and is stored on
#'   the returned object so [as_flextable.prmComparisonNlme()] defaults to the
#'   same orientation.
#' @param log_file Name of the excluded-run log written into `dir` during
#'   auto-detection (default `"Table_log.txt"`). Set to `NULL` to disable.
#'
#' @return A tibble of class `prmComparisonNlme` with a `Description` column
#'   followed by one column per run, carrying `n_header`, `n_rse`, `nModels`,
#'   `run_labels`, and `format` attributes.
#'
#' @examples
#' \dontrun{
#' # 1) Compare two already-imported runs.
#' xp1 <- xposeNlme(dir = "run001")
#' xp2 <- xposeNlme(dir = "run002")
#' compare_prmNlme(list(run001 = xp1, run002 = xp2))
#'
#' # 2) Auto-detect every completed run under the working directory,
#' #    put %RSE on its own rows, and write a CSV.
#' tbl <- compare_prmNlme(
#'   rse_separate = TRUE,
#'   output_file  = "TableofParameters.csv"
#' )
#'
#' # 3) Render the comparison as a flextable (requires the flextable and
#' #    officer packages).
#' flextable::as_flextable(tbl)
#' }
#'
#' @seealso [get_summaryNlme()], [get_prmNlme()], [get_overallNlme()],
#'   [xposeNlme()]
#' @importFrom magrittr %>%
#' @export
compare_prmNlme <- function(x = NULL,
                            dir = ".",
                            runs = NULL,
                            auto_detect = TRUE,
                            max_runs = NULL,
                            transform = c("untransformed", "sqrt_om2",
                                          "sqrt_exp_om2_minus_1"),
                            param_order = c("original", "alphabetical"),
                            rse_separate = FALSE,
                            drop_dOFV = FALSE,
                            output_file = NULL,
                            format = c("column", "row"),
                            log_file = "Table_log.txt") {
  transform <- match.arg(transform)
  param_order <- match.arg(param_order)
  format <- match.arg(format)
  rse_separate <- .as_logical_flag(rse_separate)
  drop_dOFV <- .as_logical_flag(drop_dOFV)

  xpdb_list <- .resolve_xpdb_list(
    x = x, dir = dir, runs = runs, auto_detect = auto_detect,
    max_runs = max_runs, log_file = log_file
  )

  tbl <- .build_prmComparison(
    xpdb_list,
    transform = transform,
    param_order = param_order,
    rse_separate = rse_separate,
    drop_dOFV = drop_dOFV
  )
  attr(tbl, "format") <- format

  if (!is.null(output_file)) {
    .write_prmComparison(tbl, output_file, format = format)
  }

  tbl
}


# --------------------------------------------------------------------------
# Input resolution: pre-built list, explicit runs, or auto-detection
# --------------------------------------------------------------------------

.resolve_xpdb_list <- function(x, dir, runs, auto_detect, max_runs, log_file) {
  n_max <- .resolve_max_runs(max_runs)

  if (!is.null(x)) {
    if (inherits(x, "xpose_data")) {
      xpdb_list <- list(x)
    } else if (is.list(x)) {
      is_xpdb <- vapply(x, function(e) inherits(e, "xpose_data"), logical(1))
      if (!all(is_xpdb)) {
        stop("`x` must be an xpose_data object or a list of xpose_data objects.",
             call. = FALSE)
      }
      xpdb_list <- x
    } else {
      stop("`x` must be an xpose_data object or a list of xpose_data objects.",
           call. = FALSE)
    }
    xpdb_list <- .name_xpdb_list(xpdb_list)
    if (!is.na(n_max) && length(xpdb_list) > n_max) {
      xpdb_list <- xpdb_list[seq_len(n_max)]
    }
    if (length(xpdb_list) < 1) {
      stop("No runs found to summarise.", call. = FALSE)
    }
    return(xpdb_list)
  }

  run_names <- .detect_runs(dir, runs, auto_detect)
  if (!length(run_names)) {
    stop("No runs found to summarise - check `dir`, `runs`, or `auto_detect`.",
         call. = FALSE)
  }
  .load_runs(dir, run_names, log_file,
             max_runs = n_max, explicit_runs = !is.null(runs))
}

# Ensure every run has a usable, unique column label.
.name_xpdb_list <- function(xpdb_list) {
  nm <- names(xpdb_list)
  if (is.null(nm) || any(is.na(nm)) || any(!nzchar(nm))) {
    names(xpdb_list) <- paste0("run", seq_along(xpdb_list))
  }
  .assert_unique_labels(names(xpdb_list))
  xpdb_list
}

# Names become the table's column headers, so duplicates would produce
# ambiguous columns (and a corrupt wide pivot). Error clearly instead,
# matching the MCP `xpose_compare_params` handler.
.assert_unique_labels <- function(labels) {
  if (anyDuplicated(labels)) {
    dups <- unique(labels[duplicated(labels)])
    stop("Model names must be unique (they become column headers): ",
         paste(dups, collapse = ", "), call. = FALSE)
  }
  invisible(labels)
}

# Interpret max_runs (NULL / "" / NA -> no cap; else a positive integer).
.resolve_max_runs <- function(max_runs) {
  if (is.null(max_runs) || length(max_runs) == 0) return(NA_integer_)
  if (length(max_runs) == 1 &&
      (is.na(max_runs) || identical(as.character(max_runs), ""))) {
    return(NA_integer_)
  }
  n <- suppressWarnings(as.integer(max_runs))
  if (is.na(n) || n < 1) {
    stop("`max_runs` must be a positive whole number, NULL, or blank.",
         call. = FALSE)
  }
  n
}

# Determine which run names to load: explicit `runs`, else auto-detect model
# files (`*.mdl` or `*.mmdl`, case-insensitive) whose same-named output folder
# exists (in alphanumeric order). RsNLME's textual/metamodel workflow writes
# `.mmdl` files, so both extensions must be recognised.
.detect_runs <- function(dir, runs, auto_detect) {
  if (!is.null(runs)) {
    rn <- as.character(runs)
    .assert_unique_labels(rn)
  } else if (isTRUE(auto_detect)) {
    mdl <- list.files(path = dir, pattern = "\\.mm?dl$", ignore.case = TRUE)
    run_names <- sort(unique(sub("\\.mm?dl$", "", mdl, ignore.case = TRUE)))
    rn <- character(0)
    for (r in run_names) {
      if (dir.exists(file.path(dir, r))) rn <- c(rn, r)
    }
  } else {
    stop("No runs specified: supply `x`, `runs`, or set `auto_detect = TRUE`.",
         call. = FALSE)
  }
  rn[!is.na(rn) & nzchar(rn)]
}

# Load each run via xposeNlme(); a run is only kept when its
# nlme7engine.log exists and import succeeds. Excluded runs are recorded in
# `log_file` (written into `dir`); a stale log written by a previous call is
# removed when there are no exclusions this call (only if the file still
# carries our header). `max_runs` caps successful loads, not candidates.
# Explicit `runs` exclusions are warned; auto-detect exclusions are messaged.
.load_runs <- function(dir, run_names, log_file,
                       max_runs = NA_integer_, explicit_runs = FALSE) {
  xpdb_built <- list()
  excluded <- character(0)
  reasons <- character(0)

  for (rn in run_names) {
    if (!is.na(max_runs) && length(xpdb_built) >= max_runs) break

    run_dir <- file.path(dir, rn)
    log_path <- file.path(run_dir, "nlme7engine.log")
    if (!file.exists(log_path)) {
      excluded <- c(excluded, rn)
      reasons <- c(reasons,
        "missing nlme7engine.log (run failed, incomplete, or folder not found)")
      next
    }
    xpdb_try <- tryCatch(xposeNlme(dir = run_dir), error = function(e) e)
    if (inherits(xpdb_try, "error")) {
      excluded <- c(excluded, rn)
      reasons <- c(reasons,
        paste0("xposeNlme() error: ", conditionMessage(xpdb_try)))
      next
    }
    xpdb_built[[rn]] <- xpdb_try
  }

  included <- names(xpdb_built)
  log_target <- if (is.character(log_file) && length(log_file) == 1 &&
                    nzchar(log_file)) file.path(dir, log_file) else NULL

  if (length(excluded) > 0 && !is.null(log_target)) {
    log_lines <- c(
      paste0("Table of Parameters - runs excluded (", format(Sys.time()), ")"),
      paste0("Included ", length(included), " of ", length(run_names), " run(s)."),
      "",
      sprintf("%-24s %s", "Run", "Reason"),
      sprintf("%-24s %s", excluded, reasons)
    )
    writeLines(log_lines, log_target)
    msg <- sprintf("%d run(s) excluded (see %s): %s",
                   length(excluded), log_target,
                   paste(excluded, collapse = ", "))
    if (isTRUE(explicit_runs)) {
      warning(msg, call. = FALSE)
    } else {
      message(msg)
    }
  } else if (!is.null(log_target) && file.exists(log_target)) {
    first <- tryCatch(readLines(log_target, n = 1L, warn = FALSE),
                      error = function(e) character(0))
    if (length(first) &&
        grepl("^Table of Parameters - runs excluded", first[[1]])) {
      file.remove(log_target)
    }
  }

  if (length(included) < 1) {
    stop("No runs with a valid nlme7engine.log were found - nothing to ",
         "summarise.", if (!is.null(log_target)) paste0(" See ", log_target, "."),
         call. = FALSE)
  }

  xpdb_built
}


# --------------------------------------------------------------------------
# Table construction
# --------------------------------------------------------------------------

.build_prmComparison <- function(xpdb_list, transform, param_order,
                                 rse_separate, drop_dOFV) {
  run_labels <- names(xpdb_list)

  param_long <- dplyr::bind_rows(lapply(run_labels, function(nm) {
    .prep_params_xpnlme(xpdb_list[[nm]], nm, transform)
  }))

  # Canonical (Section, Description) order for param_order = "original":
  # seed from the first run, then append each later run's new parameters in
  # that run's own order (avoids alphabetical tiebreaks across heterogeneous
  # parameter sets).
  canon_ord <- .canonical_param_ord(param_long, run_labels)

  # Collapse + pivot a chosen value column to one wide row per
  # (Section, Description). Shared by the estimate table and the optional
  # RSE table so both use identical collapsing, pivoting and ordering.
  make_wide <- function(value_col) {
    pc <- param_long %>%
      dplyr::group_by(Section, Description, run) %>%
      dplyr::summarise(
        value_str = .first_nonNA(.data[[value_col]]),
        .groups = "drop"
      )

    w <- .pivot_runs_wide(
      pc[, c("Section", "Description", "run", "value_str")], run_labels
    )
    w <- dplyr::left_join(w, canon_ord, by = c("Section", "Description"))

    sec_rank <- factor(w$Section, levels = c("TH", "OM", "SI"))
    ord <- if (param_order == "alphabetical") {
      order(sec_rank, w$Description)
    } else {
      order(sec_rank, w$param_ord)
    }
    w <- w[ord, , drop = FALSE]
    w$param_ord <- NULL
    w[, c("Section", "Description", run_labels), drop = FALSE]
  }

  # Estimate table: value only when RSE is separated out, else value (%RSE).
  out <- make_wide(if (rse_separate) "disp_val" else "display")

  # Optional RSE block: same parameters, same order, names suffixed " (RSE)".
  rse_block <- NULL
  if (rse_separate) {
    rse_block <- make_wide("disp_rse")
    rse_block$Description <- paste0(rse_block$Description, " (RSE)")
  }

  out_cols <- colnames(out)

  overall <- lapply(xpdb_list, function(x) get_overallNlme(x))
  num_field <- function(field) {
    vapply(overall, function(o) {
      v <- o[[field]]
      if (is.null(v) || length(v) == 0) return(NA_real_)
      as.numeric(v[[1]])
    }, numeric(1))
  }
  chr_field <- function(field) {
    vapply(overall, function(o) {
      v <- o[[field]]
      if (is.null(v) || length(v) == 0 || is.na(v[[1]])) return(NA_character_)
      as.character(v[[1]])
    }, character(1))
  }

  # Format helpers: NA -> blank; -2LL / OFV keep fixed decimals for diffs.
  fmt_num <- function(x, fmt = "%g") {
    ifelse(is.na(x), "", sprintf(fmt, x))
  }

  ofv_vals <- num_field("-2LL")
  ofv_row  <- .diag_row("-2LL", fmt_num(ofv_vals, "%.3f"), run_labels, out_cols)

  dofv_vals <- ofv_vals - ofv_vals[[1]]
  dofv_str <- fmt_num(dofv_vals, "%.3f")
  dofv_str[1] <- ""
  dofv_row <- .diag_row("OFV diff", dofv_str, run_labels, out_cols)

  met_vals <- vapply(xpdb_list, .compare_method, character(1))
  met_row  <- .diag_row("method", met_vals, run_labels, out_cols)

  ret_row  <- .diag_row("RetCode", fmt_num(num_field("RetCode")),
                        run_labels, out_cols)
  cn_vals  <- num_field("Condition")
  cn_row   <- .diag_row("condition", fmt_num(cn_vals), run_labels, out_cols)

  basis_vals <- chr_field("ConditionBasis")
  basis_str  <- ifelse(is.na(basis_vals), "", basis_vals)
  basis_row  <- .diag_row("condition basis", basis_str, run_labels, out_cols)
  basis_unique <- unique(basis_vals[!is.na(basis_vals) & nzchar(basis_vals)])
  if (length(basis_unique) > 1L) {
    warning("Condition bases differ across runs (",
            paste(basis_unique, collapse = "; "),
            "); condition numbers may not be comparable.", call. = FALSE)
  }

  nsub_row <- .diag_row("nSub", fmt_num(num_field("nSub")),
                        run_labels, out_cols)
  nobs_row <- .diag_row("nObs", fmt_num(num_field("nObs")),
                        run_labels, out_cols)

  rt_vals <- vapply(xpdb_list, .compare_runtime, numeric(1))
  rt_row  <- .diag_row("total runtime (sec)", fmt_num(rt_vals),
                       run_labels, out_cols)

  n_param_rows <- nrow(out)

  header_rows <- if (isTRUE(drop_dOFV)) {
    list(ofv_row, met_row, ret_row, cn_row, basis_row,
         nsub_row, nobs_row, rt_row)
  } else {
    list(ofv_row, dofv_row, met_row, ret_row, cn_row, basis_row,
         nsub_row, nobs_row, rt_row)
  }

  final <- dplyr::bind_rows(c(header_rows, list(out), list(rse_block)))
  final <- final[, setdiff(colnames(final), "Section"), drop = FALSE]
  final <- final[, c("Description", run_labels), drop = FALSE]
  final <- tibble::as_tibble(final)

  n_rse <- if (is.null(rse_block)) 0L else nrow(rse_block)
  n_header <- nrow(final) - n_param_rows - n_rse

  attr(final, "n_rse") <- n_rse
  attr(final, "n_header") <- n_header
  attr(final, "nModels") <- length(run_labels)
  attr(final, "run_labels") <- run_labels
  class(final) <- c("prmComparisonNlme", class(final))
  final
}

# First-run-seeded canonical order of (Section, Description), appending each
# later run's previously unseen parameters in that run's own order.
.canonical_param_ord <- function(param_long, run_labels) {
  sections <- character(0)
  descs <- character(0)
  seen <- character(0)
  for (rn in run_labels) {
    sub <- param_long[param_long$run == rn,
                      c("Section", "Description", "param_ord"),
                      drop = FALSE]
    if (!nrow(sub)) next
    sec_rank <- match(sub$Section, c("TH", "OM", "SI"), nomatch = 99L)
    sub <- sub[order(sec_rank, sub$param_ord), , drop = FALSE]
    for (i in seq_len(nrow(sub))) {
      key <- paste(sub$Section[i], sub$Description[i], sep = "\r")
      if (key %in% seen) next
      seen <- c(seen, key)
      sections <- c(sections, sub$Section[i])
      descs <- c(descs, sub$Description[i])
    }
  }
  data.frame(Section = sections, Description = descs,
             param_ord = seq_along(descs), stringsAsFactors = FALSE)
}

# Per-run parameter rows in long form, with Pirana-compatible OMEGA
# transforms and delta-method %RSE. Diagonal SIGMA is kept on the reported
# get_prmNlme() scale.
.prep_params_xpnlme <- function(xpdb, run_label, transform) {
  get_prmNlme(xpdb, level = NULL) %>%
    dplyr::mutate(
      Description = label,
      Section = dplyr::case_when(
        type == "the" ~ "TH",
        type == "ome" ~ "OM",
        type == "sig" ~ "SI",
        TRUE ~ ""
      ),
      run = run_label
    ) %>%
    dplyr::group_by(Section) %>%
    dplyr::mutate(param_ord = dplyr::row_number()) %>%
    dplyr::ungroup() %>%
    dplyr::rowwise() %>%
    dplyr::mutate(
      is_omega_sigma = type %in% c("ome", "sig"),
      is_diagonal = dplyr::coalesce(diagonal, FALSE),
      is_omega_diag = is_omega_sigma && is_diagonal && type == "ome",
      is_sigma_diag = is_omega_sigma && is_diagonal && type == "sig",
      # Transform only diagonal OMEGA values; keep diagonal SIGMA on the
      # reported get_prmNlme() scale (for CEps this is the SD scale, so both
      # value and %RSE remain unchanged).
      var_from_value = dplyr::case_when(is_omega_diag ~ value, TRUE ~ NA_real_),
      se_var_from_value = dplyr::case_when(is_omega_diag ~ se, TRUE ~ NA_real_),
      display_value = dplyr::if_else(
        !is_omega_sigma | !is_diagonal | is_sigma_diag,
        value,
        dplyr::if_else(
          transform == "untransformed",
          value,
          dplyr::if_else(
            transform == "sqrt_om2",
            sqrt(var_from_value),
            sqrt(exp(var_from_value) - 1)
          )
        )
      ),
      display_rse = dplyr::if_else(
        !is_omega_sigma | !is_diagonal | is_sigma_diag,
        rse,
        dplyr::if_else(
          transform == "untransformed",
          rse,
          dplyr::if_else(
            transform == "sqrt_om2",
            dplyr::if_else(!is.na(se_var_from_value) & !is.na(var_from_value) &
                             var_from_value != 0,
                           (se_var_from_value / var_from_value) * 0.5, NA_real_),
            dplyr::if_else(!is.na(se_var_from_value) & !is.na(var_from_value) &
                             (exp(var_from_value) - 1) != 0,
                           (exp(var_from_value) / 2 / (exp(var_from_value) - 1)) *
                             se_var_from_value, NA_real_)
          )
        )
      ),
      display = if (is.na(display_value)) NA_character_ else
        if (is.na(display_rse)) sprintf("%g", display_value) else
          sprintf("%g (%.1f%%)", display_value, display_rse * 100),
      disp_val = if (is.na(display_value)) NA_character_ else
        sprintf("%g", display_value),
      disp_rse = if (is.na(display_rse)) NA_character_ else
        sprintf("%.1f%%", display_rse * 100)
    ) %>%
    dplyr::ungroup() %>%
    dplyr::select(Section, Description, param_ord, run,
                  display, disp_val, disp_rse)
}


# --------------------------------------------------------------------------
# Small internal helpers
# --------------------------------------------------------------------------

# Coerce a logical-ish setting to a scalar logical. Accepts legacy
# "YES"/"Y"/"TRUE"/"ON"/"T"/"1" and "NO"/"N"/"FALSE"/"OFF"/"F"/"0" strings
# (case-insensitive). Unrecognized non-logical values error.
.as_logical_flag <- function(x) {
  if (is.logical(x)) return(isTRUE(x[1]))
  s <- toupper(as.character(x)[1])
  if (s %in% c("YES", "Y", "TRUE", "ON", "T", "1")) return(TRUE)
  if (s %in% c("NO", "N", "FALSE", "OFF", "F", "0")) return(FALSE)
  stop("Expected a logical or YES/NO-style string; got: ",
       encodeString(as.character(x)[1], quote = "\""), call. = FALSE)
}

# First non-NA character value, else NA_character_.
.first_nonNA <- function(v) {
  v <- v[!is.na(v)]
  if (length(v)) v[[1]] else NA_character_
}

# Pivot a long (Section, Description, run, value_str) frame to one wide row
# per (Section, Description), with one column per run label (missing cells
# NA). Base-R replacement for tidyr::pivot_wider() so tidyr is not needed.
.pivot_runs_wide <- function(df, run_labels) {
  keys <- unique(df[, c("Section", "Description")])
  mat <- matrix(NA_character_, nrow = nrow(keys), ncol = length(run_labels),
                dimnames = list(NULL, run_labels))
  key_id <- paste(keys$Section, keys$Description, sep = "\r")
  row_idx <- match(paste(df$Section, df$Description, sep = "\r"), key_id)
  col_idx <- match(as.character(df$run), run_labels)
  for (k in seq_len(nrow(df))) {
    if (!is.na(row_idx[k]) && !is.na(col_idx[k])) {
      mat[row_idx[k], col_idx[k]] <- df$value_str[k]
    }
  }
  cbind(keys, as.data.frame(mat, stringsAsFactors = FALSE, check.names = FALSE))
}

# Build a single-row diagnostic data.frame (Section "", given Description),
# aligned to `out_cols`. `values` is aligned to `run_labels`.
.diag_row <- function(description, values, run_labels, out_cols) {
  lst <- c(list(Section = "", Description = description),
           stats::setNames(as.list(unname(values)), run_labels))
  row <- as.data.frame(lst, stringsAsFactors = FALSE, check.names = FALSE)
  row[, out_cols, drop = FALSE]
}

# Estimation method from xpdb$summary (canonical), not files$method.
.compare_method <- function(xpdb) {
  s <- xpdb$summary
  if (is.null(s)) return(NA_character_)
  row <- s[s$label == "method" & s$problem == 1 & s$subprob == 0, ]
  if (!nrow(row)) return(NA_character_)
  as.character(row$value[[nrow(row)]])
}

# Total runtime (secs) = engine runtime + stderr/covariance runtime from
# xpdb$summary (imported from nlme7engine.log). CPU time, not wall-clock.
.compare_runtime <- function(xpdb) {
  s <- xpdb$summary
  if (is.null(s)) return(NA_real_)
  grab <- function(lab) {
    row <- s[s$label == lab & s$problem == 1 & s$subprob == 0, ]
    if (!nrow(row)) return(NA_real_)
    suppressWarnings(as.numeric(row$value[[nrow(row)]]))
  }
  vals <- c(grab("runtime"), grab("covtime"))
  if (all(is.na(vals))) return(NA_real_)
  sum(vals, na.rm = TRUE)
}


# --------------------------------------------------------------------------
# Optional CSV export
# --------------------------------------------------------------------------

# Write the comparison table to CSV in column (models as columns) or row
# (models as rows, transposed) layout. The CSV always contains all models.
.write_prmComparison <- function(tbl, output_file, format = c("column", "row")) {
  format <- match.arg(format)
  df <- as.data.frame(tbl, stringsAsFactors = FALSE, check.names = FALSE)
  # Quote fields (write.csv default) so labels/values containing commas, quotes
  # or newlines round-trip correctly.
  if (format == "column") {
    utils::write.csv(df, output_file, row.names = FALSE)
  } else {
    tbl_row <- as.data.frame(t(rbind(colnames(df), df)),
                             stringsAsFactors = FALSE)
    colnames(tbl_row) <- tbl_row[1, ]
    colnames(tbl_row)[1] <- "Run"
    tbl_row <- tbl_row[-1, , drop = FALSE]
    utils::write.csv(tbl_row, output_file, row.names = FALSE)
  }
  invisible(output_file)
}


# --------------------------------------------------------------------------
# Optional flextable rendering (flextable + officer are Suggests-only)
# --------------------------------------------------------------------------

#' Render a parameter comparison as a `flextable`
#'
#' @description `as_flextable()` method for the `prmComparisonNlme` object
#'   returned by [compare_prmNlme()]. Draws solid separators between the
#'   diagnostic block, the parameter block, and the optional RSE block, and
#'   reports the number of models in the footer. Requires the `flextable` and
#'   `officer` packages (both `Suggests`).
#'
#' @param x A `prmComparisonNlme` tibble from [compare_prmNlme()].
#' @param format One of `"column"` (models as columns) or `"row"` (models as
#'   rows, transposed). Defaults to the `format` attribute stored by
#'   [compare_prmNlme()], or `"column"` when unset.
#' @param max_show Maximum number of models to render. `NULL` or `""` shows
#'   all models. Any CSV written by [compare_prmNlme()] still contains every
#'   model regardless of `max_show`.
#' @param ... Unused; present for S3 generic compatibility.
#'
#' @return A `flextable` object.
#' @seealso [compare_prmNlme()]
#' @exportS3Method flextable::as_flextable
as_flextable.prmComparisonNlme <- function(x, format = NULL,
                                           max_show = NULL, ...) {
  if (!requireNamespace("flextable", quietly = TRUE) ||
      !requireNamespace("officer", quietly = TRUE)) {
    stop("`as_flextable()` for a prmComparisonNlme needs the 'flextable' and ",
         "'officer' packages. Install them, or use the returned tibble / ",
         "`output_file` CSV instead.", call. = FALSE)
  }
  if (is.null(format)) {
    format <- attr(x, "format")
    if (is.null(format)) format <- "column"
  }
  format <- match.arg(format, c("column", "row"))

  n_rse <- attr(x, "n_rse"); if (is.null(n_rse)) n_rse <- 0L
  n_header <- attr(x, "n_header"); if (is.null(n_header)) n_header <- 0L
  n_models <- attr(x, "nModels")
  if (is.null(n_models)) n_models <- ncol(x) - 1L

  tbl <- as.data.frame(x, stringsAsFactors = FALSE, check.names = FALSE)
  solid_line <- officer::fp_border(color = "black", width = 2)

  run_cols <- setdiff(colnames(tbl), "Description")
  n_show <- .resolve_n_show(max_show, length(run_cols))
  show_runs <- run_cols[seq_len(n_show)]
  if (n_show < length(run_cols)) {
    message(sprintf("Displaying %d of %d models in the table.",
                    n_show, length(run_cols)))
  }

  if (format == "column") {
    tbl_disp <- tbl[, c("Description", show_runs), drop = FALSE]
    ft <- flextable::as_flextable(tbl_disp, show_coltype = FALSE,
                                  max_row = nrow(tbl_disp))
    if (n_header > 0) {
      ft <- flextable::hline(ft, i = n_header, border = solid_line)
    }
    if (n_rse > 0) {
      ft <- flextable::hline(ft, i = nrow(tbl_disp) - n_rse, border = solid_line)
    }
  } else {
    tbl_row <- as.data.frame(t(rbind(colnames(tbl), tbl)),
                             stringsAsFactors = FALSE)
    colnames(tbl_row) <- tbl_row[1, ]
    colnames(tbl_row)[1] <- "Run"
    tbl_row <- tbl_row[-1, , drop = FALSE]
    # Pass the full transposed table so flextable uses the multi-row printer
    # (a single-row data.frame becomes a 2-column key/value layout). Cap
    # displayed models via max_row instead of subsetting first.
    ft <- flextable::as_flextable(tbl_row, show_coltype = FALSE,
                                  max_row = n_show)
    ft <- flextable::delete_part(ft, part = "header")
    if (n_header > 0) {
      ft <- flextable::vline(ft, j = n_header + 1, border = solid_line,
                             part = "body")
    }
    if (n_rse > 0) {
      ft <- flextable::vline(ft, j = ncol(tbl_row) - n_rse, border = solid_line,
                             part = "body")
    }
  }

  ft <- flextable::delete_part(ft, part = "footer")
  ft <- flextable::add_footer_lines(ft, values = sprintf("nModels: %d", n_models))
  ft
}

# Resolve how many models to render: NULL / "" -> all; else a positive
# integer capped at the number available.
.resolve_n_show <- function(max_show, n_available) {
  if (is.null(max_show) || length(max_show) == 0 ||
      (length(max_show) == 1 &&
       (is.na(max_show) || identical(as.character(max_show), "")))) {
    return(n_available)
  }
  n <- suppressWarnings(as.integer(max_show))
  if (is.na(n) || n < 1) {
    stop("`max_show` must be a positive whole number, NULL, or blank.",
         call. = FALSE)
  }
  min(n, n_available)
}