From b12e80a875d87f790d67a4e5a50d829060316a18 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Fri, 21 Sep 2018 17:15:06 +0200 Subject: Improve fitting the two-component error model with respect to accuracy and robustness. --- DESCRIPTION | 4 +- NAMESPACE | 3 +- NEWS.md | 6 ++- R/mkinfit.R | 130 +++++++++++++++++++++++++++++---------------- check.log | 2 +- test.log | 104 +++++++++++++++++++++++++++++++----- tests/testthat/test_irls.R | 92 ++++++++++++++++++++------------ tests/testthat/test_twa.R | 8 +-- 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 ’ 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) -- cgit v1.2.1