From 0b36fb666b84321814eabbc1a25687ee70d789e6 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Thu, 21 Apr 2016 17:34:38 +0200 Subject: Avoid unnecessary stop in edge case, improve warnings mkinfit: Do not error out in cases where the fit converges, but the Jacobian for the untransformed model cost can not be estimated. Give a warning instead and return NA for the t-test results. summary.mkinfit: Give a warning message when the covariance matrix can not be obtained. A test has been added to containing such an edge case to check that the warnings are correctly issued and the fit does not terminate. --- DESCRIPTION | 2 +- NEWS.md | 12 ++++++++-- R/mkinfit.R | 20 +++++++++++----- build.log | 2 +- test.log | 2 +- tests/testthat/test_FOMC_ill-defined.R | 42 ++++++++++++++++++++++++++++++++++ 6 files changed, 69 insertions(+), 11 deletions(-) create mode 100644 tests/testthat/test_FOMC_ill-defined.R diff --git a/DESCRIPTION b/DESCRIPTION index 6d98cd32..2b027c58 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -3,7 +3,7 @@ Type: Package Title: Routines for Fitting Kinetic Models with One or More State Variables to Chemical Degradation Data Version: 0.9.42.9000 -Date: 2016-03-25 +Date: 2016-04-21 Authors@R: c(person("Johannes", "Ranke", role = c("aut", "cre", "cph"), email = "jranke@uni-bremen.de"), person("Katrin", "Lindenberger", role = "ctb"), diff --git a/NEWS.md b/NEWS.md index c40aba39..4c3bf690 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,11 +2,19 @@ ## mkin 0.9.42-9000 +### Major changes + +- `mkinfit`: Do not error out in cases where the fit converges, but the Jacobian for the untransformed model cost can not be estimated. Give a warning instead and return NA for the t-test results. + +- `summary.mkinfit`: Give a warning message when the covariance matrix can not be obtained. + +- A test has been added to containing a corresponding edge case to check that the warnings are correctly issued and the fit does not terminate. + ### Minor changes -- Remove an outdated reference to the inline package in the compiled_models vignette +- Remove an outdated reference to the inline package in the `compiled_models` vignette -## mkin 0.9.42 +## mkin 0.9.42 (2016-03-25) ### Major changes diff --git a/R/mkinfit.R b/R/mkinfit.R index 929c73f0..2302544b 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -483,9 +483,15 @@ mkinfit <- function(mkinmod, observed, # Estimate the Hessian for the model cost without parameter transformations # to make it possible to obtain the usual t-test # Code ported from FME::modFit - Jac_notrans <- gradient(function(p, ...) cost_notrans(p)$residuals$res, - bparms.optim, centered = TRUE) - fit$hessian_notrans <- 2 * t(Jac_notrans) %*% Jac_notrans + Jac_notrans <- try(gradient(function(p, ...) cost_notrans(p)$residuals$res, + bparms.optim, centered = TRUE), silent = TRUE) + if (inherits(Jac_notrans, "try-error")) { + warning("Calculation of the Jacobian failed for the cost function of the untransformed model.\n", + "No t-test results will be available") + fit$hessian_notrans <- NA + } else { + fit$hessian_notrans <- 2 * t(Jac_notrans) %*% Jac_notrans + } # Collect initial parameter values in three dataframes fit$start <- data.frame(value = c(state.ini.optim, @@ -548,7 +554,7 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, covar_notrans <- try(solve(0.5*object$hessian_notrans), silent = TRUE) # unscaled covariance rdf <- object$df.residual resvar <- object$ssr / rdf - if (!is.numeric(covar)) { + if (!is.numeric(covar) | is.na(covar[1])) { covar <- NULL se <- lci <- uci <- rep(NA, p) } else { @@ -558,7 +564,7 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, uci <- param + qt(1-alpha/2, rdf) * se } - if (!is.numeric(covar_notrans)) { + if (!is.numeric(covar_notrans) | is.na(covar_notrans[1])) { covar_notrans <- NULL se_notrans <- tval <- pval <- rep(NA, p) } else { @@ -690,7 +696,9 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), . rownames(Corr) <- colnames(Corr) <- rownames(x$par) print(Corr, digits = digits, ...) } else { - cat("Could not estimate covariance matrix; singular system:\n") + msg <- "Could not estimate covariance matrix; singular system:\n" + warning(msg) + cat(msg) } } diff --git a/build.log b/build.log index a5b27817..c497bbbe 100644 --- a/build.log +++ b/build.log @@ -6,5 +6,5 @@ * checking for LF line-endings in source and make files * checking for empty or unneeded directories * looking to see if a ‘data/datalist’ file should be added -* building ‘mkin_0.9.42.tar.gz’ +* building ‘mkin_0.9.42.9000.tar.gz’ diff --git a/test.log b/test.log index 90389b07..9079f341 100644 --- a/test.log +++ b/test.log @@ -4,4 +4,4 @@ Lade nötiges Paket: rootSolve Lade nötiges Paket: inline Lade nötiges Paket: parallel testthat results ================================================================ -OK: 40 SKIPPED: 0 FAILED: 0 +OK: 42 SKIPPED: 0 FAILED: 0 diff --git a/tests/testthat/test_FOMC_ill-defined.R b/tests/testthat/test_FOMC_ill-defined.R new file mode 100644 index 00000000..b036a0af --- /dev/null +++ b/tests/testthat/test_FOMC_ill-defined.R @@ -0,0 +1,42 @@ +# Copyright (C) 2016 Johannes Ranke +# Contact: jranke@uni-bremen.de + +# This file is part of the R package mkin + +# mkin is free software: you can redistribute it and/or modify it under the +# terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. + +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. + +# You should have received a copy of the GNU General Public License along with +# this program. If not, see + +context("Fitting the FOMC model with large parameter correlation") + +# Dataset that I ran across during my work and for which the calculation of the +# Jacobian failed. Data were slightly fuzzed. +FOMC_test <- data.frame( + name = "test_compound", + time = c(0, 14, 31, 59, 91), + value = c(45.8, 28.0, 28.5, 35.1, 35.6)) + +test_that("Fitting with large parameter correlation gives warnings", { + + # When fitting from the maximum, the Port algorithm does not converge (with + # default settings) + expect_warning( + fit.FOMC.Port <- mkinfit("FOMC", FOMC_test, method.modFit = "Port"), + "Optimisation by method Port did not converge") + + # When we use Levenberg-Marquardt, we get a problem estimating the Jacobian + # for the untransformed model + expect_warning( + fit.FOMC.Marq <- mkinfit("FOMC", FOMC_test, method.modFit = "Marq"), + "Calculation of the Jacobian failed for the cost function of the untransformed model") + +}) -- cgit v1.2.1