diff options
Diffstat (limited to 'R/mkinfit.R')
| -rw-r--r-- | R/mkinfit.R | 130 |
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:
|