R/run_single_point_base.R

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

#' Calculate clearance using an adaptive single-point method
#'
#' Calculates clearance using an adaptive single-point pharmacokinetic method
#'
#' @param dat A data frame containing pharmacokinetic data.
#'   Required columns typically include ID, TIME, DV, tad, recent_ii,
#'   dose, routeobs, and durationobs.
#'
#' @param half_life Optional numeric value for the drug's half-life.
#'   If not provided, half-life is estimated using `get_hf()` from pooled observations.
#'
#' @param dose_type Specifies the dosing context of the pharmacokinetic
#'   observations. Required when half_life is not provided. Classified as
#'   first_dose, repeated_doses, or combined_doses based on whether observed
#'   concentrations occur following the first administration, during repeated
#'   dosing, or across both contexts.
#'
#' @param pooled_ctrl Optional list of control parameters used by `get_pooled_data()`
#'   for pooling observations. Defaults to output from `pooled_control()`.
#'
#' @param ssctrl A list of control parameters generated by
#'   `ss_control()` to guide the detection of steady-state observations.
#'
#' @details
#' Estimates individual and population clearance from steady-state
#' pharmacokinetic data. If half-life is not provided, it is estimated
#' from pooled data using `get_hf()` and pooling rules defined in
#' `pooled_control()`.
#'
#' The procedure:
#' \itemize{
#'   \item Identifies steady-state observations using `is_ss()` and
#'         `ss_control()` criteria.
#'   \item Selects peak and trough concentrations within each dose
#'         interval to represent steady-state behavior.
#'   \item Classifies concentration points as Cssmax, Cssmin, or
#'         Cssavg based on timing within the interval and decay pattern.
#'   \item Computes individual clearance as Dose / (Cssavg × tau).
#'   \item Aggregates individual clearance values using a trimmed
#'         geometric mean to obtain a population estimate.
#' }
#'
#' Supports bolus, infusion, and oral administration routes.
#'
#' @return
#' A list containing:
#'   - dat: the processed dataset with steady-state identification
#'   - cl_df: individual clearance estimates
#'   - trimmed_mean_cl: the population clearance calculated as a trimmed
#'     geometric mean with a 5 percent trimming level to reduce the impact
#'     of outliers
#'
#' @author Zhonghui Huang
#'
#' @examples
#' dat <- processData(Bolus_1CPT)$dat
#' calculate_cl(dat, get_hf(dat)$half_life_median)$trimmed_mean_cl
#'
#' @seealso
#' \code{\link{get_hf}}, \code{\link{get_pooled_data}},
#' \code{\link{pooled_control}}, \code{\link{is_ss}}, \code{\link{ss_control}},
#' \code{\link{trimmed_geom_mean}}
#'
#' @export

