aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--DESCRIPTION4
-rw-r--r--NAMESPACE3
-rw-r--r--NEWS.md6
-rw-r--r--R/mkinfit.R130
-rw-r--r--check.log2
-rw-r--r--test.log104
-rw-r--r--tests/testthat/test_irls.R92
-rw-r--r--tests/testthat/test_twa.R8
8 files changed, 244 insertions, 105 deletions
diff --git a/DESCRIPTION b/DESCRIPTION
index 53285036..6f22f986 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -2,7 +2,7 @@ Package: mkin
Type: Package
Title: Kinetic Evaluation of Chemical Degradation Data
Version: 0.9.47.6
-Date: 2018-09-18
+Date: 2018-09-21
Authors@R: c(person("Johannes", "Ranke", role = c("aut", "cre", "cph"),
email = "jranke@uni-bremen.de",
comment = c(ORCID = "0000-0003-4371-6538")),
@@ -18,7 +18,7 @@ Description: Calculation routines based on the FOCUS Kinetics Report (2006,
warranty is implied for correctness of results or fitness for a particular
purpose.
Imports: stats, graphics, methods, FME, deSolve, R6, minpack.lm, rootSolve,
- inline, parallel
+ inline, parallel, plyr
Suggests: knitr, rbenchmark, tikzDevice, testthat
License: GPL
LazyLoad: yes
diff --git a/NAMESPACE b/NAMESPACE
index cfbe2851..764b4aff 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,5 +1,5 @@
# Export all names
-exportPattern(".")
+exportPattern("^[^\\.]")
S3method(print, mkinds)
S3method(print, mkinmod)
S3method(plot, mkinfit)
@@ -23,3 +23,4 @@ importFrom(deSolve, ode)
importFrom(methods, signature)
importFrom(R6, R6Class)
importFrom(grDevices, dev.cur)
+importFrom(plyr, join)
diff --git a/NEWS.md b/NEWS.md
index 8372a8b5..cce1bea6 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,4 +1,4 @@
-# mkin 0.9.47.6 (2018-09-18)
+# mkin 0.9.47.6 (2018-09-21)
- 'add_err': Respect the argument giving the number of replicates in the synthetic dataset
@@ -6,6 +6,10 @@
- 'mkinpredict': Make the function generic and create a method for mkinfit objects
+- 'mkinfit': Improve the correctness of the fitted two component error model by fitting the mean absolute deviance at each observation against the observed values, weighting with the current two-component error model
+
+- 'tests/testthat/test_irls.R': Test if the components of the error model used to generate the data can be reproduced with moderate accuracy
+
# mkin 0.9.47.5 (2018-09-14)
- Make the two-component error model stop in cases where it is inadequate to avoid nls crashes on windows
diff --git a/R/mkinfit.R b/R/mkinfit.R
index d56c663a..8c7549ad 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -405,6 +405,7 @@ mkinfit <- function(mkinmod, observed,
if (! reweight.method %in% c("obs", "tc")) stop("Only reweighting methods 'obs' and 'tc' are implemented")
if (reweight.method == "obs") {
+ tc_fit <- NA
if(!quiet) {
cat("IRLS based on variance estimates for each observed variable\n")
cat("Initial variance estimates are:\n")
@@ -412,32 +413,20 @@ mkinfit <- function(mkinmod, observed,
}
}
if (reweight.method == "tc") {
- # We need unweighted residuals to update the weighting
- tmp_res <- cost(fit$par)$residuals
-
- mad_agg <- aggregate(tmp_res$res.unweighted,
- by = list(name = tmp_res$name, x_res = tmp_res$x),
- FUN = function(x) mad(x, center = 0))
- names(mad_agg) <- c("name", "x", "mad")
- tmp_res_mad <- merge(tmp_res, mad_agg)
-
- tc_fit <- try(
- nls(mad ~ sigma_twocomp(mod, sigma_low, rsd_high),
- start = list(sigma_low = tc["sigma_low"], rsd_high = tc["rsd_high"]),
- data = tmp_res_mad,
- lower = 0,
- algorithm = "port"))
-
- if (inherits(tc_fit, "try-error")) {
- stop("Estimation of the two error model components failed for the initial fit.\n",
- "Try without reweighting or with reweight.method = 'obs'.")
- }
+ tc_fit <- .fit_error_model_mad_obs(cost(fit$par)$residuals, tc, 0)
- tc_fitted <- coef(tc_fit)
- if(!quiet) {
- cat("IRLS based on variance estimates according to the two component error model\n")
- cat("Initial variance components are:\n")
- print(signif(tc_fitted))
+ if (is.character(tc_fit)) {
+ if (!quiet) {
+ cat(tc_fit, ".\n", "No reweighting will be performed.")
+ }
+ tc_fitted <- c(sigma_low = NA, rsd_high = NA)
+ } else {
+ tc_fitted <- coef(tc_fit)
+ if(!quiet) {
+ cat("IRLS based on variance estimates according to the two component error model\n")
+ cat("Initial variance components are:\n")
+ print(signif(tc_fitted))
+ }
}
}
reweight.diff = 1
@@ -445,7 +434,9 @@ mkinfit <- function(mkinmod, observed,
if (!is.null(err)) observed$err.ini <- observed[[err]]
err = "err.irls"
- while (reweight.diff > reweight.tol & n.iter < reweight.max.iter) {
+ while (reweight.diff > reweight.tol &
+ n.iter < reweight.max.iter &
+ !is.character(tc_fit)) {
n.iter <- n.iter + 1
# Store squared residual predictors used for weighting in sr_old and define new weights
if (reweight.method == "obs") {
@@ -454,7 +445,12 @@ mkinfit <- function(mkinmod, observed,
}
if (reweight.method == "tc") {
sr_old <- tc_fitted
- observed[err] <- predict(tc_fit)
+
+ tmp_predicted <- mkin_wide_to_long(out_predicted, time = "time")
+ tmp_data <- suppressMessages(join(observed, tmp_predicted, by = c("time", "name")))
+
+ #observed[err] <- predict(tc_fit, newdata = data.frame(mod = tmp_data[[4]]))
+ observed[err] <- predict(tc_fit, newdata = data.frame(obs = observed$value))
}
fit <- modFit(cost, fit$par, method = method.modFit,
@@ -464,27 +460,17 @@ mkinfit <- function(mkinmod, observed,
sr_new <- fit$var_ms_unweighted
}
if (reweight.method == "tc") {
- tmp_res <- cost(fit$par)$residuals
- mad_agg <- aggregate(tmp_res$res.unweighted,
- by = list(name = tmp_res$name, x_res = tmp_res$x),
- FUN = function(x) mad(x, center = 0))
- names(mad_agg) <- c("name", "x", "mad")
- tmp_res_mad <- merge(tmp_res, mad_agg)
-
- tc_fit <- try(
- nls(mad ~ sigma_twocomp(mod, sigma_low, rsd_high),
- start = list(sigma_low = tc["sigma_low"], rsd_high = tc["rsd_high"]),
- data = tmp_res_mad,
- lower = 0,
- algorithm = "port"))
-
- if (inherits(tc_fit, "try-error")) {
- stop("Estimation of the two error model components failed during reweighting.\n",
- "Try without reweighting or with reweight.method = 'obs'.")
+ tc_fit <- .fit_error_model_mad_obs(cost(fit$par)$residuals, tc_fitted, n.iter)
+
+ if (is.character(tc_fit)) {
+ if (!quiet) {
+ cat(tc_fit, ".\n")
+ }
+ break
+ } else {
+ tc_fitted <- coef(tc_fit)
+ sr_new <- tc_fitted
}
-
- tc_fitted <- coef(tc_fit)
- sr_new <- tc_fitted
}
reweight.diff = sum((sr_new - sr_old)^2)
@@ -872,4 +858,54 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), .
invisible(x)
}
+
+# Fit the mean absolute deviance against the observed values,
+# using the current error model for weighting
+.fit_error_model_mad_obs <- function(tmp_res, tc, iteration) {
+ mad_agg <- aggregate(tmp_res$res.unweighted,
+ by = list(name = tmp_res$name, time = tmp_res$x),
+ FUN = function(x) mad(x, center = 0))
+ names(mad_agg) <- c("name", "time", "mad")
+ error_data <- suppressMessages(
+ join(data.frame(name = tmp_res$name,
+ time = tmp_res$x,
+ obs = tmp_res$obs),
+ mad_agg))
+ error_data_complete <- na.omit(error_data)
+
+ tc_fit <- tryCatch(
+ nls(mad ~ sigma_twocomp(obs, sigma_low, rsd_high),
+ start = list(sigma_low = tc["sigma_low"], rsd_high = tc["rsd_high"]),
+ weights = 1/sigma_twocomp(error_data_complete$obs,
+ tc["sigma_low"],
+ tc["rsd_high"])^2,
+ data = error_data_complete,
+ lower = 0,
+ algorithm = "port"),
+ error = function(e) paste("Fitting the error model failed in iteration", iteration))
+ return(tc_fit)
+}
+# Alternative way to fit the error model, fitting to modelled instead of
+# observed values
+.fit_error_model_mad_mod <- function(tmp_res, tc) {
+ mad_agg <- aggregate(tmp_res$res.unweighted,
+ by = list(name = tmp_res$name, time = tmp_res$x),
+ FUN = function(x) mad(x, center = 0))
+ names(mad_agg) <- c("name", "time", "mad")
+ mod_agg <- aggregate(tmp_res$mod,
+ by = list(name = tmp_res$name, time = tmp_res$x),
+ FUN = mean)
+ names(mod_agg) <- c("name", "time", "mod")
+ mod_mad <- merge(mod_agg, mad_agg)
+
+ tc_fit <- tryCatch(
+ nls(mad ~ sigma_twocomp(mod, sigma_low, rsd_high),
+ start = list(sigma_low = tc["sigma_low"], rsd_high = tc["rsd_high"]),
+ data = mod_mad,
+ weights = 1/mod_mad$mad,
+ lower = 0,
+ algorithm = "port"),
+ error = "Fitting the error model failed in iteration")
+ return(tc_fit)
+}
# vim: set ts=2 sw=2 expandtab:
diff --git a/check.log b/check.log
index d630f59b..3afdb1b6 100644
--- a/check.log
+++ b/check.log
@@ -5,7 +5,7 @@
* using options ‘--no-tests --as-cran’
* checking for file ‘mkin/DESCRIPTION’ ... OK
* checking extension type ... Package
-* this is package ‘mkin’ version ‘0.9.47.5’
+* this is package ‘mkin’ version ‘0.9.47.6’
* package encoding: UTF-8
* checking CRAN incoming feasibility ... Note_to_CRAN_maintainers
Maintainer: ‘Johannes Ranke <jranke@uni-bremen.de>’
diff --git a/test.log b/test.log
index cdd3bc8c..7754d6ad 100644
--- a/test.log
+++ b/test.log
@@ -2,23 +2,101 @@ Loading mkin
Loading required package: testthat
Testing mkin
✔ | OK F W S | Context
- ⠏ | 0 | Calculation of FOCUS chi2 error levels ⠋ | 1 | Calculation of FOCUS chi2 error levels ⠙ | 2 | Calculation of FOCUS chi2 error levels ✔ | 2 | Calculation of FOCUS chi2 error levels [2.1 s]
- ⠏ | 0 | Results for FOCUS D established in expertise for UBA (Ranke 2014) ⠋ | 1 | Results for FOCUS D established in expertise for UBA (Ranke 2014) ⠙ | 2 | Results for FOCUS D established in expertise for UBA (Ranke 2014) ⠹ | 3 | Results for FOCUS D established in expertise for UBA (Ranke 2014) ⠸ | 4 | Results for FOCUS D established in expertise for UBA (Ranke 2014) ⠼ | 5 | Results for FOCUS D established in expertise for UBA (Ranke 2014) ⠴ | 6 | Results for FOCUS D established in expertise for UBA (Ranke 2014) ⠦ | 7 | Results for FOCUS D established in expertise for UBA (Ranke 2014) ⠧ | 8 | Results for FOCUS D established in expertise for UBA (Ranke 2014) ✔ | 8 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [6.1 s]
- ⠏ | 0 | Iteratively reweighted least squares (IRLS) fitting ⠋ | 1 | Iteratively reweighted least squares (IRLS) fitting ⠙ | 1 1 | Iteratively reweighted least squares (IRLS) fitting ✔ | 1 1 | Iteratively reweighted least squares (IRLS) fitting [7.6 s]
+ ⠏ | 0 | Calculation of FOCUS chi2 error levels ⠋ | 1 | Calculation of FOCUS chi2 error levels ⠙ | 2 | Calculation of FOCUS chi2 error levels ✔ | 2 | Calculation of FOCUS chi2 error levels [2.2 s]
+ ⠏ | 0 | Results for FOCUS D established in expertise for UBA (Ranke 2014) ⠋ | 1 | Results for FOCUS D established in expertise for UBA (Ranke 2014) ⠙ | 2 | Results for FOCUS D established in expertise for UBA (Ranke 2014) ⠹ | 3 | Results for FOCUS D established in expertise for UBA (Ranke 2014) ⠸ | 4 | Results for FOCUS D established in expertise for UBA (Ranke 2014) ⠼ | 5 | Results for FOCUS D established in expertise for UBA (Ranke 2014) ⠴ | 6 | Results for FOCUS D established in expertise for UBA (Ranke 2014) ⠦ | 7 | Results for FOCUS D established in expertise for UBA (Ranke 2014) ⠧ | 8 | Results for FOCUS D established in expertise for UBA (Ranke 2014) ✔ | 8 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [6.5 s]
+ ⠏ | 0 | Iteratively reweighted least squares (IRLS) fitting ⠋ | 0 1 | Iteratively reweighted least squares (IRLS) fitting ⠙ | 1 1 | Iteratively reweighted least squares (IRLS) fitting ⠹ | 1 1 1 | Iteratively reweighted least squares (IRLS) fitting ⠹ | 1 1 1 | Iteratively reweighted least squares (IRLS) fitting ⠹ | 1 1 1 | Iteratively reweighted least squares (IRLS) fitting ⠸ | 1 1 2 | Iteratively reweighted least squares (IRLS) fitting ⠹ | 1 1 1 | Iteratively reweighted least squares (IRLS) fitting ⠹ | 1 1 1 | Iteratively reweighted least squares (IRLS) fitting ⠸ | 1 1 2 | Iteratively reweighted least squares (IRLS) fitting ⠹ | 1 1 1 | Iteratively reweighted least squares (IRLS) fitting ⠼ | 1 1 3 | Iteratively reweighted least squares (IRLS) fitting ⠴ | 1 1 4 | Iteratively reweighted least squares (IRLS) fitting ⠹ | 1 1 1 | Iteratively reweighted least squares (IRLS) fitting ⠹ | 2 1 | Iteratively reweighted least squares (IRLS) fitting ⠸ | 3 1 | Iteratively reweighted least squares (IRLS) fitting ⠼ | 4 1 | Iteratively reweighted least squares (IRLS) fitting ⠴ | 5 1 | Iteratively reweighted least squares (IRLS) fitting ⠦ | 6 1 | Iteratively reweighted least squares (IRLS) fitting ⠧ | 7 1 | Iteratively reweighted least squares (IRLS) fitting ✖ | 7 1 | Iteratively reweighted least squares (IRLS) fitting [172.0 s]
────────────────────────────────────────────────────────────────────────────────
-test_irls.R:45: skip: Reweighting method 'tc' works
-IRLS reweighting with method 'tc' is currently under construction
+test_irls.R:38: error: Reweighting method 'obs' works
+Objekt 'tc_fit' nicht gefunden
+1: mkinfit(m_synth_SFO_lin, SFO_lin_a, reweight.method = "obs", quiet = TRUE) at /home/jranke/git/mkin/tests/testthat/test_irls.R:38
+2: system.time({
+ fit <- modFit(cost, c(state.ini.optim, transparms.optim), method = method.modFit,
+ control = control.modFit, lower = lower, upper = upper, ...)
+ if (!is.null(reweight.method)) {
+ if (!reweight.method %in% c("obs", "tc"))
+ stop("Only reweighting methods 'obs' and 'tc' are implemented")
+ if (reweight.method == "obs") {
+ if (!quiet) {
+ cat("IRLS based on variance estimates for each observed variable\n")
+ cat("Initial variance estimates are:\n")
+ print(signif(fit$var_ms_unweighted, 8))
+ }
+ }
+ if (reweight.method == "tc") {
+ tc_fit <- fit_error_model_mad_obs(cost(fit$par)$residuals, tc, 0)
+ if (is.character(tc_fit)) {
+ if (!quiet) {
+ cat(tc_fit, ".\n", "No reweighting will be performed.")
+ }
+ tc_fitted <- c(sigma_low = NA, rsd_high = NA)
+ }
+ else {
+ tc_fitted <- coef(tc_fit)
+ if (!quiet) {
+ cat("IRLS based on variance estimates according to the two component error model\n")
+ cat("Initial variance components are:\n")
+ print(signif(tc_fitted))
+ }
+ }
+ }
+ reweight.diff = 1
+ n.iter <- 0
+ if (!is.null(err))
+ observed$err.ini <- observed[[err]]
+ err = "err.irls"
+ while (reweight.diff > reweight.tol & n.iter < reweight.max.iter & !is.character(tc_fit)) {
+ n.iter <- n.iter + 1
+ if (reweight.method == "obs") {
+ sr_old <- fit$var_ms_unweighted
+ observed[err] <- sqrt(fit$var_ms_unweighted[as.character(observed$name)])
+ }
+ if (reweight.method == "tc") {
+ sr_old <- tc_fitted
+ tmp_predicted <- mkin_wide_to_long(out_predicted, time = "time")
+ tmp_data <- suppressMessages(join(observed, tmp_predicted, by = c("time",
+ "name")))
+ observed[err] <- predict(tc_fit, newdata = data.frame(obs = observed$value))
+ }
+ fit <- modFit(cost, fit$par, method = method.modFit, control = control.modFit,
+ lower = lower, upper = upper, ...)
+ if (reweight.method == "obs") {
+ sr_new <- fit$var_ms_unweighted
+ }
+ if (reweight.method == "tc") {
+ tc_fit <- fit_error_model_mad_obs(cost(fit$par)$residuals, tc_fitted,
+ n.iter)
+ if (is.character(tc_fit)) {
+ if (!quiet) {
+ cat(tc_fit, ".\n")
+ }
+ break
+ }
+ else {
+ tc_fitted <- coef(tc_fit)
+ sr_new <- tc_fitted
+ }
+ }
+ reweight.diff = sum((sr_new - sr_old)^2)
+ if (!quiet) {
+ cat("Iteration", n.iter, "yields variance estimates:\n")
+ print(signif(sr_new, 8))
+ cat("Sum of squared differences to last variance (component) estimates:",
+ signif(reweight.diff, 2), "\n")
+ }
+ }
+ }
+ }) at /home/jranke/git/mkin/R/mkinfit.R:396
────────────────────────────────────────────────────────────────────────────────
- ⠏ | 0 | Model predictions with mkinpredict ⠋ | 1 | Model predictions with mkinpredict ⠙ | 2 | Model predictions with mkinpredict ⠹ | 3 | Model predictions with mkinpredict ✔ | 3 | Model predictions with mkinpredict [0.3 s]
- ⠏ | 0 | Fitting of parent only models ⠋ | 1 | Fitting of parent only models ⠙ | 2 | Fitting of parent only models ⠹ | 3 | Fitting of parent only models ⠸ | 4 | Fitting of parent only models ⠼ | 5 | Fitting of parent only models ⠴ | 6 | Fitting of parent only models ⠦ | 7 | Fitting of parent only models ⠧ | 8 | Fitting of parent only models ⠇ | 9 | Fitting of parent only models ⠏ | 10 | Fitting of parent only models ⠋ | 11 | Fitting of parent only models ⠙ | 12 | Fitting of parent only models ⠹ | 13 | Fitting of parent only models ⠸ | 14 | Fitting of parent only models ⠼ | 15 | Fitting of parent only models ⠴ | 16 | Fitting of parent only models ⠦ | 17 | Fitting of parent only models ⠧ | 18 | Fitting of parent only models ⠇ | 19 | Fitting of parent only models ⠏ | 20 | Fitting of parent only models ⠋ | 21 | Fitting of parent only models ✔ | 21 | Fitting of parent only models [20.9 s]
- ⠏ | 0 | Complex test case from Schaefer et al. (2007) Piacenza paper ⠋ | 1 | Complex test case from Schaefer et al. (2007) Piacenza paper ⠙ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper ✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [5.1 s]
- ⠏ | 0 | Results for synthetic data established in expertise for UBA (Ranke 2014) ⠋ | 1 | Results for synthetic data established in expertise for UBA (Ranke 2014) ⠙ | 2 | Results for synthetic data established in expertise for UBA (Ranke 2014) ⠹ | 3 | Results for synthetic data established in expertise for UBA (Ranke 2014) ⠸ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) ✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [6.4 s]
+ ⠏ | 0 | Model predictions with mkinpredict ⠋ | 1 | Model predictions with mkinpredict ⠙ | 2 | Model predictions with mkinpredict ⠹ | 3 | Model predictions with mkinpredict ✔ | 3 | Model predictions with mkinpredict [0.4 s]
+ ⠏ | 0 | Fitting of parent only models ⠋ | 1 | Fitting of parent only models ⠙ | 2 | Fitting of parent only models ⠹ | 3 | Fitting of parent only models ⠸ | 4 | Fitting of parent only models ⠼ | 5 | Fitting of parent only models ⠴ | 6 | Fitting of parent only models ⠦ | 7 | Fitting of parent only models ⠧ | 8 | Fitting of parent only models ⠇ | 9 | Fitting of parent only models ⠏ | 10 | Fitting of parent only models ⠋ | 11 | Fitting of parent only models ⠙ | 12 | Fitting of parent only models ⠹ | 13 | Fitting of parent only models ⠸ | 14 | Fitting of parent only models ⠼ | 15 | Fitting of parent only models ⠴ | 16 | Fitting of parent only models ⠦ | 17 | Fitting of parent only models ⠧ | 18 | Fitting of parent only models ⠇ | 19 | Fitting of parent only models ⠏ | 20 | Fitting of parent only models ⠋ | 21 | Fitting of parent only models ✔ | 21 | Fitting of parent only models [22.1 s]
+ ⠏ | 0 | Complex test case from Schaefer et al. (2007) Piacenza paper ⠋ | 1 | Complex test case from Schaefer et al. (2007) Piacenza paper ⠙ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper ✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [5.2 s]
+ ⠏ | 0 | Results for synthetic data established in expertise for UBA (Ranke 2014) ⠋ | 1 | Results for synthetic data established in expertise for UBA (Ranke 2014) ⠙ | 2 | Results for synthetic data established in expertise for UBA (Ranke 2014) ⠹ | 3 | Results for synthetic data established in expertise for UBA (Ranke 2014) ⠸ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) ✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [6.5 s]
⠏ | 0 | Calculation of maximum time weighted average concentrations (TWAs) ⠋ | 1 | Calculation of maximum time weighted average concentrations (TWAs) ⠙ | 2 | Calculation of maximum time weighted average concentrations (TWAs) ⠹ | 3 | Calculation of maximum time weighted average concentrations (TWAs) ⠸ | 4 | Calculation of maximum time weighted average concentrations (TWAs) ⠼ | 5 | Calculation of maximum time weighted average concentrations (TWAs) ⠴ | 6 | Calculation of maximum time weighted average concentrations (TWAs) ⠦ | 7 | Calculation of maximum time weighted average concentrations (TWAs) ⠧ | 8 | Calculation of maximum time weighted average concentrations (TWAs) ✔ | 8 | Calculation of maximum time weighted average concentrations (TWAs) [7.6 s]
══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 56.2 s
+Duration: 222.6 s
-OK: 49
-Failed: 0
+OK: 55
+Failed: 1
Warnings: 0
-Skipped: 1
+Skipped: 0
diff --git a/tests/testthat/test_irls.R b/tests/testthat/test_irls.R
index 65541fb5..5e09912f 100644
--- a/tests/testthat/test_irls.R
+++ b/tests/testthat/test_irls.R
@@ -42,62 +42,88 @@ test_that("Reweighting method 'obs' works", {
test_that("Reweighting method 'tc' works", {
skip_on_cran()
- skip("IRLS reweighting with method 'tc' is currently under construction")
+ # Check if we can approximately obtain the parameters and the error model
+ # components that were used in the data generation
+
+ # Parent only
DFOP <- mkinmod(parent = mkinsub("DFOP"))
sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120)
+ parms_DFOP <- c(k1 = 0.2, k2 = 0.02, g = 0.5)
+ parms_DFOP_optim <- c(parent_0 = 100, parms_DFOP)
d_DFOP <- mkinpredict(DFOP,
- c(k1 = 0.2, k2 = 0.02, g = 0.5),
- c(parent = 100),
+ parms_DFOP, c(parent = 100),
sampling_times)
- d_100 <- add_err(d_DFOP,
+ d_2_100 <- add_err(d_DFOP,
sdfunc = function(x) sigma_twocomp(x, 0.5, 0.07),
- n = 1, reps = 100, digits = 5, LOD = -Inf)
- d_1000 <- add_err(d_DFOP,
+ n = 100, reps = 2, digits = 5, LOD = -Inf)
+ d_100_1 <- add_err(d_DFOP,
sdfunc = function(x) sigma_twocomp(x, 0.5, 0.07),
- n = 1, reps = 1000, digits = 5, LOD = -Inf)
+ n = 1, reps = 100, digits = 5, LOD = -Inf)
+
+ f_2_100 <- mmkin("DFOP", d_2_100, quiet = TRUE,
+ cores = if (Sys.getenv("TRAVIS") != "") 1 else 15)
+ parms_2_100 <- apply(sapply(f_2_100, function(x) x$bparms.optim), 1, mean)
+ parm_errors_2_100 <- (parms_2_100 - parms_DFOP_optim) / parms_DFOP_optim
+ expect_true(all(abs(parm_errors_2_100) < 0.2))
+
+ f_2_100_tc <- mmkin("DFOP", d_2_100, reweight.method = "tc", quiet = TRUE,
+ cores = if (Sys.getenv("TRAVIS") != "") 1 else 15)
+ parms_2_100_tc <- apply(sapply(f_2_100_tc, function(x) x$bparms.optim), 1, mean)
+ parm_errors_2_100_tc <- (parms_2_100_tc - parms_DFOP_optim) / parms_DFOP_optim
+ expect_true(all(abs(parm_errors_2_100_tc) < 0.1))
+
+ tcf_2_100_tc <- apply(sapply(f_2_100_tc, function(x) x$tc_fitted), 1, mean, na.rm = TRUE)
- f_100 <- mkinfit(DFOP, d_100[[1]])
- f_100$bparms.optim
- f_tc_100 <- mkinfit(DFOP, d_100[[1]], reweight.method = "tc")
- f_tc_100$bparms.optim
- f_tc_100$tc_fitted
+ tcf_2_100_error_model_errors <- (tcf_2_100_tc - c(0.5, 0.07)) / c(0.5, 0.07)
+ expect_true(all(abs(tcf_2_100_error_model_errors) < 0.2))
- f_tc_1000 <- mkinfit(DFOP, d_1000[[1]], reweight.method = "tc")
- f_tc_1000$bparms.optim
- f_tc_1000$tc_fitted
+ f_tc_100_1 <- suppressWarnings(mkinfit(DFOP, d_100_1[[1]], reweight.method = "tc", quiet = TRUE))
+ parm_errors_100_1 <- (f_tc_100_1$bparms.optim - parms_DFOP_optim) / parms_DFOP_optim
+ expect_true(all(abs(parm_errors_100_1) < 0.05))
+ tcf_100_1_error_model_errors <- (f_tc_100_1$tc_fitted - c(0.5, 0.07)) /
+ c(0.5, 0.07)
+ # Even with 100 (or even 1000, not shown) replicates at each observation time
+ # we only get a precision of 20% for the error model components
+ expect_true(all(abs(tcf_100_1_error_model_errors) < 0.2))
+
+ # Parent and two metabolites
m_synth_DFOP_lin <- mkinmod(parent = list(type = "DFOP", to = "M1"),
M1 = list(type = "SFO", to = "M2"),
M2 = list(type = "SFO"), use_of_ff = "max",
quiet = TRUE)
sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120)
- d_synth_DFOP_lin <- mkinpredict(m_synth_DFOP_lin,
- c(k1 = 0.2, k2 = 0.02, g = 0.5,
+ parms_DFOP_lin <- c(k1 = 0.2, k2 = 0.02, g = 0.5,
f_parent_to_M1 = 0.5, k_M1 = 0.3,
- f_M1_to_M2 = 0.7, k_M2 = 0.02),
+ f_M1_to_M2 = 0.7, k_M2 = 0.02)
+ d_synth_DFOP_lin <- mkinpredict(m_synth_DFOP_lin,
+ parms_DFOP_lin,
c(parent = 100, M1 = 0, M2 = 0),
sampling_times)
+ parms_DFOP_lin_optim = c(parent_0 = 100, parms_DFOP_lin)
- d_met_100 <- add_err(d_synth_DFOP_lin,
- sdfunc = function(x) sigma_twocomp(x, 0.5, 0.07),
- n = 1, reps = 100, digits = 5, LOD = -Inf)
- d_met_1000 <- add_err(d_synth_DFOP_lin,
+ d_met_2_15 <- add_err(d_synth_DFOP_lin,
sdfunc = function(x) sigma_twocomp(x, 0.5, 0.07),
- n = 1, reps = 1000, digits = 5, LOD = -Inf)
+ n = 15, reps = 1000, digits = 5, LOD = -Inf)
- f_met_100 <- mkinfit(m_synth_DFOP_lin, d_met_100[[1]])
- summary(f_met_100)$bpar
+ time_met_2_15_tc_15 <- system.time(
+ f_met_2_15_tc_e4 <- mmkin(list(m_synth_DFOP_lin), d_met_2_15, quiet = TRUE,
+ reweight.method = "tc", reweight.tol = 1e-4,
+ cores = if (Sys.getenv("TRAVIS") != "") 1 else 15)
+ )
- f_met_100 <- mkinfit(m_synth_DFOP_lin, d_met_100[[1]], reweight.method = "tc")
- summary(f.100)$bpar
+ parms_met_2_15_tc_e4 <- apply(sapply(f_met_2_15_tc_e4, function(x) x$bparms.optim), 1, mean)
+ parm_errors_met_2_15_tc_e4 <- (parms_met_2_15_tc_e4[names(parms_DFOP_lin_optim)] -
+ parms_DFOP_lin_optim) / parms_DFOP_lin_optim
+ expect_true(all(abs(parm_errors_met_2_15_tc_e4) < 0.01))
+ tcf_met_2_15_tc <- apply(sapply(f_met_2_15_tc_e4, function(x) x$tc_fitted), 1, mean, na.rm = TRUE)
- fit_irls_2 <- mkinfit(m_synth_DFOP_par, DFOP_par_c, reweight.method = "tc", quiet = TRUE)
- parms_2 <- signif(fit_irls_2$bparms.optim, 3)
- expect_equivalent(parms_2, c(99.3, 0.041, 0.00962, 0.597, 0.393, 0.298, 0.0203, 0.707))
+ tcf_met_2_15_tc_error_model_errors <- (tcf_met_2_15_tc - c(0.5, 0.07)) /
+ c(0.5, 0.07)
- fit_irls_3 <- mkinfit("DFOP", FOCUS_2006_C, reweight.method = "tc", quiet = TRUE)
- parms_3 <- signif(fit_irls_3$bparms.optim, 3)
- expect_equivalent(parms_3, c(85.0, 0.46, 0.0178, 0.854))
+ # Here we only get a precision < 30% for retrieving the original error model components
+ # from 15 datasets
+ expect_true(all(abs(tcf_met_2_15_tc_error_model_errors) < 0.3))
})
diff --git a/tests/testthat/test_twa.R b/tests/testthat/test_twa.R
index 9151ed42..42b74a7f 100644
--- a/tests/testthat/test_twa.R
+++ b/tests/testthat/test_twa.R
@@ -22,15 +22,9 @@ context("Calculation of maximum time weighted average concentrations (TWAs)")
test_that("Time weighted average concentrations are correct", {
skip_on_cran()
twa_models <- c("SFO", "FOMC", "DFOP", "HS")
- travis_env <- Sys.getenv("TRAVIS")
- if (travis_env == "") {
- travis <- FALSE
- } else {
- travis <- TRUE
- }
fits <- mmkin(twa_models,
list(FOCUS_C = FOCUS_2006_C, FOCUS_D = FOCUS_2006_D),
- quiet = TRUE, cores = if (travis) 1 else 8)
+ quiet = TRUE, cores = if (Sys.getenv("TRAVIS") == "") 15 else 1)
outtimes_10 <- seq(0, 10, length.out = 10000)

Contact - Imprint