R/eta_shrinkage_subject.R

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

#' @importFrom magrittr %>%

.create_etaSubjectDF <- function(dmp.txt, subject_ids) {
  eta_shrinkage <- dmp.txt$eta_shrinkage
  if (is.null(eta_shrinkage) || is.null(eta_shrinkage$Subject)) {
    return(NULL)
  }

  eta_vals <- dmp.txt$coefficients$random$Subject
  omega <- dmp.txt$omega
  shrink_sub <- eta_shrinkage$Subject

  eta_vals_df <- as.data.frame(eta_vals)
  shrink_sub_df <- as.data.frame(shrink_sub)

  # `subject_ids` is built upstream by `rle(d1$ID)$values` so its
  # length equals the number of contiguous ID blocks the engine saw
  # (one eta row per block). For sorted data each ID appears in a
  # single block; for unsorted data a single ID can occur in more
  # than one block, so `subject_ids` may carry duplicates -- that's
  # the engine's view and it's preserved here.
  if (length(subject_ids) != nrow(eta_vals_df)) {
    warning(
      ".create_etaSubjectDF: subject_ids length (", length(subject_ids),
      ") does not match coefficients$random$Subject rows (",
      nrow(eta_vals_df), "); skipping the subject-level eta table.",
      call. = FALSE
    )
    return(NULL)
  }

  shrink_names <- names(shrink_sub_df)
  bare_names <- sub("_shrinkage$", "", shrink_names)
  names(shrink_sub_df) <- bare_names

  eta_names <- names(eta_vals_df)
  common_etas <- intersect(eta_names, bare_names)
  if (length(common_etas) == 0) {
    warning("No matching ETA names between coefficients$random$Subject and ",
            "eta_shrinkage$Subject after stripping '_shrinkage' suffix.",
            call. = FALSE)
    return(NULL)
  }

  omega_mat <- as.matrix(omega)
  omega_diag <- stats::setNames(diag(omega_mat), colnames(omega_mat))

  rows <- vector("list", length(common_etas))
  for (i in seq_along(common_etas)) {
    eta <- common_etas[i]
    omega_var <- omega_diag[eta]
    tibble::tibble(
      ID = subject_ids,
      Eta = eta,
      ETA_VAL = eta_vals_df[[eta]],
      ETA_SE = sqrt(shrink_sub_df[[eta]] * omega_var)
    ) -> rows[[i]]
  }

  dplyr::bind_rows(rows)
}


.empty_special <- function() {
  tibble::tibble(
    problem = integer(0),
    method = character(0),
    type = character(0),
    data = list(),
    modified = logical(0)
  )
}


.build_etaSubject_special <- function(dmp.txt, subject_ids, problem = 1) {
  eta_sub <- .create_etaSubjectDF(dmp.txt, subject_ids)
  if (is.null(eta_sub)) {
    return(.empty_special())
  }

  tibble::tibble(
    problem = problem,
    method = "nlme",
    type = "eta_subject",
    data = list(eta_sub),
    modified = FALSE
  )
}


#' @title Access subject-level ETA data
#'
#' @description Retrieves the stored subject-by-ETA table from the
#'   \code{xpdb$special} slot with \code{type == "eta_subject"}.
#'
#' @param xpdb An \code{xpose_data} object created by
#'   \code{\link{xposeNlme}} or \code{\link{xposeNlmeModel}}.
#' @param .problem The problem number to extract (default \code{1}).
#'
#' @return A tibble with columns \code{ID}, \code{Eta}, \code{ETA_VAL},
#'   \code{ETA_SE}.
#'
#' @export
get_etaSubjectNlme <- function(xpdb, .problem = 1) {
  stopifnot(inherits(xpdb, "xpose_data"))

  sp <- xpdb$special
  if (is.null(sp) || nrow(sp) == 0) {
    stop("No special datasets found in xpdb. ",
         "Subject-level ETA shrinkage data is unavailable.",
         call. = FALSE)
  }

  hit <- sp$problem == .problem & sp$type == "eta_subject"
  if (!any(hit)) {
    stop("No eta_subject data found for problem ", .problem, ". ",
         "Subject-level ETA shrinkage data is unavailable for this problem.",
         call. = FALSE)
  }

  sp$data[which(hit)[1]][[1]]
}