calculate_cl <- function(dat,
                         half_life = NULL,
                         dose_type = NULL,
                         pooled_ctrl = pooled_control(),
                         ssctrl = ss_control()) {

  # Estimate half-life if not provided
  if (is.null(half_life)) {
    # Defensive check: dose_type is required when half_life is not provided
    if (is.null(half_life) && is.null(dose_type)) {
      stop("Argument 'dose_type' must be provided when 'half_life' is NULL.")
    }

    # Step 1: Generate pooled data using dose type and binning control
    pooled_data <- get_pooled_data(dat,
                                   pooled_ctrl = pooled_ctrl)

    # Step 2: Estimate the half-life from pooled data
    half_life_out <- get_hf(dat = dat,
                            pooled = pooled_data)

    half_life <- half_life_out$half_life_median

    # Defensive check: If half-life estimation failed, stop the process
    if (is.na(half_life)) {
      stop("Half-life estimation failed. Adaptive single-point method cannot proceed.")
    }
  }

  dat <- is_ss(df = dat,
               half_life = half_life,
               ssctrl = ssctrl)

  trimmed_mean_cl <- NA
  dat.ss.obs  <- NA

  dat.ss.obs <- dat %>%
    dplyr::filter(SteadyState == TRUE, !is.na(SteadyState))

  if (nrow(dat.ss.obs) <= 2) {
    return(list(
      cl_df = dat.ss.obs,
      trimmed_mean_cl = trimmed_mean_cl,
      dat = dat
    ))
  }

  # If there are multiple points within the same dose interval,
  # only the min and max values are selected for calculation
  dat.ss.obs <- dat.ss.obs %>%
    dplyr::group_by(ID, dose_number) %>%
    dplyr::mutate(
      max_value = max(DV),
      min_value = min(DV),
      max_time = tad[which.max(DV)],
      min_time = tad[which.min(DV)],
      avg_value = (max_value + min_value) / 2,
      max_interval = ifelse(DV == max(DV), TRUE, FALSE),
      min_interval = ifelse(DV == min(DV), TRUE, FALSE)
    ) %>%
    dplyr::ungroup()

  # Extract rows to be marked in dat.ss.obs
  rows_to_mark <- dat.ss.obs %>%
    dplyr::filter(max_interval == TRUE | min_interval == TRUE) %>%
    dplyr::select(ID, dose_number, TIME)

  dat <- dat %>%
    dplyr::rowwise() %>%
    dplyr::mutate(SteadyState = ifelse(
      any(
        ID == rows_to_mark$ID &
          dose_number == rows_to_mark$dose_number &
          TIME == rows_to_mark$TIME
      ),
      TRUE,
      FALSE
    )) %>%
    dplyr::ungroup()

  # Same for dat.ss.obs
  dat.ss.obs  <- dat.ss.obs  %>%
    dplyr::rowwise() %>%
    dplyr::mutate(SteadyState = ifelse(
      any(
        ID == rows_to_mark$ID &
          dose_number == rows_to_mark$dose_number &
          TIME == rows_to_mark$TIME
      ),
      TRUE,
      FALSE
    )) %>%
    dplyr::ungroup()

  # Second selection， only select the max, min points
  dat.ss.obs <- dat.ss.obs %>%
    dplyr::filter(SteadyState == TRUE, !is.na(SteadyState))

  # avg_value type identification
  # only for fast oral absorption when most absorption happens before Tmax
  dat.ss.obs  <-   dat.ss.obs %>%
    dplyr::mutate(
      # Init 'Css_type' to "Css_avg"
      Css_type = "Css_avg",

      # Determine Css_type using conditions
      Css_type = dplyr::case_when(
        # Condition 1: Css_avg if exp decay small
        # exp(-k * tau) >= 0.6667 means:
        # - Css_max/Css_min < 1.5
        # - Css_avg close to Css_min/max
        exp(-log(2) / half_life * recent_ii) >= 0.6667 ~ "Css_avg",

        # Condition 2: Css_avg if ratio exceeds
        # model threshold from decay curve
        max_value / min_value > exp(log(2) / half_life * recent_ii) ~ "Css_avg",

        # Condition 3: Css_max if both times
        # are in first 20% of interval
        # - Uses recent_ii (not half-life) for better robustness
        (max_time <= 0.2 * recent_ii & max_time != 0) &
          (min_time <= 0.2 * recent_ii & min_time != 0) ~ "Css_max",

        # Condition 4: Css_min if times are
        # in last 20% or if time == 0
        # - Uses recent_ii (not half-life) for better robustness
        (max_time >= 0.8 * recent_ii | max_time == 0) &
          (min_time >= 0.8 * recent_ii | min_time == 0) ~ "Css_min",

        # Default: Css_avg
        TRUE ~ Css_type
      )
    )

  dat.ss.obs <- dat.ss.obs %>%
    dplyr::mutate(
      Css_min_i = dplyr::case_when(
        Css_type == "Css_max" & routeobs == "bolus" ~
          max_value * exp(-log(2) / as.numeric(half_life) * recent_ii),
        Css_type == "Css_max" & routeobs == "infusion" ~
          max_value * exp(-log(2) / as.numeric(half_life) * (recent_ii - durationobs)),
        TRUE ~ NA_real_
      ),
      Css_max_i = dplyr::case_when(
        Css_type == "Css_min" & routeobs == "bolus" ~
          min_value / exp(-log(2) / as.numeric(half_life) * recent_ii),
        Css_type == "Css_min" & routeobs == "infusion" ~
          min_value / exp(-log(2) / as.numeric(half_life) * (recent_ii - durationobs)),
        TRUE ~ NA_real_
      ),
      Css_avg_i = dplyr::case_when(
        Css_type == "Css_avg" ~ avg_value,
        Css_type == "Css_max" ~ (max_value + Css_min_i) / 2,
        Css_type == "Css_min" ~ (min_value + Css_max_i) / 2,
        TRUE ~ NA_real_
      ),
      cl = signif(dose / Css_avg_i / recent_ii, digits = 3)
    )

  # Calculate geometric mean cl for each individual
  individual_mean_cl <-
    tryCatch(
      aggregate(cl ~ ID, data = dat.ss.obs, FUN = trimmed_geom_mean),
      error = function(e) {
        NA
      }
    )

  # Calculate the trimmed mean (e.g., 10% trimmed mean to reduce outlier impact)
  trimmed_mean_cl <-
    tryCatch(
      trimmed_geom_mean(individual_mean_cl$cl, trim = 0.05, na.rm = TRUE),
      error = function(e) {
        NA
      }
    )

  return(list(
    dat = dat,
    cl_df = dat.ss.obs,
    trimmed_mean_cl = trimmed_mean_cl
  ))

}



