Avoid the call to merge for analytical solutions

This increases performance up to a factor of five!

14 files changed, 135 insertions, 112 deletions

diff --git a/R/create_deg_func.R b/R/create_deg_func.R
index fb419a98..886c5e8b 100644
--- a/R/create_deg_func.R
+++ b/R/create_deg_func.R
@@ -21,19 +21,17 @@
 create_deg_func <- function(spec, use_of_ff = c("min", "max")) {
 
   use_of_ff <- match.arg(use_of_ff)
-
   min_ff <- use_of_ff == "min"
-
   obs_vars <- names(spec)
 
-  n <- character(0)
-
   parent <- obs_vars[1]
   parent_type <- spec[[1]]$type
 
   supported <- TRUE # This may be modified below
 
-  n[1] <- paste0(parent, " = ", parent_type, ".solution(outtimes, odeini['", parent,
+  predicted_text <- character(0)
+
+  predicted_text[parent] <- paste0(parent_type, ".solution(t, odeini['", parent,
     if (parent_type == "SFORB") "_free", "'], ",
     switch(parent_type,
       SFO = paste0("k_", parent, if (min_ff) "_sink" else "", ")"),
@@ -59,25 +57,27 @@ create_deg_func <- function(spec, use_of_ff = c("min", "max")) {
       k2 <- paste0("k_", n2)
       f12 <- paste0("f_", n1, "_to_", n2)
       if (parent_type == "SFO") {
-        n[2] <- paste0(n2, " = (((", k2, "-", k1, ")*", n20, "-", f12, "*", k1, "*", n10, ")*exp(-", k2, "*outtimes)+",
-          f12, "*", k1, "*", n10, "*exp(-", k1, "*outtimes))/(", k2, "-", k1, ")")
+        predicted_text[n2] <- paste0(
+          "(((", k2, "-", k1, ")*", n20, "-", f12, "*", k1, "*", n10, ")*exp(-", k2, "*t)+",
+          f12, "*", k1, "*", n10, "*exp(-", k1, "*t))/(", k2, "-", k1, ")")
       }
     }
   }
 
   if (supported) {
-    all_n <- paste(n, collapse = ",\n")
+    deg_func <- function(observed, odeini, odeparms) {}
 
     f_body <- paste0("{\n",
-      "out <- with(as.list(odeparms), {\n",
-      "data.frame(\n",
-        "time = outtimes,\n",
-        all_n, ")\n",
+      "predicted <- numeric(0)\n",
+      "with(as.list(odeparms), {\n")
+    for (obs_var in obs_vars) {
+      f_body <- paste0(f_body,
+      "t <- observed[observed$name == '", obs_var, "', 'time']\n",
+      "predicted <<- c(predicted, ", predicted_text[obs_var], ")\n")
+    }
+    f_body <- paste0(f_body,
       "})\n",
-      "return(out)\n}\n"
-    )
-
-    deg_func <- function(odeini, odeparms, outtimes) {}
+      "return(predicted)\n}\n")
 
     body(deg_func) <- parse(text = f_body)
     return(deg_func)
diff --git a/R/mkinerrmin.R b/R/mkinerrmin.R
index 0b647b81..1994aed0 100644
--- a/R/mkinerrmin.R
+++ b/R/mkinerrmin.R
@@ -1,12 +1,12 @@
 if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "value_mean"))
 
 #' Calculate the minimum error to assume in order to pass the variance test
-#' 
+#'
 #' This function finds the smallest relative error still resulting in passing
 #' the chi-squared test as defined in the FOCUS kinetics report from 2006.
-#' 
+#'
 #' This function is used internally by \code{\link{summary.mkinfit}}.
-#' 
+#'
 #' @param fit an object of class \code{\link{mkinfit}}.
 #' @param alpha The confidence level chosen for the chi-squared test.
 #' @importFrom stats qchisq aggregate
@@ -25,43 +25,41 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "value_mean"))
 #' \url{http://esdac.jrc.ec.europa.eu/projects/degradation-kinetics}
 #' @keywords manip
 #' @examples
