Improve handling of warnings, reorganize tests

14 files changed, 354 insertions, 166 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
diff --git a/build.log b/build.log
index b94e6450..0ed7b19a 100644
--- a/build.log
+++ b/build.log
@@ -2,8 +2,24 @@
 * preparing ‘mkin’:
 * checking DESCRIPTION meta-information ... OK
 * installing the package to build vignettes
-* creating vignettes ... OK
-* checking for LF line-endings in source and make files and shell scripts
-* checking for empty or unneeded directories
-* building ‘mkin_0.9.50.3.tar.gz’
+* creating vignettes ... ERROR
+--- re-building ‘FOCUS_D.rmd’ using rmarkdown
+--- finished re-building ‘FOCUS_D.rmd’
 
+--- re-building ‘FOCUS_L.rmd’ using rmarkdown
+Quitting from lines 80-83 (FOCUS_L.rmd) 
+Error: processing vignette 'FOCUS_L.rmd' failed with diagnostics:
+argument 1 (type 'list') cannot be handled by 'cat'
+--- failed re-building ‘FOCUS_L.rmd’
+
+--- re-building ‘mkin.rmd’ using rmarkdown
+--- finished re-building ‘mkin.rmd’
+
+--- re-building ‘twa.rmd’ using rmarkdown
+--- finished re-building ‘twa.rmd’
+
+SUMMARY: processing the following file failed:
+  ‘FOCUS_L.rmd’
+
+Error: Vignette re-building failed.
+Execution halted
diff --git a/man/mkinfit.Rd b/man/mkinfit.Rd
index 748dcb50..7f5092c5 100644
--- a/man/mkinfit.Rd
+++ b/man/mkinfit.Rd
@@ -214,6 +214,71 @@ When using the "IORE" submodel for metabolites, fitting with
 numerical ODE solver. In this situation it may help to switch off the
 internal rate transformation.
 }
