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.

6 files changed, 69 insertions, 11 deletions

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 <http://www.gnu.org/licenses/>
+
+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")
+
+})