Improve fitting the two-component error model

with respect to accuracy and robustness.

1 files changed, 83 insertions, 47 deletions

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: