Improve fitting the two-component error model

with respect to accuracy and robustness.

1 files changed, 91 insertions, 13 deletions

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