#' Calculates volume of distribution from concentration data
#'
#' Calculates the volume of distribution (Vd) using an adaptive single-point approach
#'
#' @param dat A data frame containing raw time–concentration data in the
#'   standard nlmixr2 format.
#'
#' @param half_life Optional numeric value for the drug's half-life.
#'   If not provided, it will be estimated using `get_hf()` from pooled
#'   observations.
#'
#' @param dose_type Specifies the dosing context of the pharmacokinetic
#'   observations. Required when half_life is not provided. Classified as
#'   first_dose, repeated_doses, or combined_doses based on whether observed
#'   concentrations occur following the first administration, during repeated
#'   dosing, or across both contexts.
#'
#' @param pooled_ctrl Optional list of control parameters used by `get_pooled_data()`
#'   for pooling observations. Defaults to output from `pooled_control()`.
#'
#' @param route Character string specifying the route of administration. Must be
#'   one of bolus, oral, or infusion. Currently, oral is not implemented.
#'
#' @details
#' The function uses a concentration observed within the first 20% of the elimination
#' half-life after dosing as the early point for estimating the volume of distribution.
#'
#' \deqn{Vd = \frac{\text{Dose}}{C_0}}
#' For infusion:
#' \deqn{Vd = \frac{\text{Rate} \times \min(\text{TIME}, \text{durationobs})}{C_0}}
#'
#' Here, \eqn{C_0} represents the early concentration observed within the first 20%
#' of the elimination half-life after dosing, which is used as an approximation of
#' the initial concentration for estimating volume of distribution (Vd).
#'  `TIME` refers to time after dose; `durationobs` is the actual infusion duration.
#'
#' When half_life is not provided, it is estimated from pooled data using
#' the functions `get_pooled_data()` and `get_hf()`.
#'
#' @return
#' A list with two elements:
#'   - vd_df: individual volume of distribution estimates
#'   - trimmed_mean_vd: population volume of distribution estimated as a
#'     trimmed geometric mean using a 5 percent trimming level

#' @author Zhonghui Huang
#'
#' @seealso \link{get_pooled_data}, \link{get_hf}, \link{trimmed_geom_mean}
#'
#' @examples
#'
#' dat <- Bolus_1CPT
#' out <- processData(dat)
#' fdat<- out$dat
#' froute <-out$Datainfo$Value[out$Datainfo$Infometrics == "Dose Route"]
#' half_life <- get_hf(dat = fdat)$half_life_median
#' calculate_vd(dat = fdat, half_life = half_life,route=froute)$trimmed_mean_vd
#'
#' @export

calculate_vd <- function(dat,
                         half_life = NULL,
                         dose_type = NULL,
                         pooled_ctrl = pooled_control(),
                         route = c("bolus", "oral", "infusion")) {

  # Estimate half-life if not provided
  if (is.null(half_life)) {
    # Defensive check: dose_type is required when half_life is not provided
    if (is.null(half_life) && is.null(dose_type)) {
      stop("Argument 'dose_type' must be provided when 'half_life' is NULL.")
    }

    # Step 1: Generate pooled data using dose type and binning control
    pooled_data <- get_pooled_data(dat,
                                   pooled_ctrl  = pooled_ctrl)

    # Step 2: Estimate the half-life from pooled data
    half_life_out <- get_hf(dat = dat,
                            pooled = pooled_data)

    half_life <- half_life_out$half_life_median

    # Defensive check: If half-life estimation failed, stop the process
    if (is.na(half_life)) {
      stop("Half-life estimation failed. Adaptive single-point method cannot proceed.")
    }
  }

  # ---- Defensive check: ensure route is valid ----
  route <- tryCatch(
    match.arg(route, choices = c("bolus", "oral", "infusion")),
    error = function(e) {
      stop(sprintf(
        "Invalid `route`: '%s'. Must be one of: %s.",
        as.character(route),
        paste(shQuote(c(
          "bolus", "oral", "infusion"
        )), collapse = ", ")
      ),
      call. = FALSE)
    }
  )

  trimmed_mean_vd <- NA
  dat.fd.obs <- NA

  dat <- dat %>%
    dplyr::mutate(
      C_first_flag = ifelse(
        EVID == 0 &
          dose_number == 1 &
          tad < half_life * 0.2 &
          iiobs == 0,
        1,
        0
      )
    ) %>%
    dplyr::group_by(ID) %>%
    dplyr::mutate(
      # Check if any rows meet the initial flagging condition
      has_any_C_first = any(C_first_flag == 1),
      # Compute min_time only if there are flagged rows, else NA
      min_time = ifelse(
        has_any_C_first,
        min(TIME[C_first_flag == 1], na.rm = TRUE),
        NA_real_
      ),
      # Update C_first_flag to 1 only for the earliest TIME in flagged rows
      C_first_flag = ifelse(
        C_first_flag == 1 & TIME == min_time,
        1,
        0
      )
    ) %>%
    dplyr::select(-has_any_C_first, -min_time) %>%
    dplyr::ungroup()

  dat.fd.obs <- dat %>%
    dplyr::filter(C_first_flag == 1)

  if (route == "bolus") {
    dat.fd.obs <- dat.fd.obs %>%
      dplyr::mutate(vd = signif(dose / DV, 3))

    trimmed_mean_vd <-
      trimmed_geom_mean(dat.fd.obs$vd, trim = 0.05, na.rm = TRUE)

  } else if (route == "infusion") {
    dat.fd.obs <- dat.fd.obs %>%
      dplyr::mutate(vd = signif(pmin(TIME, durationobs) * rateobs / DV, 3))

    individual_mean_vd <- tryCatch(
      dat.fd.obs %>%
        dplyr::group_by(ID) %>%
        dplyr::summarise(
          vd = trimmed_geom_mean(vd, trim = 0.05, na.rm = TRUE),
          .groups = "drop"
        ),
      error = function(e)
        NA
    )

    trimmed_mean_vd <- tryCatch(
      trimmed_geom_mean(individual_mean_vd$vd, trim = 0.05, na.rm = TRUE),
      error = function(e)
        NA
    )
  }

  return(list(vd_df = dat.fd.obs,
              trimmed_mean_vd = trimmed_mean_vd))
}



