From 4ba08177631866e8a910f561800f8a45b23e5cb9 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Thu, 20 Sep 2018 18:19:40 +0200
Subject: Use the median absolute deviance for fitting

as the absolute value is a biased estimator for the standard deviation
---
 R/mkinfit.R | 21 ++++++++++++++++-----
 1 file changed, 16 insertions(+), 5 deletions(-)

diff --git a/R/mkinfit.R b/R/mkinfit.R
index 2f05e81a..d56c663a 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -415,10 +415,16 @@ mkinfit <- function(mkinmod, observed,
         # 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(abs(res.unweighted) ~ sigma_twocomp(obs, sigma_low, rsd_high),
+          nls(mad ~ sigma_twocomp(mod, sigma_low, rsd_high),
             start = list(sigma_low = tc["sigma_low"], rsd_high = tc["rsd_high"]),
-            data = tmp_res,
+            data = tmp_res_mad,
             lower = 0,
             algorithm = "port"))
 
@@ -429,7 +435,7 @@ mkinfit <- function(mkinmod, observed,
 
         tc_fitted <- coef(tc_fit)
         if(!quiet) {
-          cat("IRLS based on variance estimates according to the Rocke and Lorenzato model\n")
+          cat("IRLS based on variance estimates according to the two component error model\n")
           cat("Initial variance components are:\n")
           print(signif(tc_fitted))
         }
@@ -459,11 +465,16 @@ mkinfit <- function(mkinmod, observed,
         }
         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(abs(res.unweighted) ~ sigma_twocomp(obs, sigma_low, rsd_high),
+            nls(mad ~ sigma_twocomp(mod, sigma_low, rsd_high),
               start = list(sigma_low = tc["sigma_low"], rsd_high = tc["rsd_high"]),
-              data = tmp_res,
+              data = tmp_res_mad,
               lower = 0,
               algorithm = "port"))
 
-- 
cgit v1.2.1