Nothing
R/sim_sens_params.R
In nlmixr2autoinit: Automatic Generation of Initial Estimates for Population Pharmacokinetic Modeling
Defines functions
sim_sens_1cmpt_mm
Documented in
sim_sens_1cmpt_mm
#' Parameter sweeping for a one-compartment Michaelis-Menten model
#'
#' Performs parameter sweeping by varying pharmacokinetic parameters in a
#' one-compartment model with Michaelis-Menten elimination.
#'
#' @param dat A data frame containing raw time–concentration data in the
#' standard nlmixr2 format.
#' @param sim_vmax List specifying Vmax:
#' - mode: "manual" or "auto"
#' - values: numeric vector if mode = "manual"
#' - est.cl: required if mode = "auto"
#' @param sim_km List specifying Km:
#' - mode: "manual" or "auto"
#' - values: numeric vector if mode = "manual"
#' @param sim_vd List specifying Vd:
#' - mode: must be "manual"
#' - values: numeric vector
#' @param sim_ka List specifying Ka (oral route only):
#' - mode: must be "manual"
#' - values: numeric vector
#' @param route Dosing route, either "iv" or "oral". Default is "iv".
#'
#' @param verbose Logical (default = TRUE).
#' Controls whether progress information is displayed during parameter sweeping.
#' When TRUE, a dynamic progress bar is shown using the `progressr` package to
#' indicate simulation status and elapsed time.
#' When FALSE, progress output is suppressed and the function runs silently.
#'
#' @details
#' The function generates a parameter grid and performs model fitting for each
#' combination using `Fit_1cmpt_mm_iv`. Parameters can be specified manually or
#' automatically derived. Model predictions and fit metrics are computed for each
#' simulation to assess parameter sensitivity.
#'
#' @return A data frame containing parameter combinations and model fit metrics.
#'
#' @author Zhonghui Huang
#'
#' @seealso \code{\link{Fit_1cmpt_mm_iv}}, \code{\link{Fit_1cmpt_mm_oral}}
#' @examples
#' \donttest{
#' # Example: IV dosing scenario with automatic Vmax and Km
#' out <- sim_sens_1cmpt_mm(
#' dat = Bolus_1CPTMM[Bolus_1CPTMM$ID<50,],
#' sim_vmax = list(mode = "auto", est.cl = 4),
#' sim_km = list(mode = "auto"),
#' sim_vd = list(mode = "manual", values = 70),
#' sim_ka = list(mode = "manual", values = NA),
#' route = "iv"
#' )
#' head(out[out$rRMSE2==min(out$rRMSE2),])
#'}
#' @export
sim_sens_1cmpt_mm <- function(dat,
sim_vmax = list(mode = "auto",
values = NULL,
est.cl = NULL),
sim_km = list(mode = "auto",
values = NULL),
sim_vd = list(mode = "manual", values = NULL),
sim_ka = list(mode = "manual", values = NULL),
route = c("iv", "oral"),
verbose= TRUE) {
`%>%` <- magrittr::`%>%`
route <- tryCatch(
match.arg(route, choices = c("iv", "oral")),
error = function(e) {
stop(
sprintf(
"Invalid value for `route`: '%s'. Must be 'iv' or 'oral'.",
as.character(route)
),
call. = FALSE
)
}
)
sim_vmax$mode <- tryCatch(
match.arg(sim_vmax$mode, choices = c("manual", "auto")),
error = function(e) {
stop("Invalid `sim_vmax$mode`. Must be 'manual' or 'auto'.",
call. = FALSE)
}
)
sim_km$mode <- tryCatch(
match.arg(sim_km$mode, choices = c("manual", "auto")),
error = function(e) {
stop("Invalid `sim_km$mode`. Must be 'manual' or 'auto'.", call. = FALSE)
}
)
sim_vd$mode <- tryCatch(
match.arg(sim_vd$mode, choices = c("manual")),
error = function(e) {
stop("Invalid `sim_vd$mode`. Must be 'manual'.", call. = FALSE)
}
)
vd_values <- sim_vd$values
if (is.null(vd_values) || all(is.na(vd_values))) {
stop("No candidate Vd values available for parameter sweeping.",
call. = FALSE)
}
# Deduplicate similar values if values are close
vd_values <- stats::na.omit(vd_values) %>%
sort() %>%
tibble::tibble(value = .) %>%
dplyr::mutate(prev = dplyr::lag(value),
rel_diff = abs(value - prev) / prev) %>%
dplyr::filter(is.na(rel_diff) | rel_diff > 0.2) %>%
dplyr::pull(value)
if (route == "oral") {
sim_ka$mode <- tryCatch(
match.arg(sim_ka$mode, choices = c("manual")),
error = function(e) {
stop("Invalid `sim_ka$mode`. Must be 'manual'.", call. = FALSE)
}
)
if (is.null(sim_ka$values) || all(is.na(sim_ka$values))) {
stop("No candidate Ka values available for parameter sweeping (oral route).",
call. = FALSE)
}
}
# Handle Vmax/Km grid
if (sim_vmax$mode == "auto" || sim_km$mode == "auto") {
if (is.null(sim_vmax$est.cl)) {
stop("Auto mode for Vmax/Km requires `sim_vmax$est.cl`.",
call. = FALSE)
}
cl_values <- sim_vmax$est.cl
# Deduplicate similar values if valus are close
cl_values <- stats::na.omit(cl_values) %>%
sort() %>%
tibble::tibble(value = .) %>%
dplyr::mutate(prev = dplyr::lag(value),
rel_diff = abs(value - prev) / prev) %>%
dplyr::filter(is.na(rel_diff) | rel_diff > 0.2) %>%
dplyr::pull(value)
# obtain the observed cmax
dat.obs <- dat[dat$EVID == 0, ]
pop.cmax <- aggregate(dat.obs$DV,
list(dat.obs$ID),
FUN = max,
na.rm = TRUE)
cmax <- mean(pop.cmax$x, trim = 0.05, na.rm = TRUE)
km_range <- c(4, 2, 1, 0.5, 0.25, 0.125, 0.1, 0.05) * cmax
conc_range <- c(0.05, 0.1, 0.25, 0.5, 0.75) * cmax
combs <-
expand.grid(Km = km_range, Conc = conc_range, cl = cl_values)
combs$Vmax <- (combs$Km + combs$Conc) * combs$cl
param_grid <- combs %>% dplyr::select(Vmax, Km)
# Filter combinations similar to each other
keep <- rep(TRUE, nrow(param_grid))
for (i in 1:(nrow(param_grid) - 1)) {
if (!keep[i])
next
for (j in (i + 1):nrow(param_grid)) {
if (!keep[j])
next
vmax_diff <-
abs((param_grid$Vmax[i] - param_grid$Vmax[j]) / param_grid$Vmax[j])
km_diff <-
abs((param_grid$Km[i] - param_grid$Km[j]) / param_grid$Km[j])
if (vmax_diff <= 0.2 && km_diff <= 0.2) {
keep[j] <- FALSE
}
}
}
param_grid <- param_grid[keep, ]
} else {
if (is.null(sim_vmax$values) || all(is.na(sim_vmax$values))) {
stop("No candidate Vmax values available for parameter sweeping.",
call. = FALSE)
}
if (is.null(sim_km$values) || all(is.na(sim_km$values))) {
stop("No candidate Km values available for parameter sweeping.",
call. = FALSE)
}
param_grid <-
tidyr::crossing(Vmax = sim_vmax$values, Km = sim_km$values)
}
if (route == "oral") {
param_grid <-
tidyr::crossing(param_grid, Vd = vd_values, Ka = sim_ka$values)
} else {
param_grid <- tidyr::crossing(param_grid, Vd = vd_values) %>%
dplyr::mutate(Ka = NA_real_)
}
start_time <- Sys.time()
if (isTRUE(verbose)) {
handler.lists <- list(
progressr::handler_txtprogressbar(
format = ":message [:bar] :percent (:current/:total)",
width = 80
)
)
} else {
handler.lists <- list(progressr::handler_void())
}
sim.1cmpt.mm.results.all <- progressr::with_progress({
p <- progressr::progressor(steps = nrow(param_grid))
param_grid %>%
dplyr::mutate(row = dplyr::row_number()) %>%
purrr::pmap_dfr(function(Vmax, Km, Vd, Ka, row) {
p(sprintf("Running simulation: Vmax=%.2f, Km=%.2f", Vmax, Km))
sim_out <- if (route == "iv") {
suppressMessages(suppressWarnings(Fit_1cmpt_mm_iv(
data = dat[dat$EVID != 2,],
est.method = "rxSolve",
input.vmax = Vmax,
input.km = Km,
input.vd = Vd,
input.add = 0
)))
} else {
suppressMessages(suppressWarnings(Fit_1cmpt_mm_oral(
data = dat[dat$EVID != 2,],
est.method = "rxSolve",
input.ka = Ka,
input.vmax = Vmax,
input.km = Km,
input.vd = Vd,
input.add = 0
)))
}
met <-
metrics.(pred.x = sim_out$cp, obs.y = dat[dat$EVID == 0,]$DV)
elapsed <-
round(difftime(Sys.time(), start_time, units = "secs"), 2)
# rm(sim_out)
# gc()
tibble::tibble(
Vmax = Vmax,
Km = Km,
Vd = Vd,
Ka = if (route == "oral")
Ka
else
NA_real_,
APE = round(met[1], 2),
MAE = round(met[2], 2),
MAPE = round(met[3], 2),
RMSE = round(met[4], 2),
rRMSE1 = round(met[5], 2),
rRMSE2 = round(met[6], 2),
Cumulative.Time.Sec = as.numeric(elapsed)
)
})
},handlers = handler.lists)
return(sim.1cmpt.mm.results.all)
}
#' Parameter sweeping for a two-compartment pharmacokinetic model
#'
#' Performs parameter sweeping by varying pharmacokinetic parameters in a
#' two-compartment model under IV or oral dosing. Model fit is evaluated across
#' combinations of CL, Vc, Vp, Q, and Ka (oral only).
#'
#' @details
#' The function generates a parameter grid and performs model fitting for each
#' combination using `Fit_2cmpt_iv` or `Fit_2cmpt_oral`. Parameters can be
#' specified manually or automatically derived. Model predictions and fit metrics
#' are computed for each simulation to assess parameter sensitivity.
#'
#' @param dat Pharmacokinetic dataset.
#' @param sim_ka List specifying Ka (oral route only):
#' - mode: must be "manual"
#' - values: numeric vector
#' @param sim_cl List specifying clearance (CL):
#' - mode: must be "manual"
#' - values: numeric vector
#' @param sim_vc List specifying central volume (Vc):
#' - mode: must be "manual"
#' - values: numeric vector
#' @param sim_vp List specifying peripheral volume (Vp):
#' - mode: "manual" or "auto"
#' - values: numeric vector if manual
#' @param sim_q List specifying inter-compartmental clearance (Q):
#' - mode: "manual" or "auto"
#' - values: numeric vector if manual
#' @param route Dosing route, either "iv" or "oral". Default is "iv".
#'
#' @param verbose Logical (default = TRUE).
#' Controls whether progress information is displayed during parameter sweeping.
#' When TRUE, a dynamic progress bar is shown using the `progressr` package to
#' indicate simulation status and elapsed time.
#' When FALSE, progress output is suppressed and the function runs silently.
#'
#' @return A data frame containing parameter combinations with model fit metrics.
#'
#' @author Zhonghui Huang
#'
#' @seealso \code{\link{Fit_2cmpt_iv}}, \code{\link{Fit_2cmpt_oral}}
#' @examples
#' \donttest{
#' out <- sim_sens_2cmpt(
#' dat = Bolus_2CPT[Bolus_2CPT$ID<50,],
#' sim_cl = list(mode = "manual", values = 4),
#' sim_vc = list(mode = "manual", values = 50),
#' sim_vp = list(mode = "auto"),
#' sim_q = list(mode = "auto"),
#' sim_ka = list(mode = "manual", values = NA),
#' route = "iv",verbose=FALSE
#' )
#' head(out[out$rRMSE2==min(out$rRMSE2),])
#' }
#' @export
sim_sens_2cmpt <- function(dat,
sim_ka = list(mode = "manual", values = NULL),
sim_cl = list(mode = "manual", values = NULL),
sim_vc = list(mode = "manual", values = NULL),
sim_vp = list(mode = c("auto", "manual"), values = NULL),
sim_q = list(
mode = c("auto", "manual"),
values = NULL,
auto.strategy = c("scaled", "fixed")
),
route = c("iv", "oral"),
verbose=TRUE) {
`%>%` <- magrittr::`%>%`
# Safe check
# Defensive validation for `route`
route <- tryCatch(
match.arg(route, choices = c("iv", "oral")),
error = function(e) {
stop(
sprintf(
"Invalid value for `route`: '%s'. Must be one of: %s.",
as.character(route),
paste(shQuote(c("iv", "oral")), collapse = ", ")
),
call. = FALSE
)
}
)
# Safe check and set sim ka
sim_ka$mode <- tryCatch(
match.arg(sim_ka$mode, choices = c("manual")),
error = function(e) {
stop("Invalid `sim_ka$mode`. Must be 'manual'.", call. = FALSE)
}
)
if (route == "oral" &&
(is.null(sim_ka$values) ||
length(sim_ka$values) == 0 || all(is.na(sim_ka$values)))) {
stop("No candidate Ka values available for parameter sweeping.",
call. = FALSE)
}
# Defensive checks for sim_cl
sim_cl$mode <- tryCatch(
match.arg(sim_cl$mode, choices = c("manual")),
error = function(e) {
stop("Invalid `sim_cl$mode`. Must be 'manual'.", call. = FALSE)
}
)
if (is.null(sim_cl$values) ||
length(sim_cl$values) == 0 || all(is.na(sim_cl$values))) {
stop("No candidate CL values available for parameter sweeping.",
call. = FALSE)
}
# Safe check and set sim vc
sim_vc$mode <- tryCatch(
match.arg(sim_vc$mode, choices = c("manual")),
error = function(e) {
stop("Invalid `sim_vc$mode`. Must be 'manual'.", call. = FALSE)
}
)
vc_values <- sim_vc$values
if (is.null(vc_values) ||
length(vc_values) == 0 || all(is.na(vc_values))) {
stop("No candidate Vc values available for parameter sweeping.",
call. = FALSE)
}
# Deduplicate similar values if values are close
vc_values <- stats::na.omit(vc_values)%>%
sort() %>%
tibble::tibble(value = .) %>%
dplyr::mutate(prev = dplyr::lag(value),
rel_diff = abs(value - prev) / prev) %>%
dplyr::filter(is.na(rel_diff) | rel_diff > 0.2) %>%
dplyr::pull(value)
# Safe check and set sim vp
sim_vp$mode <- tryCatch(
match.arg(sim_vp$mode, choices = c("manual", "auto")),
error = function(e) {
stop("Invalid `sim_vp$mode`. Must be 'manual' or 'auto'.", call. = FALSE)
}
)
if (sim_vp$mode == "auto") {
vp_ratios <- c(10, 5, 2, 1, 0.5, 0.2, 0.1)
vp_combs <- expand.grid(vc = vc_values, ratio = vp_ratios)
vp_combs$vp <- vp_combs$vc / vp_combs$ratio
} else {
if (is.null(sim_vp$values) ||
length(sim_vp$values) == 0 || all(is.na(sim_vp$values))) {
stop("No candidate Vp values available for parameter sweeping.",
call. = FALSE)
}
vp_combs <- expand.grid(vc = vc_values, vp = sim_vp$values)
}
# Safe check and set sim q
sim_q$mode <- tryCatch(
match.arg(sim_q$mode, choices = c("manual", "auto")),
error = function(e) {
stop("Invalid `sim_q$mode`. Must be 'manual' or 'auto'.", call. = FALSE)
}
)
if (sim_q$mode == "auto") {
strategy <- tryCatch(
match.arg(sim_q$auto.strategy, choices = c("scaled", "fixed")),
error = function(e) {
stop("Invalid `sim_q$auto.strategy`. Must be 'scaled' or 'fixed'.",
call. = FALSE)
}
)
cl_value <- sim_cl$values[1]
if (strategy == "scaled") {
q_values <- c(cl_value / 4, cl_value / 2, cl_value, cl_value * 2)
} else {
q_values <- c(1, 10, 100)
}
} else {
if (is.null(sim_q$values) ||
length(sim_q$values) == 0 || all(is.na(sim_q$values))) {
stop("No candidate Q values available for parameter sweeping.",
call. = FALSE)
}
q_values <- sim_q$values
}
vp_combs <- vp_combs %>%
dplyr::rename(Vc = vc, Vp = vp)
if (route == "oral") {
param_grid <- tidyr::crossing(vp_combs,
Q = q_values,
CL = sim_cl$values,
Ka = sim_ka$values)
} else {
param_grid <- tidyr::crossing(vp_combs,
Q = q_values,
CL = sim_cl$values) %>%
dplyr::mutate(Ka = NA_real_)
}
param_grid <- param_grid %>%
dplyr::select(Vc, Vp, Q, CL, Ka)
# k21 boundary was set (0.01,5)
param_grid <- param_grid%>%
dplyr::filter(
(Vp < 1000 & Q / Vp > 0.01 & Q / Vp <= 5) |
(Vp >= 1000 & Q / Vp > 0.001 & Q / Vp <= 5)
)
# === Begin Simulation ===
start_time <- Sys.time()
if (isTRUE(verbose)) {
handler.lists <- list(
progressr::handler_progress(
format = ":message [:bar] :percent (:current/:total)",
width = 80
)
)
} else {
handler.lists <- list(progressr::handler_void())
}
sim.2cmpt.results.all <- progressr::with_progress({
p <- progressr::progressor(steps = nrow(param_grid))
param_grid %>%
dplyr::mutate(row = dplyr::row_number()) %>%
purrr::pmap_dfr(function(Vc, Vp, Q, CL, Ka, row) {
p(sprintf("Running simulation: Vc=%.2f, Vp=%.2f", Vc, Vp))
sim_out <- if (route == "iv") {
suppressMessages(suppressWarnings(Fit_2cmpt_iv(
data = dat[dat$EVID != 2,],
est.method = "rxSolve",
input.cl = CL,
input.vc2cmpt = Vc,
input.vp2cmpt = Vp,
input.q2cmpt = Q,
input.add = 0
)))
} else {
suppressMessages(suppressWarnings(Fit_2cmpt_oral(
data = dat[dat$EVID != 2,],
est.method = "rxSolve",
input.ka = Ka,
input.cl = CL,
input.vc2cmpt = Vc,
input.vp2cmpt = Vp,
input.q2cmpt = Q,
input.add = 0
)))
}
met <-
metrics.(pred.x = sim_out$cp, obs.y = dat[dat$EVID == 0,]$DV)
elapsed <-
round(difftime(Sys.time(), start_time, units = "secs"), 2)
# sim_out<-NULL
tibble::tibble(
Vc = Vc,
Vp = Vp,
Q = Q,
CL = CL,
Ka = if (route == "oral")
Ka
else
NA,
APE = round(met[1], 2),
MAE = round(met[2], 2),
MAPE = round(met[3], 2),
RMSE = round(met[4], 2),
rRMSE1 = round(met[5], 2),
rRMSE2 = round(met[6], 2),
Cumulative.Time.Sec = as.numeric(elapsed)
)
})
},handlers = handler.lists)
return(sim.2cmpt.results.all)
}
#' Parameter sweeping for a three-compartment pharmacokinetic model
#'
#' Performs parameter sweeping by varying pharmacokinetic parameters in a
#' three-compartment model under IV or oral dosing. Parameter combinations
#' include Vc, Vp1, Vp2, Q1, Q2, CL, and Ka (oral only).
#'
#' @details
#' The function generates a parameter grid and evaluates each combination
#' using `Fit_3cmpt_iv` or `Fit_3cmpt_oral`. Model predictions and fit metrics
#' are calculated for each simulation to assess parameter influence and
#' identify optimal regions of the parameter space. Parameters can be provided
#' manually or derived automatically.
#'
#' @param dat Pharmacokinetic dataset.
#' @param sim_vc List specifying Vc:
#' - mode: must be "manual"
#' - values: numeric vector
#' @param sim_vp List specifying Vp1:
#' - mode: "manual" or "auto"
#' - values: numeric vector if manual
#' @param sim_vp2 List specifying Vp2:
#' - mode: "manual" or "auto"
#' - values: numeric vector if manual
#' @param sim_q List specifying Q1:
#' - mode: "manual" or "auto"
#' - values: numeric vector if manual
#' - auto.strategy: "scaled" or "fixed" when auto
#' @param sim_q2 List specifying Q2:
#' - mode: "manual" or "auto"
#' - values: numeric vector if manual
#' - auto.strategy: "scaled" or "fixed" when auto
#' @param sim_cl List specifying CL:
#' - mode: must be "manual"
#' - values: numeric vector
#' @param sim_ka List specifying Ka (oral route only):
#' - mode: must be "manual"
#' - values: numeric vector
#' @param route Dosing route, either "iv" or "oral". Default is "iv".
#'
#' @param verbose Logical (default = TRUE).
#' Controls whether progress information is displayed during parameter sweeping.
#' When TRUE, a dynamic progress bar is shown using the `progressr` package to
#' indicate simulation status and elapsed time.
#' When FALSE, progress output is suppressed and the function runs silently.
#'
#' @return A data frame containing parameter combinations with model fit metrics.
#'
#' @author Zhonghui Huang
#'
#' @seealso \code{\link{Fit_3cmpt_iv}}, \code{\link{Fit_3cmpt_oral}}
#' @examples
#' \donttest{
#' out <- sim_sens_3cmpt(
#' dat = Bolus_2CPT,
#' sim_cl = list(mode = "manual", values = 4),
#' sim_vc = list(mode = "manual", values = 50),
#' sim_vp = list(mode = "auto"),
#' sim_vp2 = list(mode = "auto"),
#' sim_q = list(mode = "auto", auto.strategy = "scaled"),
#' sim_q2 = list(mode = "auto", auto.strategy = "scaled"),
#' route = "iv",verbose=FALSE
#' )
#' head(out[out$rRMSE2==min(out$rRMSE2),])
#' }
#'
#' @export
sim_sens_3cmpt <- function(dat,
sim_vc = list(mode = "manual", values = NULL),
sim_vp = list(mode = c("auto", "manual"), values = NULL),
sim_vp2 = list(mode = c("auto", "manual"), values = NULL),
sim_q = list(
mode = c("auto", "manual"),
values = NULL,
auto.strategy = c("scaled", "fixed")
),
sim_q2 = list(
mode = c("auto", "manual"),
values = NULL,
auto.strategy = c("scaled", "fixed")
),
sim_cl = list(mode = "manual", values = NULL),
sim_ka = list(mode = "manual", values = NULL),
route = c("iv", "oral"),
verbose=TRUE) {
# --- Route check ---
route <- tryCatch(
match.arg(route, choices = c("iv", "oral")),
error = function(e) {
stop(sprintf(
"Invalid `route`: '%s'. Must be one of: %s.",
as.character(route),
paste(shQuote(c("iv", "oral")), collapse = ", ")
),
call. = FALSE)
}
)
# --- sim_ka ---
sim_ka$mode <- tryCatch(
match.arg(sim_ka$mode, choices = c("manual")),
error = function(e) {
stop("Invalid `sim_ka$mode`. Must be 'manual'.", call. = FALSE)
}
)
if (route == "oral" &&
(is.null(sim_ka$values) ||
length(sim_ka$values) == 0 || all(is.na(sim_ka$values)))) {
stop("No candidate Ka values available for parameter sweeping.",
call. = FALSE)
}
# --- sim_cl ---
sim_cl$mode <- tryCatch(
match.arg(sim_cl$mode, choices = c("manual")),
error = function(e) {
stop("Invalid `sim_cl$mode`. Must be 'manual'.", call. = FALSE)
}
)
if (is.null(sim_cl$values) ||
length(sim_cl$values) == 0 || all(is.na(sim_cl$values))) {
stop("No candidate CL values available for parameter sweeping.",
call. = FALSE)
}
# --- sim_vc ---
sim_vc$mode <- tryCatch(
match.arg(sim_vc$mode, choices = c("manual")),
error = function(e) {
stop("Invalid `sim_vc$mode`. Must be 'manual'.", call. = FALSE)
}
)
vc_values <- sim_vc$values
if (is.null(vc_values) ||
length(vc_values) == 0 || all(is.na(vc_values))) {
stop("No candidate Vc values available for parameter sweeping.",
call. = FALSE)
}
vc_values <- stats::na.omit(vc_values) %>%
sort() %>%
tibble::tibble(value = .) %>%
dplyr::mutate(prev = dplyr::lag(value),
rel_diff = abs(value - prev) / prev) %>%
dplyr::filter(is.na(rel_diff) | rel_diff > 0.2) %>%
dplyr::pull(value)
# --- sim_vp & sim_vp2 ---
sim_vp$mode <- tryCatch(
match.arg(sim_vp$mode, choices = c("auto", "manual")),
error = function(e) {
stop("Invalid `sim_vp$mode`. Must be 'auto' or 'manual'.", call. = FALSE)
}
)
sim_vp2$mode <- tryCatch(
match.arg(sim_vp2$mode, choices = c("auto", "manual")),
error = function(e) {
stop("Invalid `sim_vp2$mode`. Must be 'auto' or 'manual'.", call. = FALSE)
}
)
if (sim_vp$mode == "manual" &&
(is.null(sim_vp$values) || all(is.na(sim_vp$values)))) {
stop("No candidate Vp1 values in manual mode.", call. = FALSE)
}
if (sim_vp2$mode == "manual" &&
(is.null(sim_vp2$values) || all(is.na(sim_vp2$values)))) {
stop("No candidate Vp2 values in manual mode.", call. = FALSE)
}
if (sim_vp$mode == "auto" && sim_vp2$mode == "auto") {
vc_vp_ratio_range <- c(10, 5, 2, 1, 0.5, 0.2, 0.1)
combs <-
expand.grid(estvc = vc_values, vc_vp_ratio = vc_vp_ratio_range)
combs$estvp <- signif(combs$estvc / combs$vc_vp_ratio)
combs2 <- combs
combs$estvp2 <- combs$estvc
combs2$estvp2 <- combs2$estvp
diagonal_combs <-
combs %>% dplyr::transmute(estvc = estvc,
estvp = estvp,
estvp2 = estvp)
combs_df <- dplyr::bind_rows(combs, combs2, diagonal_combs) %>%
dplyr::distinct(estvc, estvp, estvp2, .keep_all = TRUE) %>%
dplyr::transmute(Vc = estvc, Vp1 = estvp, Vp2 = estvp2)
} else if (sim_vp$mode == "manual" && sim_vp2$mode == "manual") {
combs_df <-
expand.grid(Vc = vc_values,
Vp1 = sim_vp$values,
Vp2 = sim_vp2$values)
} else if (sim_vp$mode == "manual" && sim_vp2$mode == "auto") {
vp2_df <-
expand.grid(vc = vc_values, ratio = c(10, 5, 2, 1, 0.5, 0.2, 0.1)) %>%
dplyr::mutate(vp2 = vc / ratio)
combs_df <- expand.grid(Vc = vc_values, Vp1 = sim_vp$values) %>%
dplyr::left_join(vp2_df %>% dplyr::rename(Vc = vc), by = "Vc") %>%
dplyr::rename(Vp2 = vp2)
} else if (sim_vp$mode == "auto" && sim_vp2$mode == "manual") {
vp1_df <-
expand.grid(vc = vc_values, ratio = c(10, 5, 2, 1, 0.5, 0.2, 0.1)) %>%
dplyr::mutate(vp1 = vc / ratio)
combs_df <-
expand.grid(Vc = vc_values, Vp2 = sim_vp2$values) %>%
dplyr::left_join(vp1_df %>% dplyr::rename(Vc = vc), by = "Vc") %>%
dplyr::rename(Vp1 = vp1)
}
# --- sim_q ---
sim_q$mode <- tryCatch(
match.arg(sim_q$mode, choices = c("manual", "auto")),
error = function(e) {
stop("Invalid `sim_q$mode`. Must be 'manual' or 'auto'.", call. = FALSE)
}
)
sim_q$auto.strategy <- tryCatch(
match.arg(sim_q$auto.strategy, choices = c("scaled", "fixed")),
error = function(e) {
stop("Invalid `sim_q$auto.strategy`. Must be 'scaled' or 'fixed'.",
call. = FALSE)
}
)
q1_values <- if (sim_q$mode == "auto") {
cl_val <- sim_cl$values[1]
if (sim_q$auto.strategy == "scaled")
c(cl_val / 4, cl_val / 2, cl_val, cl_val * 2)
else
c(1, 10, 100)
} else {
if (is.null(sim_q$values) || all(is.na(sim_q$values))) {
stop("No candidate Q values available for parameter sweeping.",
call. = FALSE)
}
sim_q$values
}
# --- sim_q2 ---
sim_q2$mode <- tryCatch(
match.arg(sim_q2$mode, choices = c("manual", "auto")),
error = function(e) {
stop("Invalid `sim_q2$mode`. Must be 'manual' or 'auto'.", call. = FALSE)
}
)
sim_q2$auto.strategy <- tryCatch(
match.arg(sim_q2$auto.strategy, choices = c("scaled", "fixed")),
error = function(e) {
stop("Invalid `sim_q2$auto.strategy`. Must be 'scaled' or 'fixed'.",
call. = FALSE)
}
)
q2_strategy_scaled <- FALSE
q2_values <- if (sim_q2$mode == "auto") {
if (sim_q2$auto.strategy == "scaled") {
q2_strategy_scaled <- TRUE
NULL
} else {
c(1, 10, 100)
}
} else {
if (is.null(sim_q2$values) || all(is.na(sim_q2$values))) {
stop("No candidate Q2 values available for parameter sweeping.",
call. = FALSE)
}
sim_q2$values
}
# --- Param grid ---
if (q2_strategy_scaled) {
param_grid <- tidyr::crossing(
combs_df,
CL = sim_cl$values,
Q1 = q1_values,
Ka = if (route == "oral")
sim_ka$values
else
NA_real_
)
param_grid$Q2 <- param_grid$Q1
param_grid <-
param_grid %>% dplyr::select(Vc, Vp1, Vp2, Q1, Q2, CL, Ka)
} else {
param_grid <- tidyr::crossing(
combs_df,
CL = sim_cl$values,
Q1 = q1_values,
Q2 = q2_values,
Ka = if (route == "oral")
sim_ka$values
else
NA_real_
) %>% dplyr::select(Vc, Vp1, Vp2, Q1, Q2, CL, Ka)
}
# k21,k31 boundary was set (0.01,5)
param_grid <- param_grid%>%
dplyr::filter((
(Vp1 < 1000 & Q1 / Vp1 > 0.01 & Q1 / Vp1 <= 5) |
(Vp1 >= 1000 & Q1 / Vp1 > 0.001 & Q1 / Vp1 <= 5)) &
((Vp2 < 1000 & Q2 / Vp2 > 0.01 & Q2 / Vp2 <= 5) |
(Vp2 >= 1000 & Q2 / Vp2 > 0.001 & Q2 / Vp2 <= 5)))
# --- Simulations ---
start_time <- Sys.time()
if (isTRUE(verbose)) {
handler.lists <- list(
progressr::handler_progress(
format = ":message [:bar] :percent (:current/:total)",
width = 80
)
)
} else {
handler.lists <- list(progressr::handler_void())
}
results <- progressr::with_progress({
p <- progressr::progressor(steps = nrow(param_grid))
param_grid %>%
dplyr::mutate(row = dplyr::row_number()) %>%
purrr::pmap_dfr(function(Vc, Vp1, Vp2, Q1, Q2, CL, Ka, row) {
p(sprintf(
"Running simulation: Vc=%.0f Vp1=%.0f Vp2=%.0f",
Vc,
Vp1,
Vp2
))
sim_out <- if (route == "iv") {
suppressMessages(suppressWarnings(Fit_3cmpt_iv(dat[dat$EVID != 2, ], "rxSolve", CL, Vc, Vp1, Vp2, Q1, Q2, input.add = 0)))
} else {
suppressMessages(suppressWarnings(Fit_3cmpt_oral(dat[dat$EVID != 2, ], "rxSolve", Ka, CL, Vc, Vp1, Vp2, Q1, Q2, input.add = 0)))
}
met <-
metrics.(pred.x = sim_out$cp, obs.y = dat[dat$EVID == 0,]$DV)
elapsed <-
round(difftime(Sys.time(), start_time, units = "secs"), 2)
tibble::tibble(
Vc = Vc,
Vp1 = Vp1,
Vp2 = Vp2,
Q1 = Q1,
Q2 = Q2,
CL = CL,
Ka = if (route == "oral")
Ka
else
NA_real_,
APE = round(met[1], 2),
MAE = round(met[2], 2),
MAPE = round(met[3], 2),
RMSE = round(met[4], 2),
rRMSE1 = round(met[5], 2),
rRMSE2 = round(met[6], 2),
Cumulative.Time.Sec = as.numeric(elapsed)
)
})
},handlers = handler.lists)
return(results)
}
Try the nlmixr2autoinit package in your browser
Any scripts or data that you put into this service are public.
nlmixr2autoinit documentation built on Nov. 14, 2025, 1:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions sim_sens_1cmpt_mm
Documented in sim_sens_1cmpt_mm
#' Parameter sweeping for a one-compartment Michaelis-Menten model
#'
#' Performs parameter sweeping by varying pharmacokinetic parameters in a
#' one-compartment model with Michaelis-Menten elimination.
#'
#' @param dat A data frame containing raw time–concentration data in the
#' standard nlmixr2 format.
#' @param sim_vmax List specifying Vmax:
#' - mode: "manual" or "auto"
#' - values: numeric vector if mode = "manual"
#' - est.cl: required if mode = "auto"
#' @param sim_km List specifying Km:
#' - mode: "manual" or "auto"
#' - values: numeric vector if mode = "manual"
#' @param sim_vd List specifying Vd:
#' - mode: must be "manual"
#' - values: numeric vector
#' @param sim_ka List specifying Ka (oral route only):
#' - mode: must be "manual"
#' - values: numeric vector
#' @param route Dosing route, either "iv" or "oral". Default is "iv".
#'
#' @param verbose Logical (default = TRUE).
#' Controls whether progress information is displayed during parameter sweeping.
#' When TRUE, a dynamic progress bar is shown using the `progressr` package to
#' indicate simulation status and elapsed time.
#' When FALSE, progress output is suppressed and the function runs silently.
#'
#' @details
#' The function generates a parameter grid and performs model fitting for each
#' combination using `Fit_1cmpt_mm_iv`. Parameters can be specified manually or
#' automatically derived. Model predictions and fit metrics are computed for each
#' simulation to assess parameter sensitivity.
#'
#' @return A data frame containing parameter combinations and model fit metrics.
#'
#' @author Zhonghui Huang
#'
#' @seealso \code{\link{Fit_1cmpt_mm_iv}}, \code{\link{Fit_1cmpt_mm_oral}}
#' @examples
#' \donttest{
#' # Example: IV dosing scenario with automatic Vmax and Km
#' out <- sim_sens_1cmpt_mm(
#' dat = Bolus_1CPTMM[Bolus_1CPTMM$ID<50,],
#' sim_vmax = list(mode = "auto", est.cl = 4),
#' sim_km = list(mode = "auto"),
#' sim_vd = list(mode = "manual", values = 70),
#' sim_ka = list(mode = "manual", values = NA),
#' route = "iv"
#' )
#' head(out[out$rRMSE2==min(out$rRMSE2),])
#'}
#' @export
sim_sens_1cmpt_mm <- function(dat,
sim_vmax = list(mode = "auto",
values = NULL,
est.cl = NULL),
sim_km = list(mode = "auto",
values = NULL),
sim_vd = list(mode = "manual", values = NULL),
sim_ka = list(mode = "manual", values = NULL),
route = c("iv", "oral"),
verbose= TRUE) {
`%>%` <- magrittr::`%>%`
route <- tryCatch(
match.arg(route, choices = c("iv", "oral")),
error = function(e) {
stop(
sprintf(
"Invalid value for `route`: '%s'. Must be 'iv' or 'oral'.",
as.character(route)
),
call. = FALSE
)
}
)
sim_vmax$mode <- tryCatch(
match.arg(sim_vmax$mode, choices = c("manual", "auto")),
error = function(e) {
stop("Invalid `sim_vmax$mode`. Must be 'manual' or 'auto'.",
call. = FALSE)
}
)
sim_km$mode <- tryCatch(
match.arg(sim_km$mode, choices = c("manual", "auto")),
error = function(e) {
stop("Invalid `sim_km$mode`. Must be 'manual' or 'auto'.", call. = FALSE)
}
)
sim_vd$mode <- tryCatch(
match.arg(sim_vd$mode, choices = c("manual")),
error = function(e) {
stop("Invalid `sim_vd$mode`. Must be 'manual'.", call. = FALSE)
}
)
vd_values <- sim_vd$values
if (is.null(vd_values) || all(is.na(vd_values))) {
stop("No candidate Vd values available for parameter sweeping.",
call. = FALSE)
}
# Deduplicate similar values if values are close
vd_values <- stats::na.omit(vd_values) %>%
sort() %>%
tibble::tibble(value = .) %>%
dplyr::mutate(prev = dplyr::lag(value),
rel_diff = abs(value - prev) / prev) %>%
dplyr::filter(is.na(rel_diff) | rel_diff > 0.2) %>%
dplyr::pull(value)
if (route == "oral") {
sim_ka$mode <- tryCatch(
match.arg(sim_ka$mode, choices = c("manual")),
error = function(e) {
stop("Invalid `sim_ka$mode`. Must be 'manual'.", call. = FALSE)
}
)
if (is.null(sim_ka$values) || all(is.na(sim_ka$values))) {
stop("No candidate Ka values available for parameter sweeping (oral route).",
call. = FALSE)
}
}
# Handle Vmax/Km grid
if (sim_vmax$mode == "auto" || sim_km$mode == "auto") {
if (is.null(sim_vmax$est.cl)) {
stop("Auto mode for Vmax/Km requires `sim_vmax$est.cl`.",
call. = FALSE)
}
cl_values <- sim_vmax$est.cl
# Deduplicate similar values if valus are close
cl_values <- stats::na.omit(cl_values) %>%
sort() %>%
tibble::tibble(value = .) %>%
dplyr::mutate(prev = dplyr::lag(value),
rel_diff = abs(value - prev) / prev) %>%
dplyr::filter(is.na(rel_diff) | rel_diff > 0.2) %>%
dplyr::pull(value)
# obtain the observed cmax
dat.obs <- dat[dat$EVID == 0, ]
pop.cmax <- aggregate(dat.obs$DV,
list(dat.obs$ID),
FUN = max,
na.rm = TRUE)
cmax <- mean(pop.cmax$x, trim = 0.05, na.rm = TRUE)
km_range <- c(4, 2, 1, 0.5, 0.25, 0.125, 0.1, 0.05) * cmax
conc_range <- c(0.05, 0.1, 0.25, 0.5, 0.75) * cmax
combs <-
expand.grid(Km = km_range, Conc = conc_range, cl = cl_values)
combs$Vmax <- (combs$Km + combs$Conc) * combs$cl
param_grid <- combs %>% dplyr::select(Vmax, Km)
# Filter combinations similar to each other
keep <- rep(TRUE, nrow(param_grid))
for (i in 1:(nrow(param_grid) - 1)) {
if (!keep[i])
next
for (j in (i + 1):nrow(param_grid)) {
if (!keep[j])
next
vmax_diff <-
abs((param_grid$Vmax[i] - param_grid$Vmax[j]) / param_grid$Vmax[j])
km_diff <-
abs((param_grid$Km[i] - param_grid$Km[j]) / param_grid$Km[j])
if (vmax_diff <= 0.2 && km_diff <= 0.2) {
keep[j] <- FALSE
}
}
}
param_grid <- param_grid[keep, ]
} else {
if (is.null(sim_vmax$values) || all(is.na(sim_vmax$values))) {
stop("No candidate Vmax values available for parameter sweeping.",
call. = FALSE)
}
if (is.null(sim_km$values) || all(is.na(sim_km$values))) {
stop("No candidate Km values available for parameter sweeping.",
call. = FALSE)
}
param_grid <-
tidyr::crossing(Vmax = sim_vmax$values, Km = sim_km$values)
}
if (route == "oral") {
param_grid <-
tidyr::crossing(param_grid, Vd = vd_values, Ka = sim_ka$values)
} else {
param_grid <- tidyr::crossing(param_grid, Vd = vd_values) %>%
dplyr::mutate(Ka = NA_real_)
}
start_time <- Sys.time()
if (isTRUE(verbose)) {
handler.lists <- list(
progressr::handler_txtprogressbar(
format = ":message [:bar] :percent (:current/:total)",
width = 80
)
)
} else {
handler.lists <- list(progressr::handler_void())
}
sim.1cmpt.mm.results.all <- progressr::with_progress({
p <- progressr::progressor(steps = nrow(param_grid))
param_grid %>%
dplyr::mutate(row = dplyr::row_number()) %>%
purrr::pmap_dfr(function(Vmax, Km, Vd, Ka, row) {
p(sprintf("Running simulation: Vmax=%.2f, Km=%.2f", Vmax, Km))
sim_out <- if (route == "iv") {
suppressMessages(suppressWarnings(Fit_1cmpt_mm_iv(
data = dat[dat$EVID != 2,],
est.method = "rxSolve",
input.vmax = Vmax,
input.km = Km,
input.vd = Vd,
input.add = 0
)))
} else {
suppressMessages(suppressWarnings(Fit_1cmpt_mm_oral(
data = dat[dat$EVID != 2,],
est.method = "rxSolve",
input.ka = Ka,
input.vmax = Vmax,
input.km = Km,
input.vd = Vd,
input.add = 0
)))
}
met <-
metrics.(pred.x = sim_out$cp, obs.y = dat[dat$EVID == 0,]$DV)
elapsed <-
round(difftime(Sys.time(), start_time, units = "secs"), 2)
# rm(sim_out)
# gc()
tibble::tibble(
Vmax = Vmax,
Km = Km,
Vd = Vd,
Ka = if (route == "oral")
Ka
else
NA_real_,
APE = round(met[1], 2),
MAE = round(met[2], 2),
MAPE = round(met[3], 2),
RMSE = round(met[4], 2),
rRMSE1 = round(met[5], 2),
rRMSE2 = round(met[6], 2),
Cumulative.Time.Sec = as.numeric(elapsed)
)
})
},handlers = handler.lists)
return(sim.1cmpt.mm.results.all)
}
#' Parameter sweeping for a two-compartment pharmacokinetic model
#'
#' Performs parameter sweeping by varying pharmacokinetic parameters in a
#' two-compartment model under IV or oral dosing. Model fit is evaluated across
#' combinations of CL, Vc, Vp, Q, and Ka (oral only).
#'
#' @details
#' The function generates a parameter grid and performs model fitting for each
#' combination using `Fit_2cmpt_iv` or `Fit_2cmpt_oral`. Parameters can be
#' specified manually or automatically derived. Model predictions and fit metrics
#' are computed for each simulation to assess parameter sensitivity.
#'
#' @param dat Pharmacokinetic dataset.
#' @param sim_ka List specifying Ka (oral route only):
#' - mode: must be "manual"
#' - values: numeric vector
#' @param sim_cl List specifying clearance (CL):
#' - mode: must be "manual"
#' - values: numeric vector
#' @param sim_vc List specifying central volume (Vc):
#' - mode: must be "manual"
#' - values: numeric vector
#' @param sim_vp List specifying peripheral volume (Vp):
#' - mode: "manual" or "auto"
#' - values: numeric vector if manual
#' @param sim_q List specifying inter-compartmental clearance (Q):
#' - mode: "manual" or "auto"
#' - values: numeric vector if manual
#' @param route Dosing route, either "iv" or "oral". Default is "iv".
#'
#' @param verbose Logical (default = TRUE).
#' Controls whether progress information is displayed during parameter sweeping.
#' When TRUE, a dynamic progress bar is shown using the `progressr` package to
#' indicate simulation status and elapsed time.
#' When FALSE, progress output is suppressed and the function runs silently.
#'
#' @return A data frame containing parameter combinations with model fit metrics.
#'
#' @author Zhonghui Huang
#'
#' @seealso \code{\link{Fit_2cmpt_iv}}, \code{\link{Fit_2cmpt_oral}}
#' @examples
#' \donttest{
#' out <- sim_sens_2cmpt(
#' dat = Bolus_2CPT[Bolus_2CPT$ID<50,],
#' sim_cl = list(mode = "manual", values = 4),
#' sim_vc = list(mode = "manual", values = 50),
#' sim_vp = list(mode = "auto"),
#' sim_q = list(mode = "auto"),
#' sim_ka = list(mode = "manual", values = NA),
#' route = "iv",verbose=FALSE
#' )
#' head(out[out$rRMSE2==min(out$rRMSE2),])
#' }
#' @export
sim_sens_2cmpt <- function(dat,
sim_ka = list(mode = "manual", values = NULL),
sim_cl = list(mode = "manual", values = NULL),
sim_vc = list(mode = "manual", values = NULL),
sim_vp = list(mode = c("auto", "manual"), values = NULL),
sim_q = list(
mode = c("auto", "manual"),
values = NULL,
auto.strategy = c("scaled", "fixed")
),
route = c("iv", "oral"),
verbose=TRUE) {
`%>%` <- magrittr::`%>%`
# Safe check
# Defensive validation for `route`
route <- tryCatch(
match.arg(route, choices = c("iv", "oral")),
error = function(e) {
stop(
sprintf(
"Invalid value for `route`: '%s'. Must be one of: %s.",
as.character(route),
paste(shQuote(c("iv", "oral")), collapse = ", ")
),
call. = FALSE
)
}
)
# Safe check and set sim ka
sim_ka$mode <- tryCatch(
match.arg(sim_ka$mode, choices = c("manual")),
error = function(e) {
stop("Invalid `sim_ka$mode`. Must be 'manual'.", call. = FALSE)
}
)
if (route == "oral" &&
(is.null(sim_ka$values) ||
length(sim_ka$values) == 0 || all(is.na(sim_ka$values)))) {
stop("No candidate Ka values available for parameter sweeping.",
call. = FALSE)
}
# Defensive checks for sim_cl
sim_cl$mode <- tryCatch(
match.arg(sim_cl$mode, choices = c("manual")),
error = function(e) {
stop("Invalid `sim_cl$mode`. Must be 'manual'.", call. = FALSE)
}
)
if (is.null(sim_cl$values) ||
length(sim_cl$values) == 0 || all(is.na(sim_cl$values))) {
stop("No candidate CL values available for parameter sweeping.",
call. = FALSE)
}
# Safe check and set sim vc
sim_vc$mode <- tryCatch(
match.arg(sim_vc$mode, choices = c("manual")),
error = function(e) {
stop("Invalid `sim_vc$mode`. Must be 'manual'.", call. = FALSE)
}
)
vc_values <- sim_vc$values
if (is.null(vc_values) ||
length(vc_values) == 0 || all(is.na(vc_values))) {
stop("No candidate Vc values available for parameter sweeping.",
call. = FALSE)
}
# Deduplicate similar values if values are close
vc_values <- stats::na.omit(vc_values)%>%
sort() %>%
tibble::tibble(value = .) %>%
dplyr::mutate(prev = dplyr::lag(value),
rel_diff = abs(value - prev) / prev) %>%
dplyr::filter(is.na(rel_diff) | rel_diff > 0.2) %>%
dplyr::pull(value)
# Safe check and set sim vp
sim_vp$mode <- tryCatch(
match.arg(sim_vp$mode, choices = c("manual", "auto")),
error = function(e) {
stop("Invalid `sim_vp$mode`. Must be 'manual' or 'auto'.", call. = FALSE)
}
)
if (sim_vp$mode == "auto") {
vp_ratios <- c(10, 5, 2, 1, 0.5, 0.2, 0.1)
vp_combs <- expand.grid(vc = vc_values, ratio = vp_ratios)
vp_combs$vp <- vp_combs$vc / vp_combs$ratio
} else {
if (is.null(sim_vp$values) ||
length(sim_vp$values) == 0 || all(is.na(sim_vp$values))) {
stop("No candidate Vp values available for parameter sweeping.",
call. = FALSE)
}
vp_combs <- expand.grid(vc = vc_values, vp = sim_vp$values)
}
# Safe check and set sim q
sim_q$mode <- tryCatch(
match.arg(sim_q$mode, choices = c("manual", "auto")),
error = function(e) {
stop("Invalid `sim_q$mode`. Must be 'manual' or 'auto'.", call. = FALSE)
}
)
if (sim_q$mode == "auto") {
strategy <- tryCatch(
match.arg(sim_q$auto.strategy, choices = c("scaled", "fixed")),
error = function(e) {
stop("Invalid `sim_q$auto.strategy`. Must be 'scaled' or 'fixed'.",
call. = FALSE)
}
)
cl_value <- sim_cl$values[1]
if (strategy == "scaled") {
q_values <- c(cl_value / 4, cl_value / 2, cl_value, cl_value * 2)
} else {
q_values <- c(1, 10, 100)
}
} else {
if (is.null(sim_q$values) ||
length(sim_q$values) == 0 || all(is.na(sim_q$values))) {
stop("No candidate Q values available for parameter sweeping.",
call. = FALSE)
}
q_values <- sim_q$values
}
vp_combs <- vp_combs %>%
dplyr::rename(Vc = vc, Vp = vp)
if (route == "oral") {
param_grid <- tidyr::crossing(vp_combs,
Q = q_values,
CL = sim_cl$values,
Ka = sim_ka$values)
} else {
param_grid <- tidyr::crossing(vp_combs,
Q = q_values,
CL = sim_cl$values) %>%
dplyr::mutate(Ka = NA_real_)
}
param_grid <- param_grid %>%
dplyr::select(Vc, Vp, Q, CL, Ka)
# k21 boundary was set (0.01,5)
param_grid <- param_grid%>%
dplyr::filter(
(Vp < 1000 & Q / Vp > 0.01 & Q / Vp <= 5) |
(Vp >= 1000 & Q / Vp > 0.001 & Q / Vp <= 5)
)
# === Begin Simulation ===
start_time <- Sys.time()
if (isTRUE(verbose)) {
handler.lists <- list(
progressr::handler_progress(
format = ":message [:bar] :percent (:current/:total)",
width = 80
)
)
} else {
handler.lists <- list(progressr::handler_void())
}
sim.2cmpt.results.all <- progressr::with_progress({
p <- progressr::progressor(steps = nrow(param_grid))
param_grid %>%
dplyr::mutate(row = dplyr::row_number()) %>%
purrr::pmap_dfr(function(Vc, Vp, Q, CL, Ka, row) {
p(sprintf("Running simulation: Vc=%.2f, Vp=%.2f", Vc, Vp))
sim_out <- if (route == "iv") {
suppressMessages(suppressWarnings(Fit_2cmpt_iv(
data = dat[dat$EVID != 2,],
est.method = "rxSolve",
input.cl = CL,
input.vc2cmpt = Vc,
input.vp2cmpt = Vp,
input.q2cmpt = Q,
input.add = 0
)))
} else {
suppressMessages(suppressWarnings(Fit_2cmpt_oral(
data = dat[dat$EVID != 2,],
est.method = "rxSolve",
input.ka = Ka,
input.cl = CL,
input.vc2cmpt = Vc,
input.vp2cmpt = Vp,
input.q2cmpt = Q,
input.add = 0
)))
}
met <-
metrics.(pred.x = sim_out$cp, obs.y = dat[dat$EVID == 0,]$DV)
elapsed <-
round(difftime(Sys.time(), start_time, units = "secs"), 2)
# sim_out<-NULL
tibble::tibble(
Vc = Vc,
Vp = Vp,
Q = Q,
CL = CL,
Ka = if (route == "oral")
Ka
else
NA,
APE = round(met[1], 2),
MAE = round(met[2], 2),
MAPE = round(met[3], 2),
RMSE = round(met[4], 2),
rRMSE1 = round(met[5], 2),
rRMSE2 = round(met[6], 2),
Cumulative.Time.Sec = as.numeric(elapsed)
)
})
},handlers = handler.lists)
return(sim.2cmpt.results.all)
}
#' Parameter sweeping for a three-compartment pharmacokinetic model
#'
#' Performs parameter sweeping by varying pharmacokinetic parameters in a
#' three-compartment model under IV or oral dosing. Parameter combinations
#' include Vc, Vp1, Vp2, Q1, Q2, CL, and Ka (oral only).
#'
#' @details
#' The function generates a parameter grid and evaluates each combination
#' using `Fit_3cmpt_iv` or `Fit_3cmpt_oral`. Model predictions and fit metrics
#' are calculated for each simulation to assess parameter influence and
#' identify optimal regions of the parameter space. Parameters can be provided
#' manually or derived automatically.
#'
#' @param dat Pharmacokinetic dataset.
#' @param sim_vc List specifying Vc:
#' - mode: must be "manual"
#' - values: numeric vector
#' @param sim_vp List specifying Vp1:
#' - mode: "manual" or "auto"
#' - values: numeric vector if manual
#' @param sim_vp2 List specifying Vp2:
#' - mode: "manual" or "auto"
#' - values: numeric vector if manual
#' @param sim_q List specifying Q1:
#' - mode: "manual" or "auto"
#' - values: numeric vector if manual
#' - auto.strategy: "scaled" or "fixed" when auto
#' @param sim_q2 List specifying Q2:
#' - mode: "manual" or "auto"
#' - values: numeric vector if manual
#' - auto.strategy: "scaled" or "fixed" when auto
#' @param sim_cl List specifying CL:
#' - mode: must be "manual"
#' - values: numeric vector
#' @param sim_ka List specifying Ka (oral route only):
#' - mode: must be "manual"
#' - values: numeric vector
#' @param route Dosing route, either "iv" or "oral". Default is "iv".
#'
#' @param verbose Logical (default = TRUE).
#' Controls whether progress information is displayed during parameter sweeping.
#' When TRUE, a dynamic progress bar is shown using the `progressr` package to
#' indicate simulation status and elapsed time.
#' When FALSE, progress output is suppressed and the function runs silently.
#'
#' @return A data frame containing parameter combinations with model fit metrics.
#'
#' @author Zhonghui Huang
#'
#' @seealso \code{\link{Fit_3cmpt_iv}}, \code{\link{Fit_3cmpt_oral}}
#' @examples
#' \donttest{
#' out <- sim_sens_3cmpt(
#' dat = Bolus_2CPT,
#' sim_cl = list(mode = "manual", values = 4),
#' sim_vc = list(mode = "manual", values = 50),
#' sim_vp = list(mode = "auto"),
#' sim_vp2 = list(mode = "auto"),
#' sim_q = list(mode = "auto", auto.strategy = "scaled"),
#' sim_q2 = list(mode = "auto", auto.strategy = "scaled"),
#' route = "iv",verbose=FALSE
#' )
#' head(out[out$rRMSE2==min(out$rRMSE2),])
#' }
#'
#' @export
sim_sens_3cmpt <- function(dat,
sim_vc = list(mode = "manual", values = NULL),
sim_vp = list(mode = c("auto", "manual"), values = NULL),
sim_vp2 = list(mode = c("auto", "manual"), values = NULL),
sim_q = list(
mode = c("auto", "manual"),
values = NULL,
auto.strategy = c("scaled", "fixed")
),
sim_q2 = list(
mode = c("auto", "manual"),
values = NULL,
auto.strategy = c("scaled", "fixed")
),
sim_cl = list(mode = "manual", values = NULL),
sim_ka = list(mode = "manual", values = NULL),
route = c("iv", "oral"),
verbose=TRUE) {
# --- Route check ---
route <- tryCatch(
match.arg(route, choices = c("iv", "oral")),
error = function(e) {
stop(sprintf(
"Invalid `route`: '%s'. Must be one of: %s.",
as.character(route),
paste(shQuote(c("iv", "oral")), collapse = ", ")
),
call. = FALSE)
}
)
# --- sim_ka ---
sim_ka$mode <- tryCatch(
match.arg(sim_ka$mode, choices = c("manual")),
error = function(e) {
stop("Invalid `sim_ka$mode`. Must be 'manual'.", call. = FALSE)
}
)
if (route == "oral" &&
(is.null(sim_ka$values) ||
length(sim_ka$values) == 0 || all(is.na(sim_ka$values)))) {
stop("No candidate Ka values available for parameter sweeping.",
call. = FALSE)
}
# --- sim_cl ---
sim_cl$mode <- tryCatch(
match.arg(sim_cl$mode, choices = c("manual")),
error = function(e) {
stop("Invalid `sim_cl$mode`. Must be 'manual'.", call. = FALSE)
}
)
if (is.null(sim_cl$values) ||
length(sim_cl$values) == 0 || all(is.na(sim_cl$values))) {
stop("No candidate CL values available for parameter sweeping.",
call. = FALSE)
}
# --- sim_vc ---
sim_vc$mode <- tryCatch(
match.arg(sim_vc$mode, choices = c("manual")),
error = function(e) {
stop("Invalid `sim_vc$mode`. Must be 'manual'.", call. = FALSE)
}
)
vc_values <- sim_vc$values
if (is.null(vc_values) ||
length(vc_values) == 0 || all(is.na(vc_values))) {
stop("No candidate Vc values available for parameter sweeping.",
call. = FALSE)
}
vc_values <- stats::na.omit(vc_values) %>%
sort() %>%
tibble::tibble(value = .) %>%
dplyr::mutate(prev = dplyr::lag(value),
rel_diff = abs(value - prev) / prev) %>%
dplyr::filter(is.na(rel_diff) | rel_diff > 0.2) %>%
dplyr::pull(value)
# --- sim_vp & sim_vp2 ---
sim_vp$mode <- tryCatch(
match.arg(sim_vp$mode, choices = c("auto", "manual")),
error = function(e) {
stop("Invalid `sim_vp$mode`. Must be 'auto' or 'manual'.", call. = FALSE)
}
)
sim_vp2$mode <- tryCatch(
match.arg(sim_vp2$mode, choices = c("auto", "manual")),
error = function(e) {
stop("Invalid `sim_vp2$mode`. Must be 'auto' or 'manual'.", call. = FALSE)
}
)
if (sim_vp$mode == "manual" &&
(is.null(sim_vp$values) || all(is.na(sim_vp$values)))) {
stop("No candidate Vp1 values in manual mode.", call. = FALSE)
}
if (sim_vp2$mode == "manual" &&
(is.null(sim_vp2$values) || all(is.na(sim_vp2$values)))) {
stop("No candidate Vp2 values in manual mode.", call. = FALSE)
}
if (sim_vp$mode == "auto" && sim_vp2$mode == "auto") {
vc_vp_ratio_range <- c(10, 5, 2, 1, 0.5, 0.2, 0.1)
combs <-
expand.grid(estvc = vc_values, vc_vp_ratio = vc_vp_ratio_range)
combs$estvp <- signif(combs$estvc / combs$vc_vp_ratio)
combs2 <- combs
combs$estvp2 <- combs$estvc
combs2$estvp2 <- combs2$estvp
diagonal_combs <-
combs %>% dplyr::transmute(estvc = estvc,
estvp = estvp,
estvp2 = estvp)
combs_df <- dplyr::bind_rows(combs, combs2, diagonal_combs) %>%
dplyr::distinct(estvc, estvp, estvp2, .keep_all = TRUE) %>%
dplyr::transmute(Vc = estvc, Vp1 = estvp, Vp2 = estvp2)
} else if (sim_vp$mode == "manual" && sim_vp2$mode == "manual") {
combs_df <-
expand.grid(Vc = vc_values,
Vp1 = sim_vp$values,
Vp2 = sim_vp2$values)
} else if (sim_vp$mode == "manual" && sim_vp2$mode == "auto") {
vp2_df <-
expand.grid(vc = vc_values, ratio = c(10, 5, 2, 1, 0.5, 0.2, 0.1)) %>%
dplyr::mutate(vp2 = vc / ratio)
combs_df <- expand.grid(Vc = vc_values, Vp1 = sim_vp$values) %>%
dplyr::left_join(vp2_df %>% dplyr::rename(Vc = vc), by = "Vc") %>%
dplyr::rename(Vp2 = vp2)
} else if (sim_vp$mode == "auto" && sim_vp2$mode == "manual") {
vp1_df <-
expand.grid(vc = vc_values, ratio = c(10, 5, 2, 1, 0.5, 0.2, 0.1)) %>%
dplyr::mutate(vp1 = vc / ratio)
combs_df <-
expand.grid(Vc = vc_values, Vp2 = sim_vp2$values) %>%
dplyr::left_join(vp1_df %>% dplyr::rename(Vc = vc), by = "Vc") %>%
dplyr::rename(Vp1 = vp1)
}
# --- sim_q ---
sim_q$mode <- tryCatch(
match.arg(sim_q$mode, choices = c("manual", "auto")),
error = function(e) {
stop("Invalid `sim_q$mode`. Must be 'manual' or 'auto'.", call. = FALSE)
}
)
sim_q$auto.strategy <- tryCatch(
match.arg(sim_q$auto.strategy, choices = c("scaled", "fixed")),
error = function(e) {
stop("Invalid `sim_q$auto.strategy`. Must be 'scaled' or 'fixed'.",
call. = FALSE)
}
)
q1_values <- if (sim_q$mode == "auto") {
cl_val <- sim_cl$values[1]
if (sim_q$auto.strategy == "scaled")
c(cl_val / 4, cl_val / 2, cl_val, cl_val * 2)
else
c(1, 10, 100)
} else {
if (is.null(sim_q$values) || all(is.na(sim_q$values))) {
stop("No candidate Q values available for parameter sweeping.",
call. = FALSE)
}
sim_q$values
}
# --- sim_q2 ---
sim_q2$mode <- tryCatch(
match.arg(sim_q2$mode, choices = c("manual", "auto")),
error = function(e) {
stop("Invalid `sim_q2$mode`. Must be 'manual' or 'auto'.", call. = FALSE)
}
)
sim_q2$auto.strategy <- tryCatch(
match.arg(sim_q2$auto.strategy, choices = c("scaled", "fixed")),
error = function(e) {
stop("Invalid `sim_q2$auto.strategy`. Must be 'scaled' or 'fixed'.",
call. = FALSE)
}
)
q2_strategy_scaled <- FALSE
q2_values <- if (sim_q2$mode == "auto") {
if (sim_q2$auto.strategy == "scaled") {
q2_strategy_scaled <- TRUE
NULL
} else {
c(1, 10, 100)
}
} else {
if (is.null(sim_q2$values) || all(is.na(sim_q2$values))) {
stop("No candidate Q2 values available for parameter sweeping.",
call. = FALSE)
}
sim_q2$values
}
# --- Param grid ---
if (q2_strategy_scaled) {
param_grid <- tidyr::crossing(
combs_df,
CL = sim_cl$values,
Q1 = q1_values,
Ka = if (route == "oral")
sim_ka$values
else
NA_real_
)
param_grid$Q2 <- param_grid$Q1
param_grid <-
param_grid %>% dplyr::select(Vc, Vp1, Vp2, Q1, Q2, CL, Ka)
} else {
param_grid <- tidyr::crossing(
combs_df,
CL = sim_cl$values,
Q1 = q1_values,
Q2 = q2_values,
Ka = if (route == "oral")
sim_ka$values
else
NA_real_
) %>% dplyr::select(Vc, Vp1, Vp2, Q1, Q2, CL, Ka)
}
# k21,k31 boundary was set (0.01,5)
param_grid <- param_grid%>%
dplyr::filter((
(Vp1 < 1000 & Q1 / Vp1 > 0.01 & Q1 / Vp1 <= 5) |
(Vp1 >= 1000 & Q1 / Vp1 > 0.001 & Q1 / Vp1 <= 5)) &
((Vp2 < 1000 & Q2 / Vp2 > 0.01 & Q2 / Vp2 <= 5) |
(Vp2 >= 1000 & Q2 / Vp2 > 0.001 & Q2 / Vp2 <= 5)))
# --- Simulations ---
start_time <- Sys.time()
if (isTRUE(verbose)) {
handler.lists <- list(
progressr::handler_progress(
format = ":message [:bar] :percent (:current/:total)",
width = 80
)
)
} else {
handler.lists <- list(progressr::handler_void())
}
results <- progressr::with_progress({
p <- progressr::progressor(steps = nrow(param_grid))
param_grid %>%
dplyr::mutate(row = dplyr::row_number()) %>%
purrr::pmap_dfr(function(Vc, Vp1, Vp2, Q1, Q2, CL, Ka, row) {
p(sprintf(
"Running simulation: Vc=%.0f Vp1=%.0f Vp2=%.0f",
Vc,
Vp1,
Vp2
))
sim_out <- if (route == "iv") {
suppressMessages(suppressWarnings(Fit_3cmpt_iv(dat[dat$EVID != 2, ], "rxSolve", CL, Vc, Vp1, Vp2, Q1, Q2, input.add = 0)))
} else {
suppressMessages(suppressWarnings(Fit_3cmpt_oral(dat[dat$EVID != 2, ], "rxSolve", Ka, CL, Vc, Vp1, Vp2, Q1, Q2, input.add = 0)))
}
met <-
metrics.(pred.x = sim_out$cp, obs.y = dat[dat$EVID == 0,]$DV)
elapsed <-
round(difftime(Sys.time(), start_time, units = "secs"), 2)
tibble::tibble(
Vc = Vc,
Vp1 = Vp1,
Vp2 = Vp2,
Q1 = Q1,
Q2 = Q2,
CL = CL,
Ka = if (route == "oral")
Ka
else
NA_real_,
APE = round(met[1], 2),
MAE = round(met[2], 2),
MAPE = round(met[3], 2),
RMSE = round(met[4], 2),
rRMSE1 = round(met[5], 2),
rRMSE2 = round(met[6], 2),
Cumulative.Time.Sec = as.numeric(elapsed)
)
})
},handlers = handler.lists)
return(results)
}
Try the nlmixr2autoinit package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.