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. --- R/mkinfit.R | 130 ++++++++++++++++++++++++++++++++++++++---------------------- 1 file changed, 83 insertions(+), 47 deletions(-) (limited to 'R') 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: -- cgit v1.2.1