From 9f8e1eb33b586beb7e889212bdababa081b6ff67 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 15 Jul 2020 12:30:39 +0200 Subject: Improve handling of (partially) failing fits --- NEWS.md | 4 ++++ R/mkinfit.R | 46 +++++++++++++++++++++++------------------ R/mmkin.R | 5 +++-- check.log | 12 ++++++++--- test.log | 10 ++++----- tests/testthat/FOCUS_2006_D.csf | 2 +- 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 ’ + +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 ’ * 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] -- cgit v1.2.1