From 0b36fb666b84321814eabbc1a25687ee70d789e6 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
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.
---
 R/mkinfit.R | 20 ++++++++++++++------
 1 file changed, 14 insertions(+), 6 deletions(-)

(limited to 'R')

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)
     }
   }
 
-- 
cgit v1.2.1