RXR-2888: skip in case of fit error and add options to reduce precision

NAMESPACE

@@ -15,6 +15,7 @@ export(calculate_shrinkage)
 export(calculate_stats)
 export(compare_psn_execute_results)
 export(compare_psn_proseval_results)
+export(fit_options)
 export(group_by_dose)
 export(group_by_time)
 export(install_default_literature_model)

R/calculate_stats.R

@@ -20,6 +20,16 @@ calculate_stats <- function(
   if(inherits(.res, "mipdeval_results")) {
     .res <- .res$results
   }
+  ## Check for errors during fits / predictions
+  errors <- dplyr::filter(
+    .res,
+    is.na(.data$pred) |
+    (is.na(.data$map_ipred) & !.data$apriori) |
+    is.na(.data$iter_ipred)
+  )
+  if(nrow(errors) > 0) {
+    cli::cli_warn("Errors were encountered in {nrow(errors)} out of {nrow(.res)} evaluated predictions. The problems occurred in patient(s) {unique(errors$id)}.")
+  }
   out <- .res |>
     tidyr::pivot_longer(
       cols = c("pred", "map_ipred", "iter_ipred"), names_to = "type"

R/run_eval.R

@@ -32,6 +32,8 @@
 #'   result from a call to [stats_summ_options()].
 #' @param .vpc_options Options for VPC simulations. This must be the result from
 #'   a call to [vpc_options()].
+#' @param .fit_options Options for controlling MAP Bayesian fit. This must be
+#'   the result from a call to [fit_options()].
 #' @param threads number of threads to divide computations on. Default is 1,
 #'   i.e. no parallel execution
 #' @param ruv residual error variability magnitude, specified as list.
@@ -63,6 +65,7 @@ run_eval <- function(
   incremental = FALSE,
   .stats_summ_options = stats_summ_options(),
   .vpc_options = vpc_options(),
+  .fit_options = fit_options(),
   threads = 1,
   progress = TRUE,
   verbose = TRUE
@@ -127,6 +130,7 @@ run_eval <- function(
     weight_prior = weight_prior,
     incremental = incremental,
     progress_function = p,
+    .fit_options = .fit_options,
     .threads = threads,
     .skip = .vpc_options$vpc_only
   )
@@ -189,3 +193,26 @@ run_eval <- function(
   ## 5. Return results
   out
 }
+
+#' Options for controlling MAP Bayesian fit
+#'
+#' @param ... These dots are reserved for future extensibility and must be empty.
+#' @param reltol Relative convergence tolerance. `reltol = 1e-04` will perform a
+#'   slightly coarser but more stable fit, which can be useful in some case.
+#'
+#' @returns A list.
+#' @export
+fit_options <- function(
+    ...,
+    reltol = 1e-05
+) {
+  rlang::check_dots_empty()
+  out <- list(
+    control = list(
+      reltol = vctrs::vec_assert(
+        reltol, ptype = numeric(), size = 1L, arg = "reltol"
+      )
+    )
+  )
+  structure(out, class = "mipdeval_fit_options")
+}

R/run_eval_core.R

@@ -14,7 +14,8 @@ run_eval_core <- function(
   weight_prior = 1,
   censor_covariates = TRUE,
   incremental = FALSE,
-  progress_function = function() {}
+  progress_function = function() {},
+  .fit_options = NULL
 ) {
 
   progress_function()
@@ -52,20 +53,26 @@ run_eval_core <- function(
       mod_upd$parameters <- fit$parameters # take params from previous fit
       mod_upd$omega <- fit$vcov
     }
-    fit <- PKPDmap::get_map_estimates(
-      model = mod_obj$model,
-      parameters = mod_upd$parameters,
-      omega = mod_upd$omega,
-      error = mod_obj$ruv,
-      fixed = mod_obj$fixed,
-      as_eta = mod_obj$kappa,
-      data = data$observations,
-      covariates = cov_data,
-      regimen = data$regimen,
-      weight_prior = weight_prior,
-      weights = weights,
-      iov_bins = mod_obj$bins,
-      verbose = FALSE
+    fit <- do.call(
+      PKPDmap::get_map_estimates,
+      c(
+        list(
+          model = mod_obj$model,
+          parameters = mod_upd$parameters,
+          omega = mod_upd$omega,
+          error = mod_obj$ruv,
+          fixed = mod_obj$fixed,
+          as_eta = mod_obj$kappa,
+          data = data$observations,
+          covariates = cov_data,
+          regimen = data$regimen,
+          weight_prior = weight_prior,
+          weights = weights,
+          iov_bins = mod_obj$bins,
+          verbose = FALSE
+        ),
+        .fit_options
+      )
     )
 
     ## Data frame with predictive data
@@ -85,9 +92,44 @@ run_eval_core <- function(
       `_iteration` = iterations[i],
       `_grouper` = obs_data$`_grouper`
     )
+    if(inherits(fit, "error")) {
+      ## create NA records for this fit
+      pred_data <- tibble::tibble(
+        id = obs_data$id,
+        t = obs_data$t,
+        dv = NA,
+        ipred = NA,
+        pred = NA,
+        ofv = NA,
+        ss_w = NA,
+        `_iteration` = iterations[i],
+        `_grouper` = obs_data$`_grouper`
+      )
+      par_dummy <- as.data.frame(mod_upd$parameters)
+      par_dummy[, 1:ncol(par_dummy)] <- NA
+      fit_pars <- dplyr::mutate(as.data.frame(par_dummy), id = obs_data$id[1])
+    } else {
+      ## Data frame with predictive data
+      pred_data <- tibble::tibble(
+        id = obs_data$id,
+        t = obs_data$t,
+        dv = fit$dv,
+        ipred = fit$ipred,
+        ires = fit$ires,
+        iwres = fit$iwres,
+        pred = fit$pred,
+        res = fit$res,
+        wres = fit$wres,
+        cwres = fit$cwres,
+        ofv = fit$fit$value,
+        ss_w = ss(fit$dv, fit$ipred, weights),
+        `_iteration` = iterations[i],
+        `_grouper` = obs_data$`_grouper`
+      )
+      ## Add parameter estimates
+      fit_pars <- dplyr::mutate(as.data.frame(fit$parameters), id = obs_data$id[1])
+    }
 
-    ## Add parameter estimates
-    fit_pars <- dplyr::mutate(as.data.frame(fit$parameters), id = obs_data$id[1])
     comb <- dplyr::bind_rows(
       comb, dplyr::left_join(pred_data, fit_pars, by = "id")
     )