#' Run adaptive single-point pharmacokinetic analysis
#'
#' Implements adaptive single-point pharmacokinetic analysis to calculate
#' clearance and volume of distribution.
#'
#' @param dat A data frame containing raw time–concentration data in the
#'   standard nlmixr2 format.
#' @param route Route of administration. Must be one of bolus, oral, or infusion.
#' @param half_life Optional numeric value for drug half-life. If not provided,
#'   it is estimated from the dataset.
#' @param dose_type Specifies the dosing context of the pharmacokinetic
#'   observations. Required when half_life is not provided. Classified as
#'   first_dose, repeated_doses, or combined_doses based on whether observed
#'   concentrations occur following the first administration, during repeated
#'   dosing, or across both contexts.
#' @param pooled_ctrl Optional list of control parameters used by `get_pooled_data()`
#'   for pooling observations. Defaults to output from `pooled_control()`.
#' @param ssctrl A list of control parameters generated by
#'   `ss_control()` to guide the detection of steady-state observations.
#'
#' @return
#' A list containing:
#'   - summary: a data frame with trimmed mean clearance and volume of
#'     distribution, and run time information
#'   - dat: the processed dataset used for analysis
#'   - cl_df: individual clearance estimates
#'   - vd_df: individual volume of distribution estimates
#'
#' @details
#' This function integrates clearance and volume estimation into a unified
#' adaptive workflow, using steady-state pharmacokinetic observations and
#' trimmed mean statistics to reduce the influence of outliers.
#'
#' @seealso \link{calculate_cl}, \link{calculate_vd},
#' \link{pooled_control}, \link{ss_control}
#'
#' @examples
#' dat <- Bolus_1CPT
#' out <- processData(dat)
#' fdat <- out$dat
#' route <- out$Datainfo$Value[out$Datainfo$Infometrics == "Dose Route"]
#' half_life <- get_hf(dat = fdat)$half_life_median
#' run_single_point_base(dat = fdat, half_life = half_life, route = route)$summary
#'
#' @export

run_single_point_base <- function(dat,
                             route = c("bolus", "oral", "infusion"),
                             half_life = NULL,
                             dose_type = NULL,
                             pooled_ctrl = pooled_control(),
                             ssctrl = ss_control()) {

  start.time <- Sys.time()

  # Defensive check for route
  route <- tryCatch(
    match.arg(route, choices = c("bolus", "oral", "infusion")),
    error = function(e) {
      stop(sprintf(
        "Invalid `route`: '%s'. Must be one of: %s.",
        as.character(route),
        paste(shQuote(c(
          "bolus", "oral", "infusion"
        )), collapse = ", ")
      ),
      call. = FALSE)
    }
  )

  # Calculate CL
  cl_out <- calculate_cl(
    dat = dat,
    half_life = half_life,
    dose_type = dose_type,
    pooled_ctrl = pooled_ctrl,
    ssctrl = ssctrl
  )
  dat <- cl_out$dat

  # Calculate Vd
  vd_out <- calculate_vd(
    dat = dat,
    half_life = half_life,
    dose_type = dose_type,
    pooled_ctrl = pooled_ctrl,
    route = route
  )

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

  #  Return summary and details
  summary_df <- data.frame(
    cl = signif(cl_out$trimmed_mean_cl, 3),
    vd = signif(vd_out$trimmed_mean_vd, 3),
    start.time =  start.time,
    time.spent = time.spent
  )

  return(list(
    summary = summary_df,
    dat = dat,
    cl_df = cl_out$cl_df,
    vd_df = vd_out$vd_df
  ))
}

#' Computes the trimmed geometric mean
#'
#' Computes the trimmed geometric mean of a numeric vector
#'
#' @param x A numeric vector containing the values for geometric mean
#'   calculation.
#' @param trim A numeric value between 0 and 0.5 indicating the proportion of
#'   values to be trimmed from each end of the vector. Default is 0.
#' @param na.rm Logical value indicating whether missing values should be
#'   removed before computation. Default is TRUE.
#'
#' @return
#' A numeric value representing the trimmed geometric mean. Returns NA if no
#' values remain after trimming.
#'
#' @examples
#' x <- c(1, 2, 3, 4, 5, 100)
#' trimmed_geom_mean(x, trim = 0.05)
#'
#' @export
#'
trimmed_geom_mean <- function(x, trim = 0, na.rm = TRUE) {
  if (na.rm)
    x <- stats::na.omit(x)  # Remove NA values if na.rm is TRUE

  # Sort the vector and apply trimming
  x <- sort(x)
  n <- length(x)
  lower <- floor(n * trim) + 1  # Index of first value to keep
  upper <- n - floor(n * trim)  # Index of last value to keep

  # Keep only trimmed range
  x_trimmed <- x[lower:upper]

  # Calculate the geometric mean of the trimmed values
  exp(mean(log(x_trimmed), na.rm = FALSE))
}