Use pharmpy create_basic_pk_model() as default instead of templates

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: pharmr.extra
 Title: Extension of pharmr (Pharmpy) functionality
-Version: 0.0.0.9034
+Version: 0.0.0.9035
 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"),

R/create_model.R

@@ -127,6 +127,7 @@ create_model <- function(
     auto_stack_encounters = TRUE,
     drop_input = NULL,
     mu_reference = "auto",
+    use_template = FALSE,
     settings = list(), # TBD
     verbose = FALSE
 ) {
@@ -158,17 +159,28 @@ create_model <- function(
   }
 
   ## Read base model
-  if(verbose) cli::cli_alert_info("Reading base model")
-  mod <- pharmr::read_model(
-    path = get_template_modelfile(route, n_cmt, force_ode)
-  )
-  if(!is.logical(force_ode)) {
-    if(is.numeric(force_ode) || is.integer(force_ode) || is.character(force_ode)) {
-      advan <- as.integer(force_ode)
-      if(advan == 9L || advan == 13L) { # default ADVAN for templates is 6, update if needed:
-        mod <- update_advan(mod, advan)
+  if(use_template) { ## Use built-in templates
+    if(verbose) cli::cli_alert_info("Reading base model")
+    mod <- pharmr::read_model(
+      path = get_template_modelfile(route, n_cmt, force_ode)
+    )
+    if(!is.logical(force_ode)) {
+      if(is.numeric(force_ode) || is.integer(force_ode) || is.character(force_ode)) {
+        advan <- as.integer(force_ode)
+        if(advan == 9L || advan == 13L) { # default ADVAN for templates is 6, update if needed:
+          mod <- update_advan(mod, advan)
+        }
       }
     }
+  } else { ## Use pharmr to create the base model (default)
+    mod <- pharmr::create_basic_pk_model(
+      administration = route
+    )
+    ## Convert to NONMEM format early so downstream functions (e.g. set_iiv_block)
+    ## can work with NONMEM control stream syntax
+    if(tool != "nlmixr") {
+      mod <- pharmr::convert_model(mod, "nonmem")
+    }
   }
 
   ## Absorption
@@ -257,6 +269,23 @@ create_model <- function(
       cli::cli_alert_warning("Could not compute initial estimates automatically, please check manually.")
     } else {
       inits <- stats::setNames(inits, paste0("POP_", names(inits)))
+      ## Map init names to actual model parameter names (e.g. POP_V -> POP_VC)
+      model_params <- mod$parameters$names
+      init_aliases <- list(
+        POP_V = c("POP_VC"),
+        POP_CL = c("POP_CLMM")
+      )
+      for(i in seq_along(inits)) {
+        nm <- names(inits)[i]
+        if(!nm %in% model_params && nm %in% names(init_aliases)) {
+          for(alias in init_aliases[[nm]]) {
+            if(alias %in% model_params) {
+              names(inits)[i] <- alias
+              break
+            }
+          }
+        }
+      }
       mod <- pharmr::set_initial_estimates(
         model = mod,
         inits = inits
@@ -271,6 +300,8 @@ create_model <- function(
       model = mod,
       to_format = "nlmixr"
     )
+  } else {
+    mod <- pharmr::convert_model(mod, "nonmem")
   }
 
   ## Estimation method
@@ -288,6 +319,7 @@ create_model <- function(
       estimation_method = first_method, 
       estimation_options = per_step_options[[1]]
     )
+    if(!is.null(tool_options$MAXEVAL)) tool_options$MAXEVAL <- NULL
   } else {
     tool_options <- list()
     cli::cli_alert_warning(paste0("Skipping estimation options for ", tool, ", since not supported by Pharmpy. Please set manually"))

R/create_pharmpy_model_from_list.R

@@ -8,12 +8,16 @@ create_pharmpy_model_from_list <- function(model_obj) {
   ## Pharmpy bug: datainfo not updated when using pharmar::set_dataset()
   ## So need to make sure the dataset is on file when loading the model
   code <- model_obj$code
-  tmpfile <- tempfile()
-  write.csv(model_obj$dataset, tmpfile, quote=F, row.names=F)
-  code <- stringr::str_replace(
-    code,
-    "\\$DATA ([\\/a-zA-Z0-9\\.]*)",
-    paste0("$DATA ", tmpfile)
-  )
+  ## Strip trailing blank lines/whitespace that cause pharmpy DatasetError
+  code <- sub("[\\s\\n]+$", "", code, perl = TRUE)
+  if(!is.null(model_obj$dataset) && nrow(model_obj$dataset) > 0) {
+    tmpfile <- tempfile()
+    write.csv(model_obj$dataset, tmpfile, quote=F, row.names=F)
+    code <- stringr::str_replace(
+      code,
+      "\\$DATA ([\\/a-zA-Z0-9\\.]*)",
+      paste0("$DATA ", tmpfile)
+    )
+  }
   model <- pharmr::read_model_from_string(code)
 }

R/find_pk_parameter.R

@@ -23,15 +23,24 @@ find_pk_parameter <- function(parameter, model) {
   ## then, try to find depending on advan
   advan <- get_advan(model)
   if(advan %in% c(1, 3, 11)) {
-    map <- list("V" = "V1", "Q" = "QP1", "V2" = "VP1", "V3" = "VP2")
+    map <- list("V" = c("V1", "VC"), "V1" = c("VC"), "Q" = c("QP1"), "V2" = c("VP1"), "V3" = c("VP2"))
   } else {
-    map <- list("V" = "V2", "Q" = "QP1", "V3" = "VP1", "V4" = "VP2")
+    map <- list("V" = c("V2", "VC"), "V1" = c("VC"), "Q" = c("QP1"), "V3" = c("VP1"), "V4" = c("VP2"))
   }
   if(is.null(map[[parameter]])) {
     cli::cli_warn("Could not find parameter {parameter} in model as {parameter}, nor under different name.")
     return(parameter)
   } else {
-    cli::cli_alert_info("Found parameter {parameter} in model as {map[[parameter]]}.")
-    return(map[[parameter]])
+    ## Try each candidate; prefer one that exists in model
+    for(candidate in map[[parameter]]) {
+      if(candidate %in% model_params) {
+        cli::cli_alert_info("Found parameter {parameter} in model as {candidate}.")
+        return(candidate)
+      }
+    }
+    ## If none found in model, return the first candidate as best guess
+    candidate <- map[[parameter]][1]
+    cli::cli_alert_info("Found parameter {parameter} in model as {candidate}.")
+    return(candidate)
   }
 }

R/set_iiv.R

@@ -27,16 +27,6 @@ set_iiv <- function(mod, iiv, iiv_type = "exp") {
     }
   }
 
-  ## Make sure iiv_type is a list
-  if(inherits(iiv_type, "character")) {
-    iiv_type_list <- list()
-    for(key in names(iiv)) {
-      iiv_type_list[[key]] <- iiv_type
-    }
-  } else {
-    iiv_type_list <- iiv_type
-  }
-
   if(!is.null(iiv)) {
     if(!inherits(iiv, "list")) {
       stop("`iiv` parameter should be a `list` or a `character` object.")
@@ -46,6 +36,83 @@ set_iiv <- function(mod, iiv, iiv_type = "exp") {
     ## Then, add univariate IIV (no BLOCKs yet)
     all_params <- get_defined_pk_parameters(mod)
     current <- get_parameters_with_iiv(mod)
+
+    ## Map user-provided parameter names to actual model parameter names
+    ## before computing set differences. Pharmpy's create_basic_pk_model uses
+    ## different naming than the templates (e.g. VC instead of V, QP1 instead
+    ## of Q, VP1 instead of V2, CLMM instead of CL for MM models).
+    param_aliases <- list(
+      V   = c("V1", "VC"),
+      V1  = c("VC"),
+      V2  = c("VP1"),
+      V3  = c("VP2"),
+      Q   = c("QP1"),
+      Q2  = c("QP1"),
+      Q3  = c("QP2"),
+      CL  = c("CLMM")
+    )
+    ## Get all model assignments for fallback matching
+    all_assignments <- tryCatch({
+      stmts <- mod$statements$to_dict()$statements
+      vapply(
+        Filter(function(s) s$class == "Assignment", stmts),
+        function(s) gsub("(Symbol\\(\\'|\\'\\))", "", s$symbol),
+        character(1)
+      )
+    }, error = function(e) character(0))
+
+    iiv_name_map <- stats::setNames(names(iiv), names(iiv))
+    for(nm in names(iiv_name_map)) {
+      if(stringr::str_detect(nm, "~")) next
+      if(nm %in% names(param_aliases)) {
+        ## Check if the original name is a real parameter (not just an alias)
+        has_real_param <- nm %in% current ||
+          (nm %in% all_assignments && length(mod$statements$find_assignment(nm)) > 0 &&
+           !as.character(mod$statements$find_assignment(nm)$expression) %in% param_aliases[[nm]])
+        if(!has_real_param) {
+          for(alias in param_aliases[[nm]]) {
+            if(alias %in% current || alias %in% all_assignments) {
+              iiv_name_map[nm] <- alias
+              break
+            }
+          }
+        }
+      }
+    }
+    ## Apply the mapping to iiv names (including correlation entries like "CL~V1")
+    new_iiv_names <- names(iiv)
+    for(i in seq_along(new_iiv_names)) {
+      nm <- new_iiv_names[i]
+      if(stringr::str_detect(nm, "~")) {
+        parts <- stringr::str_split(nm, "~")[[1]]
+        mapped_parts <- vapply(parts, function(p) {
+          if(p %in% names(iiv_name_map)) iiv_name_map[p] else p
+        }, character(1))
+        new_iiv_names[i] <- paste(mapped_parts, collapse = "~")
+      } else if(nm %in% names(iiv_name_map)) {
+        new_iiv_names[i] <- iiv_name_map[nm]
+      }
+    }
+    names(iiv) <- new_iiv_names
+
+    ## Make sure iiv_type is a list (built after name mapping so keys match)
+    if(inherits(iiv_type, "character")) {
+      iiv_type_list <- list()
+      for(key in names(iiv)) {
+        iiv_type_list[[key]] <- iiv_type
+      }
+    } else {
+      ## Map user-provided iiv_type keys to match mapped iiv names
+      iiv_type_list <- iiv_type
+      for(orig_nm in names(iiv_name_map)) {
+        mapped_nm <- iiv_name_map[orig_nm]
+        if(orig_nm != mapped_nm && orig_nm %in% names(iiv_type_list)) {
+          iiv_type_list[[mapped_nm]] <- iiv_type_list[[orig_nm]]
+          iiv_type_list[[orig_nm]] <- NULL
+        }
+      }
+    }
+
     iiv_goal <- names(iiv)[!stringr::str_detect(names(iiv), "~")]
     iiv_corr <- names(iiv)[stringr::str_detect(names(iiv), "~")]
     has_corr <- unique(unlist(stringr::str_split(iiv_corr, "~")))
@@ -59,16 +126,6 @@ set_iiv <- function(mod, iiv, iiv_type = "exp") {
       dplyr::mutate(parameter = .data$name) |>
       dplyr::mutate(correlation = .data$name %in% has_corr) |>
       dplyr::arrange(.data$reset, .data$correlation) # make sure to first do the parameters that don't need a reset, to avoid creating DUMMYOMEGA
-    for(i in seq_along(map$name)) {
-      key <- map$name[i]
-      if(key == "V" && (! "V" %in% all_params) && "V1" %in% all_params) {
-        map$parameter[i] <- "V1"
-      }
-      if(key == "Q" && (! "QP1" %in% all_params) && "QP1" %in% all_params) {
-        map$parameter[i] <- "QP1"
-      }
-      names(iiv)[key == names(iiv)] <- map$parameter[i]
-    }
     for(i in seq_along(map$parameter)) {
       key <- map$name[i]
       par <- map$parameter[i]
@@ -145,6 +202,10 @@ set_iiv_block <- function(
     omega,
     code[(max(omega_idx)+1):length(code)]
   )
+  ## Remove trailing blank lines that cause pharmpy DatasetError
+  while(length(new_code) > 0 && new_code[length(new_code)] == "") {
+    new_code <- new_code[-length(new_code)]
+  }
   temp <- list(
     code = paste0(new_code, collapse = "\n"),
     dataset = model$dataset,

tests/testthat/test-add_default_output_tables.R

@@ -8,7 +8,7 @@ test_that("adding both tables when none exist works", {
     EVID = c(1, 0, 0),
     MDV = c(1, 0, 0)
   )
-  mod <- create_model(route = "iv", data = dat, tables = NULL, verbose = FALSE)
+  mod <- create_model(route = "iv", data = dat, tables = NULL, use_template = TRUE, verbose = FALSE)
 
   # Model has no tables initially:
   expect_length(get_tables_in_model_code(mod$code), 0)
@@ -39,7 +39,7 @@ test_that("adding only parameters table works", {
     EVID = c(1, 0, 0),
     MDV = c(1, 0, 0)
   )
-  mod <- create_model(route = "iv", data = dat, tables = NULL, verbose = FALSE)
+  mod <- create_model(route = "iv", data = dat, tables = NULL, use_template = TRUE, verbose = FALSE)
 
   # Model has no tables initially:
   expect_length(get_tables_in_model_code(mod$code), 0)
@@ -67,7 +67,7 @@ test_that("adding only fit table works", {
     EVID = c(1, 0, 0),
     MDV = c(1, 0, 0)
   )
-  mod <- create_model(route = "iv", data = dat, tables = NULL, verbose = FALSE)
+  mod <- create_model(route = "iv", data = dat, tables = NULL, use_template = TRUE, verbose = FALSE)
 
   # Model has no tables initially:
   expect_length(get_tables_in_model_code(mod$code), 0)