+\examples{
+
+# Use shorthand notation for parent only degradation
+fit <- mkinfit("FOMC", FOCUS_2006_C, quiet = TRUE)
+summary(fit)
+
+# One parent compound, one metabolite, both single first order.
+# We remove zero values from FOCUS dataset D in order to avoid warnings
+FOCUS_D <- subset(FOCUS_2006_D, value != 0)
+# Use mkinsub for convenience in model formulation. Pathway to sink included per default.
+SFO_SFO <- mkinmod(
+  parent = mkinsub("SFO", "m1"),
+  m1 = mkinsub("SFO"))
+
+# Fit the model quietly to the FOCUS example dataset D using defaults
+fit <- mkinfit(SFO_SFO, FOCUS_D, quiet = TRUE)
+# Since mkin 0.9.50.3, we get a warning about non-normality of residuals,
+# so we try an alternative error model
+fit.tc <- mkinfit(SFO_SFO, FOCUS_D, quiet = TRUE, error_model = "tc")
+# This avoids the warning, and the likelihood ratio test confirms it is preferable
+lrtest(fit.tc, fit)
+# We can also allow for different variances of parent and metabolite as error model
+fit.obs <- mkinfit(SFO_SFO, FOCUS_D, quiet = TRUE, error_model = "obs")
+# This also avoids the warning about non-normality, but the two-component error model
+# has significantly higher likelihood
+lrtest(fit.obs, fit.tc)
+parms(fit.tc)
+endpoints(fit.tc)
+
+# We can show a quick (only one replication) benchmark for this case, as we
+# have several alternative solution methods for the model. We skip
+# uncompiled deSolve, as it is so slow. More benchmarks are found in the
+# benchmark vignette
+\dontrun{
+if(require(rbenchmark)) {
+  benchmark(replications = 1, order = "relative", columns = c("test", "relative", "elapsed"),
+    deSolve_compiled = mkinfit(SFO_SFO, FOCUS_D, quiet = TRUE, error_model = "tc",
+      solution_type = "deSolve", use_compiled = TRUE),
+    eigen = mkinfit(SFO_SFO, FOCUS_D, quiet = TRUE, error_model = "tc",
+      solution_type = "eigen"),
+    analytical = mkinfit(SFO_SFO, FOCUS_D, quiet = TRUE, error_model = "tc",
+      solution_type = "analytical"))
+}
+}
+
+# Use stepwise fitting, using optimised parameters from parent only fit, FOMC-SFO
+\dontrun{
+FOMC_SFO <- mkinmod(
+  parent = mkinsub("FOMC", "m1"),
+  m1 = mkinsub("SFO"))
+fit.FOMC_SFO <- mkinfit(FOMC_SFO, FOCUS_D, quiet = TRUE)
+# Again, we get a warning and try a more sophisticated error model
+fit.FOMC_SFO.tc <- mkinfit(FOMC_SFO, FOCUS_D, quiet = TRUE, error_model = "tc")
+# This model has a higher likelihood, but not significantly so
+lrtest(fit.tc, fit.FOMC_SFO.tc)
+# Also, the missing standard error for log_beta and the t-tests for alpha
+# and beta indicate overparameterisation
+summary(fit.FOMC_SFO.tc, data = FALSE)
+
+# We can easily use starting parameters from the parent only fit (only for illustration)
+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")
+}
+}
 \references{
 Rocke DM and Lorenzato S (1995) A two-component model
 for measurement error in analytical chemistry. \emph{Technometrics} 37(2), 176-184.
diff --git a/test.log b/test.log
index a1fa7dfb..c5f0e526 100644
--- a/test.log
+++ b/test.log
@@ -1,43 +1,36 @@
 Loading mkin
 Testing mkin
 ✔ |  OK F W S | Context
+✔ |   4       | AIC calculation
 ✔ |   2       | Export dataset for reading into CAKE
 ✔ |  14       | Results for FOCUS D established in expertise for UBA (Ranke 2014) [0.9 s]
 ✔ |   4       | Calculation of FOCUS chi2 error levels [0.4 s]
-✔ |   7       | Fitting the SFORB model [3.3 s]
+✔ |   7       | Fitting the SFORB model [3.4 s]
 ✔ |   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 [4.0 s]
+✔ |  10       | Confidence intervals and p-values [1.0 s]
+✔ |  14       | Error model fitting [4.1 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
 ✔ |  12       | Special cases of mkinfit calls [0.6 s]
-✔ |   8       | mkinmod model generation and printing [0.2 s]
+✔ |   8       | mkinmod model generation and printing [0.3 s]
 ✔ |   3       | Model predictions with mkinpredict [0.4 s]
 ✔ |  16       | Evaluations according to 2015 NAFTA guidance [1.5 s]
-✔ |   9       | Nonlinear mixed-effects models [7.8 s]
-✔ |   4       | Calculation of maximum time weighted average concentrations (TWAs) [2.4 s]
-✔ |   3   2   | Summary
-────────────────────────────────────────────────────────────────────────────────
-test_plots_summary_twa.R:59: warning: Summaries are reproducible
-The p-value for the Shapiro-Wilk test of normality on standardized residuals is < 0.05
-
-test_plots_summary_twa.R:72: warning: Summaries are reproducible
-The p-value for the Shapiro-Wilk test of normality on standardized residuals is < 0.05
-────────────────────────────────────────────────────────────────────────────────
+✔ |   9       | Nonlinear mixed-effects models [7.7 s]
 ✔ |  14       | Plotting [1.4 s]
-✔ |   4       | AIC calculation
 ✔ |   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.2 s]
 ✔ |   1       | Summaries of old mkinfit objects
 ✔ |   4       | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.1 s]
 ✔ |   9       | Hypothesis tests [6.6 s]
+✔ |   4       | Calculation of maximum time weighted average concentrations (TWAs) [2.4 s]
 
 ══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 38.3 s
+Duration: 38.2 s
 
-OK:       160
+OK:       161
 Failed:   0
-Warnings: 2
+Warnings: 0
 Skipped:  0
diff --git a/tests/testthat/summary_DFOP_FOCUS_D_deSolve.txt b/tests/testthat/summary_DFOP_FOCUS_D_deSolve.txt
index 9245c40b..0b31a81f 100644
--- a/tests/testthat/summary_DFOP_FOCUS_D_deSolve.txt
+++ b/tests/testthat/summary_DFOP_FOCUS_D_deSolve.txt
@@ -33,6 +33,10 @@ Fixed parameter values:
      value  type
 m1_0     0 state
 
