Improve handling of warnings, reorganize tests

3 files changed, 111 insertions, 14 deletions

diff --git a/R/mkinfit.R b/R/mkinfit.R
index d7b1b7f4..ec2d3412 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -233,7 +233,7 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "time", "value"))
 #' fit.FOMC = mkinfit("FOMC", FOCUS_2006_D, quiet = TRUE, error_model = "tc")
 #' fit.FOMC_SFO <- mkinfit(FOMC_SFO, FOCUS_D, quiet = TRUE,
 #'   parms.ini = fit.FOMC$bparms.ode, error_model = "tc")
-#'
+#' }
 #' @export
 mkinfit <- function(mkinmod, observed,
   parms.ini = "auto",
@@ -258,6 +258,8 @@ mkinfit <- function(mkinmod, observed,
 {
   call <- match.call()
 
+  summary_warnings <- character()
+
   # Derive the name used for the model
   if (is.character(mkinmod)) mkinmod_name <- mkinmod
   else mkinmod_name <- deparse(substitute(mkinmod))
@@ -289,7 +291,9 @@ mkinfit <- function(mkinmod, observed,
 
   # Also remove zero values to avoid instabilities (e.g. of the 'tc' error model)
   if (any(observed$value == 0)) {
-    warning("Observations with value of zero were removed from the data")
+    zero_warning <- "Observations with value of zero were removed from the data"
+    summary_warnings <- c(summary_warnings, zero_warning)
+    warning(zero_warning)
     observed <- subset(observed, value != 0)
   }
 
@@ -848,8 +852,9 @@ mkinfit <- function(mkinmod, observed,
   fit$error_model_algorithm <- error_model_algorithm
 
   if (fit$convergence != 0) {
-    fit$warning = paste0("Optimisation did not converge:\n", fit$message)
-    warning(fit$warning)
+    convergence_warning = paste0("Optimisation did not converge:\n", fit$message)
+    summary_warnings <- c(warnings, convergence_warning)
+    warning(convergence_warning)
   } else {
     if(!quiet) message("Optimisation successfully terminated.\n")
   }
@@ -918,14 +923,22 @@ mkinfit <- function(mkinmod, observed,
   fit$errparms <- errparms
   fit$df.residual <- n_observed - length(c(degparms, errparms))
 
+  # Assign the class here so method dispatch works for residuals
+  class(fit) <- c("mkinfit")
+
   # Check for normal distribution of residuals
-  fit$shapiro.p <- shapiro.test(residuals.mkinfit(fit, standardized = TRUE))$p.value
-  if (fit$shapiro.p < 0.05) warning("The p-value for the Shapiro-Wilk test of normality on standardized residuals is < 0.05")
+  fit$shapiro.p <- shapiro.test(residuals(fit, standardized = TRUE))$p.value
+  if (fit$shapiro.p < 0.05) {
+    shapiro_warning <- paste("Shapiro-Wilk test for standardized residuals: p = ", signif(fit$shapiro.p, 3))
+    warning(shapiro_warning)
+    summary_warnings <- c(summary_warnings, shapiro_warning)
+  }
+
+  fit$summary_warnings <- summary_warnings
 
   fit$date <- date()
   fit$version <- as.character(utils::packageVersion("mkin"))
   fit$Rversion <- paste(R.version$major, R.version$minor, sep=".")
 
-  class(fit) <- c("mkinfit")
   return(fit)
 }
diff --git a/R/summary.mkinfit.R b/R/summary.mkinfit.R
index 2c291ffd..f9858c32 100644
--- a/R/summary.mkinfit.R
+++ b/R/summary.mkinfit.R
@@ -122,7 +122,7 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05,
     date.fit = object$date,
     date.summary = date(),
     solution_type = object$solution_type,
-    warning = object$warning,
+    warnings = object$summary_warnings,
     use_of_ff = object$mkinmod$use_of_ff,
     error_model_algorithm = object$error_model_algorithm,
     df = c(p, rdf),
@@ -190,8 +190,6 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), .
   cat("Date of fit:    ", x$date.fit, "\n")
   cat("Date of summary:", x$date.summary, "\n")
 
-  if (!is.null(x$warning)) cat("\n\nWarning:", x$warning, "\n\n")
-
   cat("\nEquations:\n")
   nice_diffs <- gsub("^(d.*) =", "\\1/dt =", x[["diffs"]])
   writeLines(strwrap(nice_diffs, exdent = 11))
@@ -223,16 +221,20 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), .
   if(length(x$fixed$value) == 0) cat("None\n")
   else print(x$fixed)
 
+  # We used to only have one warning - show this for summarising old objects
+   if (!is.null(x[["warning"]])) cat("\n\nWarning:", x$warning, "\n\n")
+
+  if (length(x$warnings) > 0) {
+    cat("\n\nWarning(s):", "\n")
+    cat(x$warnings, sep = "\n")
+  }
+
   if (!is.null(x$AIC)) {
     cat("\nResults:\n\n")
     print(data.frame(AIC = x$AIC, BIC = x$BIC, logLik = x$logLik,
       row.names = " "))
   }
 
-  if (!is.null(x$shapiro.p) && x$shapiro.p < 0.05) {
-    warning("The p-value for the Shapiro-Wilk test of normality on standardized residuals is < 0.05")
-  }
-
   cat("\nOptimised, transformed parameters with symmetric confidence intervals:\n")
   print(signif(x$par, digits = digits))
 
diff --git a/R/summary_DFOP_FOCUS_C.txt b/R/summary_DFOP_FOCUS_C.txt
new file mode 100644
index 00000000..ab64a588
--- /dev/null
+++ b/R/summary_DFOP_FOCUS_C.txt
@@ -0,0 +1,82 @@
+mkin version used for fitting:    Dummy 0.0 for testing 
+R version used for fitting:       Dummy R version for testing 
+Date of fit:     Dummy date for testing 
+Date of summary: Dummy date for testing 
+
+Equations:
+d_parent/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
+           time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time)))
+           * parent
+
+Model predictions using solution type analytical 
+
+Fitted using test 0 model solutions performed in test time 0 s
+
+Error model: Constant variance 
+
+Error model algorithm: OLS 
+
+Starting values for parameters to be optimised:
+         value   type
+parent_0 85.10  state
+k1        0.10 deparm
+k2        0.01 deparm
+g         0.50 deparm
+
+Starting values for the transformed parameters actually optimised:
+             value lower upper
+parent_0 85.100000  -Inf   Inf
+log_k1   -2.302585  -Inf   Inf
+log_k2   -4.605170  -Inf   Inf
+g_ilr     0.000000  -Inf   Inf
+
+Fixed parameter values:
+None
+
+Results:
+
+       AIC      BIC    logLik
+  29.02372 30.00984 -9.511861
+
+Optimised, transformed parameters with symmetric confidence intervals:
+         Estimate Std. Error   Lower   Upper
+parent_0  85.0000    0.66620 83.1500 86.8500
+log_k1    -0.7775    0.03380 -0.8713 -0.6836
+log_k2    -4.0260    0.13100 -4.3890 -3.6620
+g_ilr      1.2490    0.05811  1.0870  1.4100
+sigma      0.6962    0.16410  0.2406  1.1520
+
+Parameter correlation:
+[1] "Correlation matrix is platform dependent, not tested"
+
+Backtransformed parameters:
+Confidence intervals for internally transformed parameters are asymmetric.
+t-test (unrealistically) based on the assumption of normal distribution
+for estimators of untransformed parameters.
+         Estimate t value    Pr(>t)    Lower    Upper
+parent_0 85.00000 127.600 1.131e-08 83.15000 86.85000
+k1        0.45960  29.580 3.887e-06  0.41840  0.50480
+k2        0.01785   7.636 7.901e-04  0.01241  0.02568
+g         0.85390  83.310 6.221e-08  0.82310  0.88020
+sigma     0.69620   4.243 6.618e-03  0.24060  1.15200
+
+FOCUS Chi2 error levels in percent:
+         err.min n.optim df
+All data   2.661       4  5
+parent     2.661       4  5
+
+Estimated disappearance times:
+        DT50  DT90 DT50_k1 DT50_k2
+parent 1.887 21.25   1.508   38.83
+
+Data:
+ time variable observed predicted residual
+    0   parent     85.1    85.003  0.09726
+    1   parent     57.9    58.039 -0.13912
+    3   parent     29.9    30.054 -0.15351
+    7   parent     14.6    13.866  0.73388
+   14   parent      9.7     9.787 -0.08657
+   28   parent      6.6     7.532 -0.93205
+   63   parent      4.0     4.033 -0.03269
+   91   parent      3.9     2.447  1.45348
+  119   parent      0.6     1.484 -0.88424