Improve fitting the two-component error model

with respect to accuracy and robustness.

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)