Improve handling of (partially) failing fits

6 files changed, 48 insertions, 31 deletions

diff --git a/NEWS.md b/NEWS.md
index 4a923f09..25fe7f14 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -12,6 +12,10 @@
 
 - 'mkinfit': Run 'stats::shapiro.test()' on standardized residuals and warn if p < 0.05
 
+- 'mkinfit': 'error_model_algorithm' = 'd_3' does not fail if direct fitting fails, but reports that the results for the threestep algorithm are returned
+
+- 'mmkin': Do not fail any more if one of the fits fails, but assign the try-error to the respective position in the mmkin object
+
 # mkin 0.9.50.2 (2020-05-12)
 
 - Increase tolerance for a platform specific test results on the Solaris test machine on CRAN
diff --git a/R/mkinfit.R b/R/mkinfit.R
index 154c2a18..73fe43e0 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -712,17 +712,17 @@ mkinfit <- function(mkinmod, observed,
     if (error_model_algorithm == "d_3") {
       if (!quiet) message("Directly optimising the complete model")
       parms.start <- c(degparms, errparms)
-      fit_direct <- nlminb(parms.start, cost_function,
+      fit_direct <- try(nlminb(parms.start, cost_function,
         lower = lower[names(parms.start)],
         upper = upper[names(parms.start)],
-        control = control, ...)
-      fit_direct$logLik <- - cost.current
-      if (error_model_algorithm == "direct") {
-        degparms <- fit_direct$par[degparms_index]
-        errparms <- fit_direct$par[errparms_index]
-      } else {
+        control = control, ...))
+      if (!inherits(fit_direct, "try-error")) {
+        fit_direct$logLik <- - cost.current
         cost.current <- Inf # reset to avoid conflict with the OLS step
         data_direct <- current_data # We need this later if it was better
+        direct_failed = FALSE
+      } else {
+        direct_failed = TRUE
       }
     }
     if (error_model_algorithm != "direct") {
@@ -775,24 +775,30 @@ mkinfit <- function(mkinmod, observed,
 
       if (error_model_algorithm == "d_3") {
         d_3_messages = c(
+           direct_failed = "Direct fitting failed, results of three-step fitting are returned",
            same = "Direct fitting and three-step fitting yield approximately the same likelihood",
            threestep = "Three-step fitting yielded a higher likelihood than direct fitting",
            direct = "Direct fitting yielded a higher likelihood than three-step fitting")
-        rel_diff <- abs((fit_direct$logLik - fit$logLik))/-mean(c(fit_direct$logLik, fit$logLik))
-        if (rel_diff < 0.0001) {
-          if (!quiet) message(d_3_messages["same"])
-          fit$d_3_message <- d_3_messages["same"]
+        if (direct_failed) {
+          if (!quiet) message(d_3_messages["direct_failed"])
+          fit$d_3_message <- d_3_messages["direct_failed"]
         } else {
-          if (fit$logLik > fit_direct$logLik) {
-            if (!quiet) message(d_3_messages["threestep"])
-            fit$d_3_message <- d_3_messages["threestep"]
+          rel_diff <- abs((fit_direct$logLik - fit$logLik))/-mean(c(fit_direct$logLik, fit$logLik))
+          if (rel_diff < 0.0001) {
+            if (!quiet) message(d_3_messages["same"])
+            fit$d_3_message <- d_3_messages["same"]
           } else {
-            if (!quiet) message(d_3_messages["direct"])
-            fit <- fit_direct
-            fit$d_3_message <- d_3_messages["direct"]
-            degparms <- fit$par[degparms_index]
-            errparms <- fit$par[errparms_index]
-            current_data  <- data_direct
+            if (fit$logLik > fit_direct$logLik) {
+              if (!quiet) message(d_3_messages["threestep"])
+              fit$d_3_message <- d_3_messages["threestep"]
+            } else {
+              if (!quiet) message(d_3_messages["direct"])
+              fit <- fit_direct
+              fit$d_3_message <- d_3_messages["direct"]
+              degparms <- fit$par[degparms_index]
+              errparms <- fit$par[errparms_index]
+              current_data  <- data_direct
+            }
           }
         }
       }
diff --git a/R/mmkin.R b/R/mmkin.R
index adb9f4d0..d879fec4 100644
--- a/R/mmkin.R
+++ b/R/mmkin.R
@@ -19,8 +19,9 @@
 #' @param \dots Further arguments that will be passed to \code{\link{mkinfit}}.
 #' @importFrom parallel mclapply parLapply detectCores
 #' @return A two-dimensional \code{\link{array}} of \code{\link{mkinfit}}
-#'   objects that can be indexed using the model names for the first index (row index) 
-#'   and the dataset names for the second index (column index).
+#'   objects and/or try-errors that can be indexed using the model names for the
+#'   first index (row index) and the dataset names for the second index (column
+#'   index).
 #' @author Johannes Ranke
 #' @seealso \code{\link{[.mmkin}} for subsetting, \code{\link{plot.mmkin}} for
 #'   plotting.
diff --git a/check.log b/check.log
index cae31a24..b9a17217 100644
--- a/check.log
+++ b/check.log
@@ -1,5 +1,5 @@
 * using log directory ‘/home/jranke/git/mkin/mkin.Rcheck’
-* using R version 4.0.0 (2020-04-24)
+* using R version 4.0.2 (2020-06-22)
 * using platform: x86_64-pc-linux-gnu (64-bit)
 * using session charset: UTF-8
 * using options ‘--no-tests --as-cran’
@@ -7,8 +7,10 @@
 * checking extension type ... Package
 * this is package ‘mkin’ version ‘0.9.50.3’
 * package encoding: UTF-8
-* checking CRAN incoming feasibility ... Note_to_CRAN_maintainers
+* checking CRAN incoming feasibility ... NOTE
 Maintainer: ‘Johannes Ranke <jranke@uni-bremen.de>’
+
+The Date field is over a month old.
 * checking package namespace information ... OK
 * checking package dependencies ... OK
 * checking if this is a source package ... OK
@@ -67,5 +69,9 @@ Maintainer: ‘Johannes Ranke <jranke@uni-bremen.de>’
 * checking for detritus in the temp directory ... OK
 * DONE
 
-Status: OK
+Status: 1 NOTE
+See
+  ‘/home/jranke/git/mkin/mkin.Rcheck/00check.log’
+for details.
+
 
diff --git a/test.log b/test.log
index 3f217296..e96960ef 100644
--- a/test.log
+++ b/test.log
@@ -9,7 +9,7 @@ Testing mkin
 ✔ |   5       | Analytical solutions for coupled models [3.2 s]
 ✔ |   5       | Calculation of Akaike weights
 ✔ |  10       | Confidence intervals and p-values [1.1 s]
-✔ |  14       | Error model fitting [3.8 s]
+✔ |  14       | Error model fitting [3.7 s]
 ✔ |   4       | Test fitting the decline of metabolites from their maximum [0.3 s]
 ✔ |   1       | Fitting the logistic model [0.2 s]
 ✔ |   1       | Test dataset class mkinds used in gmkin
@@ -24,22 +24,22 @@ Reason: getRversion() < "4.1.0"  is TRUE
 test_nafta.R:53: skip: Test data from Appendix D are correctly evaluated
 Reason: getRversion() < "4.1.0"  is TRUE
 ────────────────────────────────────────────────────────────────────────────────
-✔ |   9       | Nonlinear mixed-effects models [7.8 s]
+✔ |   9       | Nonlinear mixed-effects models [7.7 s]
 ✔ |   0     1 | Plotting [0.7 s]
 ────────────────────────────────────────────────────────────────────────────────
 test_plot.R:24: skip: Plotting mkinfit and mmkin objects is reproducible
 Reason: getRversion() < "4.1.0"  is TRUE
 ────────────────────────────────────────────────────────────────────────────────
 ✔ |   4       | Residuals extracted from mkinfit models
-✔ |   2       | Complex test case from Schaefer et al. (2007) Piacenza paper [1.5 s]
+✔ |   2       | Complex test case from Schaefer et al. (2007) Piacenza paper [1.4 s]
 ✔ |   4       | Summary [0.1 s]
 ✔ |   1       | Summaries of old mkinfit objects
 ✔ |   4       | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.2 s]
 ✔ |   9       | Hypothesis tests [6.8 s]
-✔ |   4       | Calculation of maximum time weighted average concentrations (TWAs) [2.5 s]
+✔ |   4       | Calculation of maximum time weighted average concentrations (TWAs) [2.4 s]
 
 ══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 37.2 s
+Duration: 36.8 s
 
 OK:       145
 Failed:   0
diff --git a/tests/testthat/FOCUS_2006_D.csf b/tests/testthat/FOCUS_2006_D.csf
index aa9fc233..232d7dad 100644
--- a/tests/testthat/FOCUS_2006_D.csf
+++ b/tests/testthat/FOCUS_2006_D.csf
@@ -5,7 +5,7 @@ Description:
 MeasurementUnits: % AR
 TimeUnits: days
 Comments: Created using mkin::CAKE_export
-Date: 2020-06-15
+Date: 2020-07-15
 Optimiser: IRLS
 
 [Data]