+
+Warning(s): 
+Shapiro-Wilk test for standardized residuals: p =  0.0165
+
 Results:
 
        AIC      BIC    logLik
diff --git a/tests/testthat/summary_DFOP_FOCUS_D_eigen.txt b/tests/testthat/summary_DFOP_FOCUS_D_eigen.txt
index 141f57a1..8ced5c4d 100644
--- a/tests/testthat/summary_DFOP_FOCUS_D_eigen.txt
+++ b/tests/testthat/summary_DFOP_FOCUS_D_eigen.txt
@@ -33,6 +33,10 @@ Fixed parameter values:
      value  type
 m1_0     0 state
 
+
+Warning(s): 
+Shapiro-Wilk test for standardized residuals: p =  0.0165
+
 Results:
 
        AIC      BIC    logLik
diff --git a/tests/testthat/test_AIC.R b/tests/testthat/test_AIC.R
new file mode 100644
index 00000000..e9698f7c
--- /dev/null
+++ b/tests/testthat/test_AIC.R
@@ -0,0 +1,12 @@
+context("AIC calculation")
+
+test_that("The AIC is reproducible", {
+  expect_equivalent(AIC(fits[["SFO", "FOCUS_C"]]), 59.3, scale = 1, tolerance = 0.1)
+  expect_equivalent(AIC(fits[, "FOCUS_C"]),
+                    data.frame(df = c(3, 4, 5, 5), AIC = c(59.3, 44.7, 29.0, 39.2)),
+                    scale = 1, tolerance = 0.1)
+  expect_error(AIC(fits["SFO", ]), "column object")
+  expect_equivalent(BIC(fits[, "FOCUS_C"]),
+                    data.frame(df = c(3, 4, 5, 5), AIC = c(59.9, 45.5, 30.0, 40.2)),
+                    scale = 1, tolerance = 0.1)
+})
diff --git a/tests/testthat/test_FOCUS_D_UBA_expertise.R b/tests/testthat/test_FOCUS_D_UBA_expertise.R
index 101c8e15..0a7e5219 100644
--- a/tests/testthat/test_FOCUS_D_UBA_expertise.R
+++ b/tests/testthat/test_FOCUS_D_UBA_expertise.R
@@ -5,7 +5,7 @@ context("Results for FOCUS D established in expertise for UBA (Ranke 2014)")
 test_that("Fits without formation fractions are correct for FOCUS D", {
   expect_warning(
     fit.noff <- mkinfit(SFO_SFO, FOCUS_D, quiet = TRUE),
-    "p-value.*Shapiro-Wilk")
+    "Shapiro-Wilk")
 
   expect_equal(round(as.numeric(endpoints(fit.noff)$distimes["parent", ]), 2),
                c(7.02, 23.33))
diff --git a/tests/testthat/test_plot.R b/tests/testthat/test_plot.R
new file mode 100644
index 00000000..a33de07f
--- /dev/null
+++ b/tests/testthat/test_plot.R
@@ -0,0 +1,42 @@
+context("Plotting")
+
+test_that("Plotting mkinfit and mmkin objects is reproducible", {
+  skip_on_cran()
+  plot_default_FOCUS_C_SFO <- function() plot(fits[["SFO", "FOCUS_C"]])
+  plot_res_FOCUS_C_SFO <- function() plot(fits[["SFO", "FOCUS_C"]], show_residuals = TRUE)
+  plot_res_FOCUS_C_SFO_2 <- function() plot_res(fits[["SFO", "FOCUS_C"]])
+  plot_sep_FOCUS_C_SFO <- function() plot_sep(fits[["SFO", "FOCUS_C"]])
+  mkinparplot_FOCUS_C_SFO <- function() mkinparplot(fits[["SFO", "FOCUS_C"]])
+  mkinerrplot_FOCUS_C_SFO <- function() mkinerrplot(fits[["SFO", "FOCUS_C"]])
+  mmkin_FOCUS_C <- function() plot(fits[, "FOCUS_C"])
+  mmkin_SFO <- function() plot(fits["SFO",])
+  fit_D_obs_eigen <- suppressWarnings(mkinfit(SFO_SFO, FOCUS_2006_D, error_model = "obs", quiet = TRUE))
+  fit_C_tc <- mkinfit("SFO", FOCUS_2006_C, error_model = "tc", quiet = TRUE)
+
+  plot_errmod_fit_D_obs_eigen <- function() plot_err(fit_D_obs_eigen, sep_obs = FALSE)
+  plot_errmod_fit_C_tc <- function() plot_err(fit_C_tc)
+
+  plot_res_sfo_sfo <- function() plot_res(f_sfo_sfo_desolve)
+  plot_err_sfo_sfo <- function() plot_err(f_sfo_sfo_desolve)
+  plot_errmod_fit_obs_1 <- function() plot_err(fit_obs_1, sep_obs = FALSE)
+  plot_errmod_fit_tc_1 <- function() plot_err(fit_tc_1, sep_obs = FALSE)
+
+  skip_if(getRversion() > 4.0)
+  vdiffr::expect_doppelganger("mkinfit plot for FOCUS C with defaults", plot_default_FOCUS_C_SFO)
+  vdiffr::expect_doppelganger("mkinfit plot for FOCUS C with residuals like in gmkin", plot_res_FOCUS_C_SFO)
+  vdiffr::expect_doppelganger("plot_res for FOCUS C", plot_res_FOCUS_C_SFO_2)
+  vdiffr::expect_doppelganger("mkinfit plot for FOCUS C with sep = TRUE", plot_sep_FOCUS_C_SFO)
+  vdiffr::expect_doppelganger("mkinparplot for FOCUS C SFO", mkinparplot_FOCUS_C_SFO)
+  vdiffr::expect_doppelganger("mkinerrplot for FOCUS C SFO", mkinerrplot_FOCUS_C_SFO)
+  vdiffr::expect_doppelganger("mmkin plot for FOCUS C", mmkin_FOCUS_C)
+  vdiffr::expect_doppelganger("mmkin plot for SFO (FOCUS C and D)", mmkin_SFO)
+  vdiffr::expect_doppelganger("plot_errmod with FOCUS C tc", plot_errmod_fit_C_tc)
+  skip_on_travis() # Still not working on Travis, maybe because of deSolve producing
+  # different results when not working with a compiled model or eigenvalues
+  vdiffr::expect_doppelganger("plot_errmod with FOCUS D obs eigen", plot_errmod_fit_D_obs_eigen)
+  vdiffr::expect_doppelganger("plot_res for FOCUS D", plot_res_sfo_sfo)
+  vdiffr::expect_doppelganger("plot_err for FOCUS D", plot_err_sfo_sfo)
+  vdiffr::expect_doppelganger("plot_errmod with SFO_lin_a_tc", plot_errmod_fit_tc_1)
+  vdiffr::expect_doppelganger("plot_errmod with SFO_lin_a_obs", plot_errmod_fit_obs_1)
+})
+
diff --git a/tests/testthat/test_plots_summary_twa.R b/tests/testthat/test_plots_summary_twa.R
deleted file mode 100644
index aedc9da3..00000000
--- a/tests/testthat/test_plots_summary_twa.R
+++ /dev/null
@@ -1,128 +0,0 @@
-context("Calculation of maximum time weighted average concentrations (TWAs)")
-
-test_that("Time weighted average concentrations are correct", {
-  skip_on_cran()
-
-  outtimes_10 <- seq(0, 10, length.out = 10000)
-
-  ds <- "FOCUS_C"
-  for (model in models) {
-    fit <- fits[[model, ds]]
-    bpar <- summary(fit)$bpar[, "Estimate"]
-    pred_10 <- mkinpredict(fit$mkinmod,
-                odeparms = bpar[2:length(bpar)],
-                odeini = c(parent = bpar[[1]]),
-                outtimes = outtimes_10)
-    twa_num <- mean(pred_10[, "parent"])
-    names(twa_num) <- 10
-    twa_ana <- max_twa_parent(fit, 10)
-
-    # Test for absolute difference (scale = 1)
-    # The tolerance can be reduced if the length of outtimes is increased,
-    # but this needs more computing time so we stay with lenght.out = 10k
-    expect_equal(twa_num, twa_ana, tolerance = 0.003, scale = 1)
-  }
-})
-
-context("Summary")
-
-test_that("Summaries are reproducible", {
-  fit <- fits[["DFOP", "FOCUS_C"]]
-  test_summary <- summary(fit)
-  test_summary$fit_version <- "Dummy 0.0 for testing"
-  test_summary$fit_Rversion <- "Dummy R version for testing"
-  test_summary$date.fit <- "Dummy date for testing"
-  test_summary$date.summary <- "Dummy date for testing"
-  test_summary$calls <- "test 0"
-  test_summary$Corr <- signif(test_summary$Corr, 1)
-  test_summary$time <- c(elapsed = "test time 0")
-  # The correlation matrix is quite platform dependent
-  # It differs between i386 and amd64 on Windows
-  # and between Travis and my own Linux system
-  test_summary$Corr <- "Correlation matrix is platform dependent, not tested"
-  expect_known_output(print(test_summary), "summary_DFOP_FOCUS_C.txt")
-
-  test_summary_2 <- summary(f_sfo_sfo_eigen)
-  test_summary_2$fit_version <- "Dummy 0.0 for testing"
-  test_summary_2$fit_Rversion <- "Dummy R version for testing"
-  test_summary_2$date.fit <- "Dummy date for testing"
-  test_summary_2$date.summary <- "Dummy date for testing"
-  test_summary_2$calls <- "test 0"
-  test_summary_2$time <- c(elapsed = "test time 0")
-  # The correlation matrix is quite platform dependent
-  # It differs between i386 and amd64 on Windows
-  # and between Travis and my own Linux system
-  # Even more so when using the Eigen method
-  test_summary_2$Corr <- "Correlation matrix is platform dependent, not tested"
-  # The residuals for this method are also platform sensitive
-  test_summary_2$data$residual <- "not tested"
-  expect_known_output(print(test_summary_2), "summary_DFOP_FOCUS_D_eigen.txt")
-
-  test_summary_3 <- summary(f_sfo_sfo_desolve)
-  test_summary_3$fit_version <- "Dummy 0.0 for testing"
-  test_summary_3$fit_Rversion <- "Dummy R version for testing"
-  test_summary_3$date.fit <- "Dummy date for testing"
-  test_summary_3$date.summary <- "Dummy date for testing"
-  test_summary_3$calls <- "test 0"
-  test_summary_3$time <- c(elapsed = "test time 0")
-  # The correlation matrix is quite platform dependent
-  # It differs between i386 and amd64 on Windows
-  # and between Travis and my own Linux system
-  test_summary_3$Corr <- "Correlation matrix is platform dependent, not tested"
-  expect_known_output(print(test_summary_3), "summary_DFOP_FOCUS_D_deSolve.txt")
-})
-
-context("Plotting")
-
-test_that("Plotting mkinfit and mmkin objects is reproducible", {
-  skip_on_cran()
-  plot_default_FOCUS_C_SFO <- function() plot(fits[["SFO", "FOCUS_C"]])
-  plot_res_FOCUS_C_SFO <- function() plot(fits[["SFO", "FOCUS_C"]], show_residuals = TRUE)
-  plot_res_FOCUS_C_SFO_2 <- function() plot_res(fits[["SFO", "FOCUS_C"]])
-  plot_sep_FOCUS_C_SFO <- function() plot_sep(fits[["SFO", "FOCUS_C"]])
-  mkinparplot_FOCUS_C_SFO <- function() mkinparplot(fits[["SFO", "FOCUS_C"]])
-  mkinerrplot_FOCUS_C_SFO <- function() mkinerrplot(fits[["SFO", "FOCUS_C"]])
-  mmkin_FOCUS_C <- function() plot(fits[, "FOCUS_C"])
-  mmkin_SFO <- function() plot(fits["SFO",])
-  fit_D_obs_eigen <- suppressWarnings(mkinfit(SFO_SFO, FOCUS_2006_D, error_model = "obs", quiet = TRUE))
-  fit_C_tc <- mkinfit("SFO", FOCUS_2006_C, error_model = "tc", quiet = TRUE)
-
-  plot_errmod_fit_D_obs_eigen <- function() plot_err(fit_D_obs_eigen, sep_obs = FALSE)
-  plot_errmod_fit_C_tc <- function() plot_err(fit_C_tc)
-
-  plot_res_sfo_sfo <- function() plot_res(f_sfo_sfo_desolve)
-  plot_err_sfo_sfo <- function() plot_err(f_sfo_sfo_desolve)
-  plot_errmod_fit_obs_1 <- function() plot_err(fit_obs_1, sep_obs = FALSE)
-  plot_errmod_fit_tc_1 <- function() plot_err(fit_tc_1, sep_obs = FALSE)
-
-  skip_if(getRversion() > 4.0)
-  vdiffr::expect_doppelganger("mkinfit plot for FOCUS C with defaults", plot_default_FOCUS_C_SFO)
-  vdiffr::expect_doppelganger("mkinfit plot for FOCUS C with residuals like in gmkin", plot_res_FOCUS_C_SFO)
-  vdiffr::expect_doppelganger("plot_res for FOCUS C", plot_res_FOCUS_C_SFO_2)
-  vdiffr::expect_doppelganger("mkinfit plot for FOCUS C with sep = TRUE", plot_sep_FOCUS_C_SFO)
-  vdiffr::expect_doppelganger("mkinparplot for FOCUS C SFO", mkinparplot_FOCUS_C_SFO)
-  vdiffr::expect_doppelganger("mkinerrplot for FOCUS C SFO", mkinerrplot_FOCUS_C_SFO)
-  vdiffr::expect_doppelganger("mmkin plot for FOCUS C", mmkin_FOCUS_C)
-  vdiffr::expect_doppelganger("mmkin plot for SFO (FOCUS C and D)", mmkin_SFO)
-  vdiffr::expect_doppelganger("plot_errmod with FOCUS C tc", plot_errmod_fit_C_tc)
-  skip_on_travis() # Still not working on Travis, maybe because of deSolve producing
-  # different results when not working with a compiled model or eigenvalues
-  vdiffr::expect_doppelganger("plot_errmod with FOCUS D obs eigen", plot_errmod_fit_D_obs_eigen)
-  vdiffr::expect_doppelganger("plot_res for FOCUS D", plot_res_sfo_sfo)
-  vdiffr::expect_doppelganger("plot_err for FOCUS D", plot_err_sfo_sfo)
-  vdiffr::expect_doppelganger("plot_errmod with SFO_lin_a_tc", plot_errmod_fit_tc_1)
-  vdiffr::expect_doppelganger("plot_errmod with SFO_lin_a_obs", plot_errmod_fit_obs_1)
-})
-
-context("AIC calculation")
-
-test_that("The AIC is reproducible", {
-  expect_equivalent(AIC(fits[["SFO", "FOCUS_C"]]), 59.3, scale = 1, tolerance = 0.1)
-  expect_equivalent(AIC(fits[, "FOCUS_C"]),
-                    data.frame(df = c(3, 4, 5, 5), AIC = c(59.3, 44.7, 29.0, 39.2)),
-                    scale = 1, tolerance = 0.1)
-  expect_error(AIC(fits["SFO", ]), "column object")
-  expect_equivalent(BIC(fits[, "FOCUS_C"]),
-                    data.frame(df = c(3, 4, 5, 5), AIC = c(59.9, 45.5, 30.0, 40.2)),
-                    scale = 1, tolerance = 0.1)
-})
diff --git a/tests/testthat/test_summary.R b/tests/testthat/test_summary.R
new file mode 100644
index 00000000..5cf6ac6b
--- /dev/null
+++ b/tests/testthat/test_summary.R
@@ -0,0 +1,57 @@
+context("Summary")
+
+test_that("Summaries are reproducible", {
+  fit <- fits[["DFOP", "FOCUS_C"]]
+  test_summary <- summary(fit)
+  test_summary$fit_version <- "Dummy 0.0 for testing"
+  test_summary$fit_Rversion <- "Dummy R version for testing"
+  test_summary$date.fit <- "Dummy date for testing"
+  test_summary$date.summary <- "Dummy date for testing"
+  test_summary$calls <- "test 0"
+  test_summary$Corr <- signif(test_summary$Corr, 1)
+  test_summary$time <- c(elapsed = "test time 0")
+  # The correlation matrix is quite platform dependent
+  # It differs between i386 and amd64 on Windows
+  # and between Travis and my own Linux system
+  test_summary$Corr <- "Correlation matrix is platform dependent, not tested"
+  expect_known_output(print(test_summary), "summary_DFOP_FOCUS_C.txt")
+
+  test_summary_2 <- summary(f_sfo_sfo_eigen)
+  test_summary_2$fit_version <- "Dummy 0.0 for testing"
+  test_summary_2$fit_Rversion <- "Dummy R version for testing"
+  test_summary_2$date.fit <- "Dummy date for testing"
+  test_summary_2$date.summary <- "Dummy date for testing"
+  test_summary_2$calls <- "test 0"
+  test_summary_2$time <- c(elapsed = "test time 0")
+  # The correlation matrix is quite platform dependent
+  # It differs between i386 and amd64 on Windows
+  # and between Travis and my own Linux system
+  # Even more so when using the Eigen method
+  test_summary_2$Corr <- "Correlation matrix is platform dependent, not tested"
+  # The residuals for this method are also platform sensitive
+  test_summary_2$data$residual <- "not tested"
+  expect_known_output(print(test_summary_2), "summary_DFOP_FOCUS_D_eigen.txt")
+
+  test_summary_3 <- summary(f_sfo_sfo_desolve)
+  test_summary_3$fit_version <- "Dummy 0.0 for testing"
+  test_summary_3$fit_Rversion <- "Dummy R version for testing"
+  test_summary_3$date.fit <- "Dummy date for testing"
+  test_summary_3$date.summary <- "Dummy date for testing"
+  test_summary_3$calls <- "test 0"
+  test_summary_3$time <- c(elapsed = "test time 0")
+  # The correlation matrix is quite platform dependent
+  # It differs between i386 and amd64 on Windows
+  # and between Travis and my own Linux system
+  test_summary_3$Corr <- "Correlation matrix is platform dependent, not tested"
+  expect_known_output(print(test_summary_3), "summary_DFOP_FOCUS_D_deSolve.txt")
+})
+
+test_that("A fit generated with mkin 0.9.48.1 can be summarised", {
+  # Generated with mkin 0.9.48.1
+  # SFO_SFO <- mkinmod(parent = list(type = "SFO", to = "m1"),
+  #                    m1 = list(type = "SFO"), quiet = TRUE)
+  # fit_old <- mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE)
+  # save(fit_old, file = "~/git/mkin/inst/testdata/fit_old_FOCUS_D.rda", version = 2 )
+  load(system.file("testdata/fit_old_FOCUS_D.rda", package = "mkin"))
+  expect_true(length(summary(fit_old)) > 0)
+})
diff --git a/tests/testthat/test_twa.R b/tests/testthat/test_twa.R
new file mode 100644
index 00000000..a8fc0bd9
--- /dev/null
+++ b/tests/testthat/test_twa.R
@@ -0,0 +1,26 @@
+context("Calculation of maximum time weighted average concentrations (TWAs)")
+
+test_that("Time weighted average concentrations are correct", {
+  skip_on_cran()
+
+  outtimes_10 <- seq(0, 10, length.out = 10000)
+
+  ds <- "FOCUS_C"
+  for (model in models) {
+    fit <- fits[[model, ds]]
+    bpar <- summary(fit)$bpar[, "Estimate"]
+    pred_10 <- mkinpredict(fit$mkinmod,
+                odeparms = bpar[2:length(bpar)],
+                odeini = c(parent = bpar[[1]]),
+                outtimes = outtimes_10)
+    twa_num <- mean(pred_10[, "parent"])
+    names(twa_num) <- 10
+    twa_ana <- max_twa_parent(fit, 10)
+
+    # Test for absolute difference (scale = 1)
+    # The tolerance can be reduced if the length of outtimes is increased,
+    # but this needs more computing time so we stay with lenght.out = 10k
+    expect_equal(twa_num, twa_ana, tolerance = 0.003, scale = 1)
+  }
+})
+