.get_omegaDiagNlme <- function(xpdb, .problem = 1) {
  prm_row <- xpdb$files[xpdb$files$name == "prmTable" &
                           xpdb$files$problem == .problem, ]
  if (nrow(prm_row) == 0) {
    stop("No prmTable found for problem ", .problem, ".", call. = FALSE)
  }

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


.resolve_eta_thresholds <- function(etas, threshold, eta_threshold) {
  stopifnot(is.numeric(threshold), length(threshold) == 1, threshold >= 0)

  resolved <- stats::setNames(rep(threshold, length(etas)), etas)

  if (!is.null(eta_threshold)) {
    if (!is.numeric(eta_threshold) || is.null(names(eta_threshold))) {
      stop("eta_threshold must be a named numeric vector.", call. = FALSE)
    }
    if (any(is.na(names(eta_threshold)) | names(eta_threshold) == "")) {
      stop("Every element of eta_threshold must be named.", call. = FALSE)
    }
    if (anyDuplicated(names(eta_threshold))) {
      stop("Duplicate names in eta_threshold: ",
           paste(names(eta_threshold)[duplicated(names(eta_threshold))],
                 collapse = ", "),
           call. = FALSE)
    }
    if (any(eta_threshold < 0)) {
      stop("All eta_threshold values must be non-negative.", call. = FALSE)
    }
    bad <- setdiff(names(eta_threshold), etas)
    if (length(bad) > 0) {
      stop("eta_threshold names not found in stored ETAs: ",
           paste(bad, collapse = ", "),
           call. = FALSE)
    }
    resolved[names(eta_threshold)] <- eta_threshold
  }

  resolved
}


.pop_sd <- function(x) sqrt(mean((x - mean(x))^2))


.compute_etaShrinkageNlme <- function(eta_subject,
                                      omega_diag,
                                      threshold,
                                      eta_threshold,
                                      eta_name = NULL) {
  all_etas <- unique(eta_subject$Eta)
  resolved <- .resolve_eta_thresholds(all_etas, threshold, eta_threshold)

  target_etas <- if (!is.null(eta_name)) {
    bad <- setdiff(eta_name, all_etas)
    if (length(bad) > 0) {
      stop("eta_name contains ETAs not in the stored subject-level table: ",
           paste(bad, collapse = ", "),
           call. = FALSE)
    }
    eta_name
  } else {
    all_etas
  }

  rows <- vector("list", length(target_etas))
  for (i in seq_along(target_etas)) {
    eta <- target_etas[i]
    thr <- resolved[eta]
    omega_var <- omega_diag[eta]
    omega_sd <- sqrt(omega_var)
    cutoff <- (1 - thr)^2

    sub <- eta_subject[eta_subject$Eta == eta, ]
    n_total <- nrow(sub)

    if (is.na(omega_sd) || omega_sd <= 0) {
      rows[[i]] <- tibble::tibble(
        Eta = eta,
        n_total = n_total,
        n_kept = 0L,
        n_removed = n_total,
        threshold_applied = unname(thr),
        original_shrinkage = NA_real_,
        filtered_shrinkage = NA_real_,
        omega_sd = unname(omega_sd)
      )
      next
    }

    shrink_var <- (sub$ETA_SE / omega_sd)^2
    finite_shrink <- is.finite(shrink_var)
    keep <- finite_shrink & (shrink_var < cutoff)
    kept <- sub[keep, ]

    n_kept <- sum(keep)
    n_removed <- n_total - n_kept

    if (n_kept > 0) {
      filtered_shrinkage <- 1 - .pop_sd(kept$ETA_VAL) / omega_sd
    } else {
      filtered_shrinkage <- NA_real_
    }

    original_shrinkage <- 1 - .pop_sd(sub$ETA_VAL) / omega_sd

    rows[[i]] <- tibble::tibble(
      Eta = eta,
      n_total = n_total,
      n_kept = n_kept,
      n_removed = n_removed,
      threshold_applied = unname(thr),
      original_shrinkage = unname(original_shrinkage),
      filtered_shrinkage = unname(filtered_shrinkage),
      omega_sd = unname(omega_sd)
    )
  }

  dplyr::bind_rows(rows)
}


.format_etaShrinkage_summary <- function(shrinkage_df,
                                         current_summary = NULL) {
  new_entries <- stats::setNames(shrinkage_df$filtered_shrinkage,
                                 shrinkage_df$Eta)

  if (!is.null(current_summary) && !is.na(current_summary)) {
    parts <- trimws(strsplit(current_summary, ",")[[1]])
    parsed <- stats::setNames(
      vapply(parts, function(p) trimws(strsplit(p, "=")[[1]][2]),
             character(1)),
      vapply(parts, function(p) trimws(strsplit(p, "=")[[1]][1]),
             character(1))
    )
    parsed[names(new_entries)] <- ifelse(
      is.na(new_entries), "NA",
      as.character(round(new_entries, 5))
    )
    all_entries <- parsed
  } else {
    all_entries <- ifelse(
      is.na(new_entries), "NA",
      as.character(round(new_entries, 5))
    )
    names(all_entries) <- names(new_entries)
  }

  sorted <- all_entries[sort(names(all_entries))]
  paste(names(sorted), "=", sorted, collapse = ", ")
}