-#' 
+#'
 #' SFO_SFO = mkinmod(parent = mkinsub("SFO", to = "m1"),
 #'                   m1 = mkinsub("SFO"),
 #'                   use_of_ff = "max")
-#' 
+#'
 #' fit_FOCUS_D = mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE)
 #' round(mkinerrmin(fit_FOCUS_D), 4)
 #' \dontrun{
 #'   fit_FOCUS_E = mkinfit(SFO_SFO, FOCUS_2006_E, quiet = TRUE)
 #'   round(mkinerrmin(fit_FOCUS_E), 4)
 #' }
-#' 
+#'
 #' @export
 mkinerrmin <- function(fit, alpha = 0.05)
 {
   parms.optim <- fit$par
 
   kinerrmin <- function(errdata, n.parms) {
-    means.mean <- mean(errdata$value_mean, na.rm = TRUE)
-    df = length(errdata$value_mean) - n.parms
+    means.mean <- mean(errdata$observed, na.rm = TRUE)
+    df = nrow(errdata) - n.parms
 
     err.min <- sqrt((1 / qchisq(1 - alpha, df)) *
-               sum((errdata$value_mean - errdata$value_pred)^2)/(means.mean^2))
+               sum((errdata$observed - errdata$predicted)^2)/(means.mean^2))
 
     return(list(err.min = err.min, n.optim = n.parms, df = df))
   }
 
-  means <- aggregate(value ~ time + name, data = fit$observed, mean, na.rm=TRUE)
-  errdata <- merge(means, fit$predicted, by = c("time", "name"),
-    suffixes = c("_mean", "_pred"))
-  errdata <- errdata[order(errdata$time, errdata$name), ]
+  errdata <- aggregate(cbind(observed, predicted) ~ time + variable, data = fit$data, mean, na.rm=TRUE)
+  errdata <- errdata[order(errdata$time, errdata$variable), ]
 
   # Remove values at time zero for variables whose value for state.ini is fixed,
   # as these will not have any effect in the optimization and should therefore not
   # be counted as degrees of freedom.
   fixed_initials = gsub("_0$", "", rownames(subset(fit$fixed, type = "state")))
-  errdata <- subset(errdata, !(time == 0 & name %in% fixed_initials))
+  errdata <- subset(errdata, !(time == 0 & variable %in% fixed_initials))
 
   n.optim.overall <- length(parms.optim) - length(fit$errparms)
 
