Refactor of main MAP fit function

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: PKPDmap
 Type: Package
 Title: MAP Bayesian estimates
-Version: 1.0.1
-Date: 2024-06-10
+Version: 1.0.3
+Date: 2025-06-18
 Author: Ron Keizer, Jasmine Hughes, Kara Woo
 Maintainer: Ron Keizer <ron@insight-rx.com>
 Description: Obtain MAP Bayesian estimates based on individual data and
@@ -12,6 +12,6 @@ Depends: utils, stats
 Imports: PKPDsim, numDeriv, mvtnorm
 Suggests: testthat (>= 3.0.0)
 Remotes: InsightRX/PKPDsim
-RoxygenNote: 7.3.1
+RoxygenNote: 7.3.2
 Encoding: UTF-8
 Config/testthat/edition: 3

NEWS.md

@@ -0,0 +1,9 @@
+# Changes in master compared to last release (20240610)
+
+- many code chunks in main estimation function now separated off to their own function
+- improved handling of floating points during warning check
+- The print method for map_estimates() now uses 
+  `sqrt(diag(PKPDsim::triangle_to_full(x$prior$omega$full)))` instead of 
+  `sqrt(diag(PKPDsim::triangle_to_full(x$prior$omega)))` for etas
+- fix for mixture models: ipreds for mixture models were found not to be correct
+- censoring information is now returned with the fit object

PKPDmap.Rproj

@@ -1,5 +1,5 @@
 Version: 1.0
-ProjectId: 42e31d31-a19a-49a8-abcd-daac065467eb
+ProjectId: 6835582e-ba67-4423-80dd-8a468b12df41
 
 RestoreWorkspace: Default
 SaveWorkspace: Default

R/calc_residuals.R

@@ -0,0 +1,110 @@
+#' Function to calculate residuals from a fit,
+#' based on a fitted parameter set, and data object.
+#' 
+#' @inheritParams get_map_estimates
+#' @param obj temporary storage object in `get_map_estimates()`
+#' @param parameters_population parameters for population predictions
+#' @param omega_full full omega matrix
+#' @param transf transformation function
+#' @param A_init_population init vector for model state (population)
+#' @param A_init_individual init vector for model state (individual)
+#' @param censoring_idx vector with indices for censoring
+#' @param data_before_init data.frame with data before initial dose
+#' 
+calc_residuals <- function(
+  obj,
+  data,
+  model,
+  parameters_population,
+  covariates,
+  regimen,
+  omega_full,
+  error,
+  weights,
+  transf,
+  A_init_population,
+  A_init_individual,
+  t_init = 0,
+  iov_bins = NULL,
+  output_include = c(),
+  int_step_size = 0.01,
+  censoring = NULL,
+  censoring_idx = NULL,
+  data_before_init = NULL,
+  ltbs = FALSE,
+  ...
+) {
+
+  ## Observation vectors
+  t_obs <- c(data_before_init$t, data$t)
+  obs_type <- c(data_before_init$obs_type, data$obs_type)
+
+  ## Perform simulations for ipred and pred
+  suppressMessages({
+    ## After fitting individual parameters, don't pass the mixture group to the simulator 
+    ## (so `mixture=NULL`), otherwise `sim()` will use the population value for the 
+    ## specified group, and not the individual fitted parameter.
+    sim_ipred <- PKPDsim::sim_ode(
+      ode = model,
+      parameters = obj$parameters,
+      mixture_group = NULL,
+      covariates = covariates,
+      n_ind = 1,
+      int_step_size = int_step_size,
+      regimen = regimen,
+      t_obs = t_obs,
+      obs_type = obs_type,
+      only_obs = TRUE,
+      checks = FALSE,
+      A_init = A_init_individual,
+      iov_bins = iov_bins,
+      output_include = output_include,
+      t_init = t_init,
+      ...
+    )
+  })
+  suppressMessages({
+    sim_pred <- PKPDsim::sim_ode(
+      ode = model,
+      parameters = parameters_population,
+      mixture_group = NULL,
+      covariates = covariates,
+      n_ind = 1,
+      int_step_size = int_step_size,
+      regimen = regimen,
+      t_obs = c(data_before_init$t, t_obs),
+      obs_type = c(data_before_init$obs_type, data$obs_type),
+      only_obs = TRUE,
+      checks = FALSE,
+      iov_bins = iov_bins,
+      A_init = A_init_population,
+      t_init = t_init,
+      ...
+    )
+  })
+
+  ## Parse the residuals from the predictions and the data, and add to obj
+  obj <- parse_residuals_from_predictions(
+    obj,
+    sim_ipred,
+    sim_pred,
+    data,
+    omega_full,
+    transf,
+    error,
+    censoring,
+    censoring_idx,
+    data_before_init,
+    weights
+  )
+
+  ## Add covariates and parameters to obj
+  if(output_include$covariates && !is.null(covariates)) {
+    obj$covariates_time <- sim_ipred[!duplicated(sim_ipred$t), names(covariates)]
+  }
+  if(output_include$parameters) {
+    obj$parameters_time <- sim_ipred[!duplicated(sim_ipred$t), names(parameters)]
+  }
+
+  obj
+}

R/check_inputs.R

@@ -0,0 +1,17 @@
+#' Some basic input checking
+#' 
+#' @inheritParams get_map_estimates
+#' 
+check_inputs <- function(model, data, parameters, omega, regimen, censoring, type) {
+  if(tolower(type) %in% c("map", "pls")) {
+    if(is.null(model) || is.null(data) || is.null(parameters) || is.null(omega) || is.null(regimen)) {
+      stop("The 'model', 'data', 'omega', 'regimen', and 'parameters' arguments are required.")
+    }
+  }
+  if(!is.null(censoring) && !inherits(censoring, "character")) {
+    stop("Censoring argument requires label specifying column in dataset with censoring info.")
+  }
+  if(!("function" %in% class(model))) {
+    stop("The 'model' argument requires a function, e.g. a model defined using the new_ode_model() function from the PKPDsim package.")
+  }
+}

R/get_individual_parameters_from_fit.R

@@ -0,0 +1,22 @@
+#' Get the individual parameters from the fitted coefficients (etas)
+#' 
+#' @param fit fit object
+#' @param parameters list of population parameters
+#'
+get_individual_parameters_from_fit <- function(
+  fit,
+  parameters,
+  nonfixed = NULL,
+  as_eta = NULL
+) {
+  par <- parameters
+  for(i in seq(nonfixed)) {
+    key <- nonfixed[i]
+    if(key %in% as_eta) {
+      par[[key]] <- as.numeric(fit$coef[i])
+    } else {
+      par[[key]] <- as.numeric(as.numeric(par[[key]]) * exp(as.numeric(fit$coef[i])))
+    }
+  }
+  par
+}