diff --git a/NAMESPACE b/NAMESPACE
index 916d326..683b32d 100644
--- a/NAMESPACE
+++ b/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)
diff --git a/R/calculate_stats.R b/R/calculate_stats.R
index bfb3e66..308d59e 100644
--- a/R/calculate_stats.R
+++ b/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"
diff --git a/R/run_eval.R b/R/run_eval.R
index 23d738f..c04c3cf 100644
--- a/R/run_eval.R
+++ b/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")
+}
diff --git a/R/run_eval_core.R b/R/run_eval_core.R
index a730c34..7f105a1 100644
--- a/R/run_eval_core.R
+++ b/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")
     )
diff --git a/man/fit_options.Rd b/man/fit_options.Rd
new file mode 100644
index 0000000..47f320c
--- /dev/null
+++ b/man/fit_options.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/run_eval.R
+\name{fit_options}
+\alias{fit_options}
+\title{Options for controlling MAP Bayesian fit}
+\usage{
+fit_options(..., reltol = 1e-05)
+}
+\arguments{
+\item{...}{These dots are reserved for future extensibility and must be empty.}
+
+\item{reltol}{Relative convergence tolerance. \code{reltol = 1e-04} will perform a
+slightly coarser but more stable fit, which can be useful in some case.}
+}
+\value{
+A list.
+}
+\description{
+Options for controlling MAP Bayesian fit
+}
diff --git a/man/run_eval.Rd b/man/run_eval.Rd
index 1c49a8d..9c0c94c 100644
--- a/man/run_eval.Rd
+++ b/man/run_eval.Rd
@@ -21,6 +21,7 @@ run_eval(
   incremental = FALSE,
   .stats_summ_options = stats_summ_options(),
   .vpc_options = vpc_options(),
+  .fit_options = fit_options(),
   threads = 1,
   progress = TRUE,
   verbose = TRUE
@@ -88,6 +89,9 @@ result from a call to \code{\link[=stats_summ_options]{stats_summ_options()}}.}
 \item{.vpc_options}{Options for VPC simulations. This must be the result from
 a call to \code{\link[=vpc_options]{vpc_options()}}.}
 
+\item{.fit_options}{Options for controlling MAP Bayesian fit. This must be
+the result from a call to \code{\link[=fit_options]{fit_options()}}.}
+
 \item{threads}{number of threads to divide computations on. Default is 1,
 i.e. no parallel execution}
 
diff --git a/man/run_eval_core.Rd b/man/run_eval_core.Rd
index 23de731..a9c6020 100644
--- a/man/run_eval_core.Rd
+++ b/man/run_eval_core.Rd
@@ -13,7 +13,8 @@ run_eval_core(
   censor_covariates = TRUE,
   incremental = FALSE,
   progress_function = function() {
- }
+ },
+  .fit_options = NULL
 )
 }
 \arguments{
@@ -46,6 +47,9 @@ approach has been called "model predictive control (MPC)"
 "regular" MAP in some scenarios. Default is \code{FALSE}.}
 
 \item{progress_function}{function to increment progress bar}
+
+\item{.fit_options}{Options for controlling MAP Bayesian fit. This must be
+the result from a call to \code{\link[=fit_options]{fit_options()}}.}
 }
 \value{
 a \code{data.frame} with individual predictions