@@ -73,7 +71,7 @@ mkinerrmin <- function(fit, alpha = 0.05)
   # The degrees of freedom are counted according to FOCUS kinetics (2011, p. 164)
   for (obs_var in fit$obs_vars)
   {
-    errdata.var <- subset(errdata, name == obs_var)
+    errdata.var <- subset(errdata, variable == obs_var)
 
     # Check if initial value is optimised
     n.initials.optim <- length(grep(paste(obs_var, ".*", "_0", sep=""), names(parms.optim)))
diff --git a/R/mkinfit.R b/R/mkinfit.R
index 6cfb88c0..c33afbf0 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -288,6 +288,10 @@ mkinfit <- function(mkinmod, observed,
     observed <- subset(observed, value != 0)
   }
 
+  # Sort observed values for efficient analytical predictions
+  observed$name <- ordered(observed$name, levels = obs_vars)
+  observed <- observed[order(observed$name, observed$time), ]
+
   # Obtain data for decline from maximum mean value if requested
   if (from_max_mean) {
     # This is only used for simple decline models
@@ -520,7 +524,7 @@ mkinfit <- function(mkinmod, observed,
   }
 
   # Unique outtimes for model solution.
-  outtimes = sort(unique(observed$time))
+  outtimes <- sort(unique(observed$time))
 
   # Define the objective function for optimisation, including (back)transformations
   cost_function <- function(P, trans = TRUE, OLS = FALSE, fixed_degparms = FALSE, fixed_errparms = FALSE, update_data = TRUE, ...)
@@ -578,47 +582,51 @@ mkinfit <- function(mkinmod, observed,
     }
 
     # Solve the system with current parameter values
-    out <- mkinpredict(mkinmod, parms,
-                       odeini, outtimes,
-                       solution_type = solution_type,
-                       use_compiled = use_compiled,
-                       method.ode = method.ode,
-                       atol = atol, rtol = rtol, ...)
-
-    out_long <- mkin_wide_to_long(out, time = "time")
-
-    # Surprisingly, the next line is the one taking the most time for one, two
-    # or three observed variables if we use compiled ODE models
-    # as evidenced by use of github:hadley/lineprof
-    cost_data <- merge(observed[c("name", "time", "value")], out_long,
-                         by = c("name", "time"), suffixes = c(".observed", ".predicted"))
+    if (solution_type == "analytical") {
+      observed$predicted <- mkinmod$deg_func(observed, odeini, parms)
+    } else {
+      out <- mkinpredict(mkinmod, parms,
+                         odeini, outtimes,
+                         solution_type = solution_type,
+                         use_compiled = use_compiled,
+                         method.ode = method.ode,
+                         atol = atol, rtol = rtol, ...)
+
+      out_long <- mkin_wide_to_long(out, time = "time")
+
+      out_merged <- merge(observed[c("name", "time")], out_long)
+      out_merged$name <- ordered(out_merged$name, levels = obs_vars)
+      out_merged <- out_merged[order(out_merged$name, out_merged$time), ]
+      observed$predicted <- out_merged$value
+    }
 
+    # Define standard deviation for each observation
     if (err_mod == "const") {
-      cost_data$std <- if (OLS) NA else cost_errparms["sigma"]
+      observed$std <- if (OLS) NA else cost_errparms["sigma"]
     }
     if (err_mod == "obs") {
-      std_names <- paste0("sigma_", cost_data$name)
-      cost_data$std <- cost_errparms[std_names]
+      std_names <- paste0("sigma_", observed$name)
+      observed$std <- cost_errparms[std_names]
     }
     if (err_mod == "tc") {
-      cost_data$std <- sqrt(cost_errparms["sigma_low"]^2 + cost_data$value.predicted^2 * cost_errparms["rsd_high"]^2)
+      observed$std <- sqrt(cost_errparms["sigma_low"]^2 + observed$predicted^2 * cost_errparms["rsd_high"]^2)
     }
 
+    # Calculate model cost
     if (OLS) {
       # Cost is the sum of squared residuals
-      cost <- with(cost_data, sum((value.observed - value.predicted)^2))
+      cost <- with(observed, sum((value - predicted)^2))
     } else {
       # Cost is the negative log-likelihood
-      cost <- - with(cost_data,
-        sum(dnorm(x = value.observed, mean = value.predicted, sd = std, log = TRUE)))
+      cost <- - with(observed,
+        sum(dnorm(x = value, mean = predicted, sd = std, log = TRUE)))
     }
 
     # We update the current cost and data during the optimisation, not
     # during hessian calculations
     if (update_data) {
 
-      assign("out_predicted", out_long, inherits = TRUE)
-      assign("current_data", cost_data, inherits = TRUE)
+      assign("current_data", observed, inherits = TRUE)
 
       if (cost < cost.current) {
         assign("cost.current", cost, inherits = TRUE)
@@ -718,7 +726,7 @@ mkinfit <- function(mkinmod, observed,
       degparms <- fit$par
 
       # Get the maximum likelihood estimate for sigma at the optimum parameter values
-      current_data$residual <- current_data$value.observed - current_data$value.predicted
+      current_data$residual <- current_data$value - current_data$predicted
       sigma_mle <- sqrt(sum(current_data$residual^2)/nrow(current_data))
 
       # Use that estimate for the constant variance, or as first guess if err_mod = "obs"
@@ -855,11 +863,7 @@ mkinfit <- function(mkinmod, observed,
   # We also need the model and a model name for summary and plotting
   fit$mkinmod <- mkinmod
   fit$mkinmod$name <- mkinmod_name
-
-  # We need data and predictions for summary and plotting
-  fit$observed <- observed
   fit$obs_vars <- obs_vars
-  fit$predicted <- out_predicted
 
   # Residual sum of squares as a function of the fitted parameters
   fit$rss <- function(P) cost_function(P, OLS = TRUE, update_data = FALSE)
@@ -886,15 +890,10 @@ mkinfit <- function(mkinmod, observed,
   fit$fixed$type = c(rep("state", length(state.ini.fixed)),
                      rep("deparm", length(parms.fixed)))
 
-  # Sort observed, predicted and residuals
-  current_data$name <- ordered(current_data$name, levels = obs_vars)
-
-  ordered_data <- current_data[order(current_data$name, current_data$time), ]
-
-  fit$data <- data.frame(time = ordered_data$time,
-                         variable = ordered_data$name,
-                         observed = ordered_data$value.observed,
-                         predicted = ordered_data$value.predicted)
+  fit$data <- data.frame(time = current_data$time,
+                         variable = current_data$name,
+                         observed = current_data$value,
+                         predicted = current_data$predicted)
 
   fit$data$residual <- fit$data$observed - fit$data$predicted
 
diff --git a/R/mkinmod.R b/R/mkinmod.R
index 21551861..7e9e9248 100644
--- a/R/mkinmod.R
+++ b/R/mkinmod.R
@@ -70,7 +70,7 @@
 #' @examples
 #'
 #' # Specify the SFO model (this is not needed any more, as we can now mkinfit("SFO", ...)
-#' SFO <- mkinmod(parent = list(type = "SFO"))
+#' SFO <- mkinmod(parent = mkinsub("SFO"))
 #'
 #' # One parent compound, one metabolite, both single first order
 #' SFO_SFO <- mkinmod(
diff --git a/R/mkinpredict.R b/R/mkinpredict.R
index 8c76ca05..75582fac 100644
--- a/R/mkinpredict.R
+++ b/R/mkinpredict.R
@@ -92,10 +92,15 @@
 #'       c(parent = 100, m1 = 0), seq(0, 20, by = 0.1),
 #'       solution_type = "analytical", use_compiled = FALSE)[201,])
 #' }
+#'     analytical = mkinpredict(SFO_SFO,
+#'       c(k_parent = 0.15, f_parent_to_m1 = 0.5, k_m1 = 0.01),
+#'       c(parent = 100, m1 = 0), seq(0, 20, by = 0.1),
+#'       solution_type = "analytical", use_compiled = FALSE)[201,]
 #'
 #' \dontrun{
 #'   # Predict from a fitted model
 #'   f <- mkinfit(SFO_SFO, FOCUS_2006_C, quiet = TRUE)
+#'   f <- mkinfit(SFO_SFO, FOCUS_2006_C, quiet = TRUE, solution_type = "analytical")
 #'   head(mkinpredict(f))
 #' }
 #'
@@ -131,8 +136,16 @@ mkinpredict.mkinmod <- function(x,
   }
 
   if (solution_type == "analytical") {
-    out <- x$deg_func(odeini = odeini,
-      odeparms = odeparms, outtimes = outtimes)
+    # This is clumsy, as we wanted fast analytical predictions for mkinfit
+    out <- data.frame(time = outtimes)
+    obs_vars <- names(x$spec)
+    pseudo_observed <-
+      data.frame(name = rep(obs_vars, each = length(outtimes)),
+      time = rep(outtimes, length(obs_vars)))
+    pseudo_observed$predicted <- x$deg_func(pseudo_observed, odeini, odeparms)
+    for (obs_var in obs_vars) {
+      out[obs_var] <- pseudo_observed[pseudo_observed$name == obs_var, "predicted"]
+    }
   }
 
   if (solution_type == "eigen") {
diff --git a/R/plot.mkinfit.R b/R/plot.mkinfit.R
index 4b5bfa3a..b2cb4890 100644
--- a/R/plot.mkinfit.R
+++ b/R/plot.mkinfit.R
@@ -69,6 +69,7 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("type", "variable", "obse
 #' \dontrun{
 #' SFO_SFO <- mkinmod(parent = mkinsub("SFO", "m1", full = "Parent"),
 #'                    m1 = mkinsub("SFO", full = "Metabolite M1" ))
+#' fit <- mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE)
 #' fit <- mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE, error_model = "tc")
 #' plot(fit)
 #' plot_res(fit)
diff --git a/man/mkinfit.Rd b/man/mkinfit.Rd
index b728aded..db3f5f3e 100644
--- a/man/mkinfit.Rd
+++ b/man/mkinfit.Rd
@@ -26,7 +26,6 @@ mkinfit(
   quiet = FALSE,
   atol = 1e-08,
   rtol = 1e-10,
-  n.outtimes = 10,
   error_model = c("const", "obs", "tc"),
   error_model_algorithm = c("auto", "d_3", "direct", "twostep", "threestep",
     "fourstep", "IRLS", "OLS"),
@@ -141,11 +140,6 @@ is 1e-8, lower than in \code{\link{lsoda}}.}
 \item{rtol}{Absolute error tolerance, passed to \code{\link{ode}}. Default
 is 1e-10, much lower than in \code{\link{lsoda}}.}
 
-\item{n.outtimes}{The length of the dataseries that is produced by the model
-prediction function \code{\link{mkinpredict}}. This impacts the accuracy
-of the numerical solver if that is used (see \code{solution_type}
-argument.}
-
 \item{error_model}{If the error model is "const", a constant standard
   deviation is assumed.
 
diff --git a/man/mkinmod.Rd b/man/mkinmod.Rd
index 2ba917d6..a5736be7 100644
--- a/man/mkinmod.Rd
+++ b/man/mkinmod.Rd
@@ -82,7 +82,7 @@ The IORE submodel is not well tested for metabolites. When using this
 \examples{
 
 # Specify the SFO model (this is not needed any more, as we can now mkinfit("SFO", ...)
-SFO <- mkinmod(parent = list(type = "SFO"))
+SFO <- mkinmod(parent = mkinsub("SFO"))
 
 # One parent compound, one metabolite, both single first order
 SFO_SFO <- mkinmod(
diff --git a/man/mkinpredict.Rd b/man/mkinpredict.Rd
index f7e4acfc..1ba5b5f8 100644
--- a/man/mkinpredict.Rd
+++ b/man/mkinpredict.Rd
@@ -91,7 +91,7 @@ FALSE).}
 solver is used.}
 }
 \value{
-A matrix in the same format as the output of \code{\link{ode}}.
+A data frame with the solution in wide format
 }
 \description{
 This function produces a time series for all the observed variables in a
@@ -131,30 +131,35 @@ mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100),
 
 # Check compiled model versions - they are faster than the eigenvalue based solutions!
 SFO_SFO = mkinmod(parent = list(type = "SFO", to = "m1"),
-                  m1 = list(type = "SFO"), use_of_ff = "min")
+                  m1 = list(type = "SFO"), use_of_ff = "max")
 if(require(rbenchmark)) {
-  benchmark(
-    eigen = mkinpredict(SFO_SFO, c(k_parent_m1 = 0.05, k_parent_sink = 0.1, k_m1_sink = 0.01),
+  benchmark(replications = 10, order = "relative", columns = c("test", "relative", "elapsed"),
+    eigen = mkinpredict(SFO_SFO,
+      c(k_parent = 0.15, f_parent_to_m1 = 0.5, k_m1 = 0.01),
       c(parent = 100, m1 = 0), seq(0, 20, by = 0.1),
       solution_type = "eigen")[201,],
-    deSolve_compiled = mkinpredict(SFO_SFO, 
-      c(k_parent_m1 = 0.05, k_parent_sink = 0.1, k_m1_sink = 0.01),
+    deSolve_compiled = mkinpredict(SFO_SFO,
+      c(k_parent = 0.15, f_parent_to_m1 = 0.5, k_m1 = 0.01),
       c(parent = 100, m1 = 0), seq(0, 20, by = 0.1),
       solution_type = "deSolve")[201,],
-    deSolve = mkinpredict(SFO_SFO, c(k_parent_m1 = 0.05, k_parent_sink = 0.1, k_m1_sink = 0.01),
-     c(parent = 100, m1 = 0), seq(0, 20, by = 0.1),
-     solution_type = "deSolve", use_compiled = FALSE)[201,],
-    replications = 10)
+    deSolve = mkinpredict(SFO_SFO,
+      c(k_parent = 0.15, f_parent_to_m1 = 0.5, k_m1 = 0.01),
+      c(parent = 100, m1 = 0), seq(0, 20, by = 0.1),
+      solution_type = "deSolve", use_compiled = FALSE)[201,],
+    analytical = mkinpredict(SFO_SFO,
+      c(k_parent = 0.15, f_parent_to_m1 = 0.5, k_m1 = 0.01),
+      c(parent = 100, m1 = 0), seq(0, 20, by = 0.1),
+      solution_type = "analytical", use_compiled = FALSE)[201,])
 }
-
-# Since mkin 0.9.49.11 we also have analytical solutions for some models, including SFO-SFO
-#    deSolve = mkinpredict(SFO_SFO, c(k_parent_m1 = 0.05, k_parent_sink = 0.1, k_m1_sink = 0.01),
-#     c(parent = 100, m1 = 0), seq(0, 20, by = 0.1),
-#     solution_type = "analytical", use_compiled = FALSE)[201,],
+    analytical = mkinpredict(SFO_SFO,
+      c(k_parent = 0.15, f_parent_to_m1 = 0.5, k_m1 = 0.01),
+      c(parent = 100, m1 = 0), seq(0, 20, by = 0.1),
+      solution_type = "analytical", use_compiled = FALSE)[201,]
 
 \dontrun{
   # Predict from a fitted model
   f <- mkinfit(SFO_SFO, FOCUS_2006_C, quiet = TRUE)
+  f <- mkinfit(SFO_SFO, FOCUS_2006_C, quiet = TRUE, solution_type = "analytical")
   head(mkinpredict(f))
 }
 
diff --git a/man/plot.mkinfit.Rd b/man/plot.mkinfit.Rd
index c3f3134a..4235557e 100644
--- a/man/plot.mkinfit.Rd
+++ b/man/plot.mkinfit.Rd
@@ -144,6 +144,7 @@ latex is being used for the formatting of the chi2 error level, if
 \dontrun{
 SFO_SFO <- mkinmod(parent = mkinsub("SFO", "m1", full = "Parent"),
                    m1 = mkinsub("SFO", full = "Metabolite M1" ))
+fit <- mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE)
 fit <- mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE, error_model = "tc")
 plot(fit)
 plot_res(fit)
diff --git a/test.log b/test.log
index 489550a2..7788f109 100644
--- a/test.log
+++ b/test.log
@@ -2,34 +2,35 @@ Loading mkin
 Testing mkin
 ✔ |  OK F W S | Context
 ✔ |   2       | Export dataset for reading into CAKE
-✔ |  13       | Results for FOCUS D established in expertise for UBA (Ranke 2014) [3.1 s]
-✔ |   4       | Calculation of FOCUS chi2 error levels [1.9 s]
-✔ |   7       | Fitting the SFORB model [9.6 s]
+✔ |  13       | Results for FOCUS D established in expertise for UBA (Ranke 2014) [3.3 s]
+✔ |   4       | Calculation of FOCUS chi2 error levels [2.0 s]
+✔ |   7       | Fitting the SFORB model [9.1 s]
+✔ |   1       | Analytical solutions for coupled models [1.1 s]
 ✔ |   5       | Calculation of Akaike weights
-✔ |  10       | Confidence intervals and p-values [8.3 s]
-✔ |  14       | Error model fitting [22.1 s]
-✔ |   6       | Test fitting the decline of metabolites from their maximum [0.7 s]
-✔ |   1       | Fitting the logistic model [0.7 s]
+✔ |  10       | Confidence intervals and p-values [1.0 s]
+✔ |  14       | Error model fitting [3.7 s]
+✔ |   6       | Test fitting the decline of metabolites from their maximum [0.2 s]
+✔ |   1       | Fitting the logistic model [0.2 s]
 ✔ |   1       | Test dataset class mkinds used in gmkin
-✔ |  12       | Special cases of mkinfit calls [2.1 s]
+✔ |  12       | Special cases of mkinfit calls [0.6 s]
 ✔ |   8       | mkinmod model generation and printing [0.2 s]
 ✔ |   3       | Model predictions with mkinpredict [0.4 s]
-✔ |  16       | Evaluations according to 2015 NAFTA guidance [3.8 s]
-✔ |   9       | Nonlinear mixed-effects models [11.9 s]
+✔ |  16       | Evaluations according to 2015 NAFTA guidance [1.4 s]
+✔ |   9       | Nonlinear mixed-effects models [12.1 s]
 ✔ |   4       | Calculation of maximum time weighted average concentrations (TWAs) [2.4 s]
 ✔ |   3       | Summary
-✔ |  14       | Plotting [4.0 s]
+✔ |  14       | Plotting [3.5 s]
 ✔ |   4       | AIC calculation
 ✔ |   4       | Residuals extracted from mkinfit models
-✔ |   2       | Complex test case from Schaefer et al. (2007) Piacenza paper [3.8 s]
+✔ |   2       | Complex test case from Schaefer et al. (2007) Piacenza paper [4.2 s]
 ✔ |   1       | Summaries of old mkinfit objects
-✔ |   4       | Results for synthetic data established in expertise for UBA (Ranke 2014) [6.0 s]
-✔ |   9       | Hypothesis tests [21.3 s]
+✔ |   4       | Results for synthetic data established in expertise for UBA (Ranke 2014) [6.4 s]
+✔ |   9       | Hypothesis tests [17.8 s]
 
 ══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 102.5 s
+Duration: 69.9 s
 
-OK:       156
+OK:       157
 Failed:   0
 Warnings: 0
 Skipped:  0
diff --git a/tests/testthat/FOCUS_2006_D.csf b/tests/testthat/FOCUS_2006_D.csf
index f27df430..f113a668 100644
--- a/tests/testthat/FOCUS_2006_D.csf
+++ b/tests/testthat/FOCUS_2006_D.csf
@@ -5,7 +5,7 @@ Description:
 MeasurementUnits: % AR
 TimeUnits: days
 Comments: Created using mkin::CAKE_export
-Date: 2020-05-08
+Date: 2020-05-09
 Optimiser: IRLS
 
 [Data]
diff --git a/tests/testthat/test_analytical.R b/tests/testthat/test_analytical.R
new file mode 100644
index 00000000..3d30e042
--- /dev/null
+++ b/tests/testthat/test_analytical.R
@@ -0,0 +1,11 @@
+context("Analytical solutions for coupled models")
+
+test_that("The analytical solution of SFO-SFO is correct", {
+  f_sfo_sfo.ff.analytical <- mkinfit(SFO_SFO.ff,
+    subset(FOCUS_2006_D, value != 0),
+    quiet = TRUE)
+  expect_equal(
+    parms(f_sfo_sfo.ff),
+    parms(f_sfo_sfo.ff.analytical)
+  )
+})
diff --git a/tests/testthat/test_tests.R b/tests/testthat/test_tests.R
index 3d1c136c..b7c3c548 100644
--- a/tests/testthat/test_tests.R
+++ b/tests/testthat/test_tests.R
@@ -49,7 +49,7 @@ test_that("Updating fitted models works", {
     error_model = "tc", quiet = TRUE)
 
   f_soil_1_nw <- update(f_soil_1_tc, error_model = "const")
-  f_soil_1_nw_A2 <- update(f_soil_1_nw, fixed_parms = c(k_A2 = 0))
+  f_soil_1_nw_A2 <- suppressWarnings(update(f_soil_1_nw, fixed_parms = c(k_A2 = 0)))
   test_nw_tc <- lrtest(f_soil_1_nw, f_soil_1_tc)
   expect_equivalent(test_nw_tc[["2", "Pr(>Chisq)"]], 2.113e-6)
   test_nw_A2 <- lrtest(f_soil_1_nw, f_soil_1_nw_A2)