aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--R/create_deg_func.R32
-rw-r--r--R/mkinerrmin.R28
-rw-r--r--R/mkinfit.R75
-rw-r--r--R/mkinmod.R2
-rw-r--r--R/mkinpredict.R17
-rw-r--r--R/plot.mkinfit.R1
-rw-r--r--man/mkinfit.Rd6
-rw-r--r--man/mkinmod.Rd2
-rw-r--r--man/mkinpredict.Rd35
-rw-r--r--man/plot.mkinfit.Rd1
-rw-r--r--test.log33
-rw-r--r--tests/testthat/FOCUS_2006_D.csf2
-rw-r--r--tests/testthat/test_analytical.R11
-rw-r--r--tests/testthat/test_tests.R2
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)

Contact - Imprint