diff --git a/DESCRIPTION b/DESCRIPTION index cb94154..e5aa79b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: pharmr.extra Title: Extension of pharmr (Pharmpy) functionality -Version: 0.0.0.9072 +Version: 0.0.0.9081 Authors@R: c( person("Ron", "Keizer", email = "ron@insight-rx.com", role = c("cre", "aut")), person("Michael", "McCarthy", email = "michael.mccarthy@insight-rx.com", role = "ctb"), @@ -35,6 +35,7 @@ Imports: tidyr, tidyselect Suggests: + callr, mockery, nlmixr2, nlmixr2data, diff --git a/R/create_model.R b/R/create_model.R index bdf4d8d..c61c301 100644 --- a/R/create_model.R +++ b/R/create_model.R @@ -404,10 +404,15 @@ create_model <- function( add_sir(options = sir_options) } - ## MU referencing? - apply_mu <- (isTRUE(mu_reference) || - (identical(mu_reference, "auto") && "saem" %in% tolower(estimation_method))) && - tool == "nonmem" + ## MU referencing? Pharmpy supports MU-referencing for both NONMEM and + ## nlmixr2 backends. For nlmixr2 SAEM the rewrite is critical for stability + ## — pharmpy emits `mu_n <- log(POP_X); X <- exp(ETA_X + mu_n)` which keeps + ## the parameter on log scale during the MCMC E-step, preventing the M-step + ## from drifting positive THETAs through zero (the source of the lsoda + ## warning cascade and POP_* divergence we used to see on the admiral + ## SAEM run). + apply_mu <- isTRUE(mu_reference) || + (identical(mu_reference, "auto") && "saem" %in% tolower(estimation_method)) if(apply_mu && !pharmr::has_mu_reference(mod)) { mod <- pharmr::mu_reference_model(mod) } @@ -563,6 +568,13 @@ create_model <- function( if(!is.null(original_data)) { attr(mod, "original_data") <- original_data } + ## For nlmixr, stash the stacked/standard-name dataset instead — rxode2 + ## references columns by name and cannot auto-stack encounters, so the raw + ## pre-stacking input would conflate periods in a crossover trial and drive + ## SAEM into the lsoda step-limit cascade. + if(tool == "nlmixr" && inherits(data, "data.frame")) { + attr(mod, "original_data") <- data + } ## For nlmixr models, cache an SAEM-safe rewrite of the generated R code as ## an attribute so run_nlme() / run_sim() can use it verbatim. Pharmpy emits @@ -572,7 +584,7 @@ create_model <- function( ## internal state). NB: subsequent pharmpy operations on this model will ## drop this attribute; re-run create_model() if the model is modified. if(tool == "nlmixr") { - nlmixr_code <- inline_nlmixr_residual_aliases(mod$code) + nlmixr_code <- make_nlmixr_saem_safe(mod$code) if(!is.null(scale_observations)) { nlmixr_code <- inject_nlmixr_scaling(nlmixr_code, scale_observations) } diff --git a/R/run_nlme.R b/R/run_nlme.R index 4e4aadc..9cc524b 100644 --- a/R/run_nlme.R +++ b/R/run_nlme.R @@ -167,6 +167,7 @@ run_nlme <- function( save_summary = save_summary, save_final = save_final, clean = clean, + mu_reference = mu_reference, verbose = verbose )) } diff --git a/R/run_nlme_nlmixr.R b/R/run_nlme_nlmixr.R index 9f90e69..3cf8916 100644 --- a/R/run_nlme_nlmixr.R +++ b/R/run_nlme_nlmixr.R @@ -25,6 +25,7 @@ run_nlme_nlmixr <- function( save_summary = TRUE, save_final = TRUE, clean = TRUE, + mu_reference = "auto", verbose = TRUE ) { if(!requireNamespace("nlmixr2", quietly = TRUE)) { @@ -36,6 +37,29 @@ run_nlme_nlmixr <- function( time_start <- Sys.time() + ## MU-referencing — critical for nlmixr2 SAEM stability. Pharmpy's + ## mu_reference_model() rewrites `CL <- POP_CL * exp(ETA_CL)` as + ## `mu_1 <- log(POP_CL); CL <- exp(ETA_CL + mu_1)`, putting the parameter + ## walk on log scale so the M-step can't drift positive THETAs negative. + ## Idempotent (the `!has_mu_reference` check guards re-application). + est_for_mu <- estimation_method %||% tryCatch( + model$execution_steps$to_dataframe()$method, + error = function(e) NULL + ) + is_saem <- !is.null(est_for_mu) && "saem" %in% tolower(est_for_mu) + is_mu_ref <- isTRUE(pharmr::has_mu_reference(model)) + if((isTRUE(mu_reference) || (identical(mu_reference, "auto") && is_saem)) && !is_mu_ref) { + if(verbose) cli::cli_alert_info("Applying mu-referencing to model.") + model <- pharmr::mu_reference_model(model) + ## The mu-ref rewrite changes model$code, so invalidate any cached + ## SAEM-safe code (it was computed against the un-mu-referenced form). + attr(model, "nlmixr_code") <- NULL + } else if(isFALSE(mu_reference) && is_saem && !is_mu_ref) { + cli::cli_warn( + "nlmixr2 SAEM benefits significantly from mu-referencing — without it the M-step can drift positive THETAs through zero. Consider {.code mu_reference = \"auto\"}." + ) + } + ## Resolve dataset: prefer explicit `data`, then `attr(model, 'original_data')`, ## then `model$dataset`. nlmixr2 wants a data.frame in memory. fit_data <- resolve_nlmixr_data(model, data) @@ -46,11 +70,15 @@ run_nlme_nlmixr <- function( ## create_vpc_data_nlmixr() pick it up when the saved fit is re-used. if(!is.null(data)) attr(model, "original_data") <- fit_data - ## Use the SAEM-safe code cached by create_model(); fall back to the - ## pharmpy-generated $code verbatim for models built outside create_model(). - ## Note: pharmpy's default residual-alias pattern is not SAEM-compatible, - ## so SAEM fits on a BYO model will hit the upstream nlmixr2 error. - model_code <- attr(model, "nlmixr_code") %||% model$code + ## Use the SAEM-safe code cached by create_model() when present; otherwise + ## re-apply the residual-alias cleanup + residual-bound enforcement on + ## $code. The cached attribute can be lost across subsequent pharmpy ops + ## (e.g. update_parameters), and models built outside create_model() never + ## had it. Without this fallback SAEM hits `endpoint 'Y' for saem cannot + ## locate the residual error(s) correctly` (alias issue) or drifts the + ## residual σ negative (unbounded-init issue). + model_code <- attr(model, "nlmixr_code") %||% + make_nlmixr_saem_safe(model$code) ## Build a fresh run folder (mirrors NONMEM layout — dataset.csv + ## run.R holding the nlmixr2 function). @@ -78,22 +106,16 @@ run_nlme_nlmixr <- function( cli::cli_process_start(paste0("Starting nlmixr2 run in ", fit_folder)) } log_path <- file.path(fit_folder, output_file) - fit_args <- list(object = nlmixr_fn, data = fit_data, est = est) - if(!is.null(control)) fit_args$control <- control - ## Send output to both a log file and the console: - log_con <- file(log_path, open = "wt") - sink(log_con, type = "output", split = TRUE) - sink(log_con, type = "message") - on.exit({ - sink(type = "message") - sink(type = "output", split = TRUE) - close(log_con) - }, add = TRUE) - raw_fit <- do.call(nlmixr2::nlmixr2, fit_args) - sink(type = "message") - sink(type = "output", split = TRUE) - close(log_con) - on.exit() + lsoda_log_path <- file.path(fit_folder, "run_lsoda.log") + raw_fit <- run_nlmixr2_in_subprocess( + fn = nlmixr_fn, + data = fit_data, + est = est, + control = control, + log_path = log_path, + lsoda_log_path = lsoda_log_path, + verbose = verbose + ) if(verbose) cli::cli_process_done() ## Build a uniform fit object (pharmpy-shaped). @@ -112,9 +134,13 @@ run_nlme_nlmixr <- function( final_model <- update_parameters(model, fit) if(!is.null(final_model)) { if(!is.null(data)) attr(final_model, "original_data") <- fit_data + ## update_parameters() returns a fresh pharmpy object — re-cache the + ## SAEM-safe code so any later run_nlme()/run_sim() on the final model + ## doesn't fall back to the raw alias pattern. + attr(final_model, "nlmixr_code") <- make_nlmixr_saem_safe(final_model$code) attr(fit, "final_model") <- final_model if(save_final) { - writeLines(final_model$code, file.path(fit_folder, "final.R")) + writeLines(attr(final_model, "nlmixr_code"), file.path(fit_folder, "final.R")) } } @@ -142,6 +168,73 @@ run_nlme_nlmixr <- function( fit } +#' Fit an nlmixr2 model in a child R process and capture all output +#' +#' `nlmixr2`/`rxode2` ODE warnings (notably the `EE:[lsoda] / intdy --` +#' cascade) are emitted by C/Fortran code that writes directly to OS-level +#' stdout/stderr, bypassing R's `sink()`. Running in-process with +#' `sink(type = "message")` only catches R-level messages and warnings, so +#' those lsoda lines stream to the user's terminal but never reach the +#' run.log. We instead fit in a child R process via [callr::r()] with +#' `stderr = "2>&1"`, which redirects the child's OS-level FDs into the +#' parent so we can route each line ourselves. +#' +#' Output routing: each captured line is sorted in real time: +#' * lsoda/intdy/`@(lsoda.c:` lines → `run_lsoda.log` (these can run into +#' the thousands during SAEM exploration and would otherwise swamp the +#' main log) +#' * everything else (iteration trace, status messages, R warnings) → +#' `run.log` +#' * if `verbose = TRUE`, non-lsoda lines are also echoed to the parent +#' console (preserves the split-to-console feel of the old sink path) +#' +#' @noRd +run_nlmixr2_in_subprocess <- function(fn, data, est, control, log_path, + lsoda_log_path, verbose) { + if(!requireNamespace("callr", quietly = TRUE)) { + cli::cli_abort( + c("Package {.pkg callr} is required to run nlmixr2 with full output capture.", + i = "Install with {.code install.packages(\"callr\")}.") + ) + } + + ## Open both log files for the lifetime of the fit. They are closed on + ## exit even if the child errors. + main_con <- file(log_path, open = "wt") + lsoda_con <- file(lsoda_log_path, open = "wt") + on.exit({ close(main_con); close(lsoda_con) }, add = TRUE) + + ## Pattern matching anything that should go to the lsoda log instead of + ## the main log. Covers all three forms we've observed: + ## "unhandled error message: EE:[lsoda] ..." + ## "intdy -- t = ... illegal ..." + ## " @(lsoda.c:" (the trailing source-location line) + lsoda_re <- "EE:\\[?lsoda|^\\s*intdy --|^\\s*@\\(lsoda\\.c" + route_line <- function(line) { + if(grepl(lsoda_re, line, perl = TRUE)) { + writeLines(line, lsoda_con) + } else { + writeLines(line, main_con) + if(isTRUE(verbose)) cat(line, "\n", sep = "") + } + } + + callr::r( + func = function(fn, data, est, control) { + ## Bring nlmixr2's DSL helpers (ini/model) into the child's search + ## path so the parsed model function evaluates the same as in-process. + suppressPackageStartupMessages(library(nlmixr2est)) + fit_args <- list(object = fn, data = data, est = est) + if(!is.null(control)) fit_args$control <- control + do.call(nlmixr2est::nlmixr2, fit_args) + }, + args = list(fn = fn, data = data, est = est, control = control), + stderr = "2>&1", + callback = route_line, + spinner = FALSE + ) +} + #' Resolve dataset for an nlmixr fit #' #' @noRd @@ -368,6 +461,33 @@ find_observation_rows <- function(data) { NULL } +#' Rewrite pharmpy-emitted nlmixr code into a SAEM-safe form +#' +#' Composes the SAEM-safety transforms applied to pharmpy's nlmixr output: +#' 1. [inline_nlmixr_residual_aliases()] — collapse `add_error`/`prop_error` +#' indirection so the residual terms in `~` reference `ini()` params +#' directly. +#' 2. [enforce_residual_bounds()] — stamp `c(0, init, Inf)` on the +#' residual-error params so SAEM can't drift them negative. +#' 3. [enforce_theta_bounds()] — stamp `c(0, init, Inf)` on unbounded +#' positive-init THETAs so SAEM's M-step can't push them through zero. +#' 4. [apply_ipred_guard()] — wire pharmpy's `IPREDADJ` floor-guard into +#' `Y` so proportional/combined residuals don't collapse at IPRED == 0. +#' +#' Each transform is a no-op when its expected pattern isn't present, so this +#' is safe to apply unconditionally on any pharmpy-emitted nlmixr code. +#' +#' @noRd +make_nlmixr_saem_safe <- function(code) { + apply_ipred_guard( + enforce_theta_bounds( + enforce_residual_bounds( + inline_nlmixr_residual_aliases(code) + ) + ) + ) +} + #' Inline pharmpy's residual-error aliases for SAEM compatibility #' #' Pharmpy's nlmixr converter always emits the residual block as @@ -399,6 +519,40 @@ inline_nlmixr_residual_aliases <- function(code) { } add_val <- trimws(sub("^\\s*add_error\\s*<-\\s*", "", lines[add_idx])) prop_val <- trimws(sub("^\\s*prop_error\\s*<-\\s*", "", lines[prop_idx])) + + ## nlmixr2's SAEM parser rejects expressions inside add()/prop() — each must + ## be a bare `ini()` parameter. The `use_template = TRUE` path emits + ## `add_error <- RUV_ADD*W` with `W <- 1` defined inside model({}), which + ## stays a multiplication after the alias substitution. Simplify + ## `*` (or `*`) when `` is assigned a numeric + ## constant in the same code (so RUV_ADD*W with W=1 collapses to RUV_ADD). + const_assigns <- regmatches( + lines, + regexec("^\\s*([A-Za-z_][A-Za-z0-9_]*)\\s*<-\\s*(-?[0-9.eE+-]+)\\s*$", lines) + ) + const_map <- list() + for(m in const_assigns) { + if(length(m) == 3) { + val <- suppressWarnings(as.numeric(m[[3]])) + if(!is.na(val)) const_map[[m[[2]]]] <- val + } + } + simplify_alias <- function(expr) { + m <- regmatches(expr, + regexec("^\\s*([A-Za-z_][A-Za-z0-9_]*)\\s*\\*\\s*([A-Za-z_][A-Za-z0-9_]*)\\s*$", + expr))[[1]] + if(length(m) != 3) return(expr) + lhs_const <- const_map[[m[[2]]]] + rhs_const <- const_map[[m[[3]]]] + if(!is.null(lhs_const) && lhs_const == 1) return(m[[3]]) + if(!is.null(rhs_const) && rhs_const == 1) return(m[[2]]) + if(!is.null(lhs_const) && lhs_const == 0) return("0") + if(!is.null(rhs_const) && rhs_const == 0) return("0") + expr + } + add_val <- simplify_alias(add_val) + prop_val <- simplify_alias(prop_val) + parts <- character(0) if(add_val != "0") parts <- c(parts, paste0("add(", add_val, ")")) if(prop_val != "0") parts <- c(parts, paste0("prop(", prop_val, ")")) @@ -409,6 +563,157 @@ inline_nlmixr_residual_aliases <- function(code) { paste(lines, collapse = "\n") } +#' Stamp `c(0, init, Inf)` bounds on residual-error params in nlmixr ini() +#' +#' Pharmpy emits THETAs with explicit `c(lower, init, upper)` bounds but +#' typically leaves the residual-error params unbounded — e.g. +#' \preformatted{ +#' sigma1 <- 0.01 # additive +#' sigma <- 0.09 # proportional +#' } +#' nlmixr2's SAEM is gradient-free and treats unbounded params as truly +#' unbounded, so `sigma1` can drift to negative values during the E-step. +#' Once σ goes negative the likelihood inverts sign, ETAs explode, and the +#' ODE solver gets driven into unphysical states (the source of the +#' `EE:[lsoda] / intdy --` warning cascade). pharmpy's downstream +#' `update_parameters()` also refuses to ingest the resulting negative inits +#' because every THETA has lower bound 0. +#' +#' This function scans the `ini({...})` block for any unbounded +#' `name <- ` line whose `name` appears inside an `add(...)` or +#' `prop(...)` call in the model formula (i.e. it's a residual-error param) +#' and rewrites it as `name <- c(0, , Inf)`. +#' +#' Run *after* [inline_nlmixr_residual_aliases()] so the residual params +#' show up as bare identifiers inside `add()`/`prop()` (not as `add_error` +#' aliases). Returns the input unchanged if no such params are found. +#' +#' @noRd +enforce_residual_bounds <- function(code) { + lines <- strsplit(code, "\n", fixed = TRUE)[[1]] + ini_start <- grep("^\\s*ini\\(", lines) + if(length(ini_start) == 0) return(code) + ini_close <- grep("^\\s*\\}\\)\\s*$", lines) + ini_close <- ini_close[ini_close > ini_start[1]] + if(length(ini_close) == 0) return(code) + ini_end <- ini_close[1] + + ## Collect identifiers used inside add(...) / prop(...) in the Y formula + y_idx <- grep("^\\s*Y\\s*~", lines) + if(length(y_idx) == 0) return(code) + residual_refs <- character(0) + call_re <- "(?:add|prop)\\(([^)]+)\\)" + for(idx in y_idx) { + calls <- regmatches(lines[idx], gregexpr(call_re, lines[idx], perl = TRUE))[[1]] + for(call in calls) { + inside <- sub("^[a-z]+\\(", "", call) + inside <- sub("\\)$", "", inside) + ## Split on operators / whitespace to grab identifier tokens + tokens <- strsplit(inside, "[*+\\-/\\s]+", perl = TRUE)[[1]] + idents <- tokens[grepl("^[A-Za-z_][A-Za-z0-9_]*$", tokens)] + residual_refs <- c(residual_refs, idents) + } + } + residual_refs <- unique(residual_refs) + if(length(residual_refs) == 0) return(code) + + ## Wrap unbounded numeric assignments inside ini({...}) + bound_re <- "^(\\s*)([A-Za-z_][A-Za-z0-9_]*)\\s*<-\\s*(-?[0-9.eE+-]+)\\s*$" + changed <- FALSE + for(i in seq(ini_start[1], ini_end)) { + m <- regmatches(lines[i], regexec(bound_re, lines[i]))[[1]] + if(length(m) == 4 && m[[3]] %in% residual_refs) { + ## guard: skip if init is already non-positive (would invert the bound) + val <- suppressWarnings(as.numeric(m[[4]])) + if(is.na(val) || val <= 0) next + lines[i] <- paste0(m[[2]], m[[3]], " <- c(0, ", m[[4]], ", Inf)") + changed <- TRUE + } + } + if(!changed) return(code) + paste(lines, collapse = "\n") +} + +#' Stamp `c(0, init, Inf)` bounds on unbounded positive-init THETAs +#' +#' Pharmpy's NONMEM backend emits `$THETA (0, init)` for population +#' parameters, but its nlmixr2 backend emits the bare `POP_CL <- init` form — +#' nlmixr2 SAEM treats unbounded params as truly unbounded, so the M-step can +#' push `POP_CL` (or any other PK theta) through zero into negative values. +#' Once that happens the structural model produces non-physical predictions +#' (e.g. negative CL means concentrations grow without bound), lsoda fails to +#' integrate, and the fit diverges (admiral SAEM on this codebase: POP_CL +#' crosses zero at iteration ~25 and reaches -7.7M by iteration 100). +#' +#' Scan the `ini({...})` block for any `name <- ` +#' line whose `name` is NOT already handled as a residual param by +#' [enforce_residual_bounds()], and rewrite as `name <- c(0, init, Inf)`. +#' Skips lines that already have bounds (`c(...)`, `fixed(...)`) and skips +#' negative inits to preserve intentional sign conventions. +#' +#' @noRd +enforce_theta_bounds <- function(code) { + lines <- strsplit(code, "\n", fixed = TRUE)[[1]] + ini_start <- grep("^\\s*ini\\(", lines) + if(length(ini_start) == 0) return(code) + ini_close <- grep("^\\s*\\}\\)\\s*$", lines) + ini_close <- ini_close[ini_close > ini_start[1]] + if(length(ini_close) == 0) return(code) + ini_end <- ini_close[1] + + bound_re <- "^(\\s*)([A-Za-z_][A-Za-z0-9_]*)\\s*<-\\s*(-?[0-9.eE+-]+)\\s*$" + changed <- FALSE + for(i in seq(ini_start[1], ini_end)) { + m <- regmatches(lines[i], regexec(bound_re, lines[i]))[[1]] + if(length(m) != 4) next + val <- suppressWarnings(as.numeric(m[[4]])) + if(is.na(val) || val <= 0) next # leave negative/zero inits alone + lines[i] <- paste0(m[[2]], m[[3]], " <- c(0, ", m[[4]], ", Inf)") + changed <- TRUE + } + if(!changed) return(code) + paste(lines, collapse = "\n") +} + +#' Wire pharmpy's `IPREDADJ` floor-guard into `Y` +#' +#' For proportional and combined residual error, pharmpy emits a floor-guard +#' block of the form +#' \preformatted{ +#' IPRED <- A_CENTRAL/VC +#' if (0 == IPRED) { +#' IPREDADJ <- 2.225e-16 +#' } else { +#' IPREDADJ <- IPRED +#' } +#' Y <- IPRED # <-- bug: should be IPREDADJ +#' Y ~ prop(sigma) +#' } +#' followed by `Y <- IPRED` — so `IPREDADJ` is declared but never used. With +#' `Y ~ prop(sigma)` against `Y = IPRED`, the residual variance collapses to +#' 0 wherever IPRED is 0 (e.g. before absorption, or for extreme ETAs during +#' SAEM's E-step). The likelihood inverts, ETAs blow up, lsoda gets driven +#' into unphysical states, and the `EE:[lsoda] / intdy --` warning cascade +#' floods stderr — what 0.0.0.9078's subprocess capture now lands in run.log. +#' +#' Rewrite `Y <- IPRED` to `Y <- IPREDADJ` when the guard block is present +#' (detected by `IPREDADJ <- IPRED` in the else branch). No-op when no guard +#' block exists (additive RUV — pharmpy correctly omits the block there). +#' +#' @noRd +apply_ipred_guard <- function(code) { + lines <- strsplit(code, "\n", fixed = TRUE)[[1]] + if(!any(grepl("^\\s*IPREDADJ\\s*<-\\s*IPRED\\s*$", lines))) return(code) + y_re <- "^(\\s*)Y\\s*<-\\s*IPRED\\s*$" + y_idx <- grep(y_re, lines) + if(length(y_idx) == 0) return(code) + for(i in y_idx) { + indent <- sub("\\S.*$", "", lines[i]) + lines[i] <- paste0(indent, "Y <- IPREDADJ") + } + paste(lines, collapse = "\n") +} + #' Inject `scale_observations` into pharmpy-generated nlmixr code #' #' Pharmpy emits the observation prediction as diff --git a/R/run_sim_nlmixr.R b/R/run_sim_nlmixr.R index b10d6ca..796ebe5 100644 --- a/R/run_sim_nlmixr.R +++ b/R/run_sim_nlmixr.R @@ -64,9 +64,11 @@ run_sim_nlmixr <- function( ## Extract the rxode2/nlmixr2 function from the pharmpy model code. Prefer ## the cached, post-processed code (set by create_model() — accounts for - ## SAEM-safe residual aliases and scale_observations) over pharmpy's - ## regenerated mod$code. - model_code <- attr(model, "nlmixr_code") %||% model$code + ## SAEM-safe residual aliases and scale_observations); otherwise re-apply + ## the residual-alias cleanup, since the cached attribute can be lost across + ## subsequent pharmpy ops (e.g. update_parameters returns a fresh object). + model_code <- attr(model, "nlmixr_code") %||% + make_nlmixr_saem_safe(model$code) nlmixr_fn <- extract_nlmixr_function(model_code) if(is.null(nlmixr_fn)) { cli::cli_abort("Could not extract an nlmixr2 model function from the model code.") diff --git a/R/set_iiv.R b/R/set_iiv.R index 986fe9c..4ec0ed9 100644 --- a/R/set_iiv.R +++ b/R/set_iiv.R @@ -119,6 +119,13 @@ set_iiv <- function(mod, iiv, iiv_type = "exp") { to_remove <- setdiff(current, iiv_goal) to_reset <- intersect(iiv_goal, current) to_add <- setdiff(iiv_goal, current) + + ## Drop IIV on parameters not in the requested spec (e.g. pharmpy's + ## create_basic_pk_model("oral") puts IIV on MAT by default; if the caller + ## only asks for IIV on CL and V, MAT should not carry one). + for(par in to_remove) { + mod <- pharmr::remove_iiv(mod, par) + } map <- data.frame( # build a map for each parameter, whether it needs to be reset or not name = c(to_add, to_reset), reset = c(rep(FALSE, length(to_add)), rep(TRUE, length(to_reset)))