From b6ea4f22fc1b6d1caea29f6b1e44774d14d6697c Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Mon, 3 Jun 2019 07:56:44 +0200
Subject: Status von Samstag morgen - untested

---
 R/mkinfit.R | 100 ++++++++++++++++++++++++++++++++++++++----------------------
 1 file changed, 63 insertions(+), 37 deletions(-)

(limited to 'R')

diff --git a/R/mkinfit.R b/R/mkinfit.R
index bc8b9d11..224dc7bd 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -34,6 +34,7 @@ mkinfit <- function(mkinmod, observed,
   quiet = FALSE,
   atol = 1e-8, rtol = 1e-10, n.outtimes = 100,
   error_model = c("const", "obs", "tc"),
+  error_model_algorithm = c("direct", "twostep", "threestep", "fourstep", "IRLS"),
   trace_parms = FALSE,
   ...)
 {
@@ -91,7 +92,7 @@ mkinfit <- function(mkinmod, observed,
   if (length(wrongpar.names) > 0) {
     warning("Initial parameter(s) ", paste(wrongpar.names, collapse = ", "),
          " not used in the model")
-    parms.ini <- parms.ini[setdiff(names(parms.ini), wrongpar.names)] 
+    parms.ini <- parms.ini[setdiff(names(parms.ini), wrongpar.names)]
   }
 
   # Warn that the sum of formation fractions may exceed one if they are not
@@ -244,6 +245,7 @@ mkinfit <- function(mkinmod, observed,
 
   # Get the error model
   err_mod <- match.arg(error_model)
+  error_model_algorithm = match.arg(error_model_algorithm)
   errparm_names <- switch(err_mod,
     "const" = "sigma",
     "obs" = paste0("sigma_", obs_vars),
@@ -276,34 +278,48 @@ mkinfit <- function(mkinmod, observed,
                                               length.out = n.outtimes))))
 
   # Define log-likelihood function for optimisation, including (back)transformations
-  nlogLik <- function(P, trans = TRUE, OLS = FALSE, local = FALSE, update_data = TRUE, ...)
+  nlogLik <- function(P, trans = TRUE, OLS = FALSE, fixed_degparms = FALSE, fixed_errparms = FALSE, update_data = TRUE, ...)
   {
     assign("calls", calls + 1, inherits = TRUE) # Increase the model solution counter
-    P.orig <- P
 
     # Trace parameter values if requested and if we are actually optimising
     if(trace_parms & update_data) cat(P, "\n")
 
-    # If we do a local optimisation of the error model, the initials
-    # for the state variabels and the parameters are given as 'local'
-    if (local[1] != FALSE) {
-      P <- local
+    if (fixed_degparms[1] != FALSE) {
+      degparms <- fixed_degparms
+      errparms <- P # This version of errparms is local to the function
+      degparms_fixed = TRUE
+    } else {
+      degparms_fixed = FALSE
+    }
+
+    if (fixed_errparms[1] != FALSE) {
+      degparms <- P
+      errparms <- fixed_errparms # Local to the function
+      errparms_fixed = TRUE
+    } else {
+      errparms_fixed = FALSE
+    }
+
+    if (OLS) {
+      degparms <- P
+    }
+
+    if (!OLS & !degparms_fixed & !errparms_fixed) {
+      degparms <- P[1:(length(P) - length(errparms))]
+      errparms <- P[(length(degparms) + 1):length(P)]
     }
 
     # Initial states for t0
     if(length(state.ini.optim) > 0) {
-      odeini <- c(P[1:length(state.ini.optim)], state.ini.fixed)
+      odeini <- c(degparms[1:length(state.ini.optim)], state.ini.fixed)
       names(odeini) <- c(state.ini.optim.boxnames, state.ini.fixed.boxnames)
     } else {
       odeini <- state.ini.fixed
       names(odeini) <- state.ini.fixed.boxnames
     }
 
-    if (OLS | identical(P, local)) {
-      odeparms.optim <- P[(length(state.ini.optim) + 1):length(P)]
-    } else {
-      odeparms.optim <- P[(length(state.ini.optim) + 1):(length(P) - length(errparms))]
-    }
+    odeparms.optim <- degparms[(length(state.ini.optim) + 1):length(degparms)]
 
     if (trans == TRUE) {
       odeparms <- c(odeparms.optim, transparms.fixed)
@@ -322,23 +338,20 @@ mkinfit <- function(mkinmod, observed,
                        method.ode = method.ode,
                        atol = atol, rtol = rtol, ...)
 
-    # Get back the original parameter vector to get the error model parameters
-    P <- P.orig
-
     out_long <- mkin_wide_to_long(out, time = "time")
 
     if (err_mod == "const") {
-      observed$std <- P["sigma"]
+      observed$std <- errparms["sigma"]
     }
     if (err_mod == "obs") {
       std_names <- paste0("sigma_", observed$name)
-      observed$std <- P[std_names]
+      observed$std <- errparms[std_names]
     }
     if (err_mod == "tc") {
       tmp <- merge(observed, out_long, by = c("time", "name"))
       tmp$name <- ordered(tmp$name, levels = obs_vars)
       tmp <- tmp[order(tmp$name, tmp$time), ]
-      observed$std <- sqrt(P["sigma_low"]^2 + tmp$value.y^2 * P["rsd_high"]^2)
+      observed$std <- sqrt(errparms["sigma_low"]^2 + tmp$value.y^2 * errparms["rsd_high"]^2)
     }
 
     data_log_lik <- merge(observed[c("name", "time", "value", "std")], out_long,
@@ -426,31 +439,44 @@ mkinfit <- function(mkinmod, observed,
     # In the inital run, we assume a constant standard deviation and do
     # not optimise it, as this is unnecessary and leads to instability of the
     # fit (most obvious when fitting the IORE model)
-    if (!quiet) message("Ordinary least squares optimisation")
-    parms.start <- c(state.ini.optim, transparms.optim)
-    fit.ols <- nlminb(parms.start, nlogLik, control = control,
-      lower = lower[names(parms.start)],
-      upper = upper[names(parms.start)], OLS = TRUE, ...)
-
-    if (err_mod == "const") {
+    degparms <- c(state.ini.optim, transparms.optim)
+
+    if (err_mod == "const" | error_model_algorithm != "direct") {
+      if (!quiet) message("Ordinary least squares optimisation")
+      fit <- nlminb(degparms, nlogLik, control = control,
+        lower = lower[names(degparms)],
+        upper = upper[names(degparms)], OLS = TRUE, ...)
+      degparms <- fit$par
       # Get the maximum likelihood estimate for sigma at the optimum parameter values
       data_errmod$residual <- data_errmod$value.observed - data_errmod$value.predicted
-      sigma_mle = sqrt(sum(data_errmod$residual^2)/nrow(data_errmod))
+      sigma_mle <- sqrt(sum(data_errmod$residual^2)/nrow(data_errmod))
 
-      errparms = c(sigma = sigma_mle)
-      fit <- fit.ols
-      fit$logLik <- - nlogLik(c(fit$par, sigma = sigma_mle), OLS = FALSE)
-    } else {
-      # After the OLS stop we fit the error model with fixed degradation model
-      # parameters
+      if (err_mod == "const") {
+        errparms <- c(sigma = sigma_mle)
+      }
+
+      nlogLik.current <- nlogLik(c(degparms, errparms), OLS = FALSE)
+      fit$logLik <- - nlogLik.current
+    }
+    if (error_model_algorithm %in% c("threestep", "fourstep")) {
       if (!quiet) message("Optimising the error model")
-      fit.err <- nlminb(errparms, nlogLik, control = control,
+      fit <- nlminb(errparms, nlogLik, control = control,
         lower = lower[names(errparms)],
         upper = upper[names(errparms)],
-        local = fit.ols$par, ...)
-      errparms.tmp <- fit.err$par
+        degparms_fixed = degparms, ...)
+      errparms <- fit$par
+    }
+    if (error_model_algorithm == "fourstep") {
+      if (!quiet) message("Optimising the degradation model")
+      fit <- nlminb(degparms, nlogLik, control = control,
+        lower = lower[names(degparms)],
+        upper = upper[names(degparms)],
+        errparms_fixed = errparms, ...)
+      degparms <- fit$par
+    }
+    if (error_model_algorithm %in% c("direct", "twostep", "threestep", "fourstep")) {
       if (!quiet) message("Optimising the complete model")
-      parms.start <- c(fit.ols$par, errparms.tmp)
+      parms.start <- c(degparms, errparms)
       fit <- nlminb(parms.start, nlogLik,
         lower = lower[names(parms.start)],
         upper = upper[names(parms.start)],
-- 
cgit v1.2.1


From 9a96391589fef9f80f9c6c4881cc48a509cb75f2 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Tue, 4 Jun 2019 11:15:52 +0200
Subject: Algorithms direct, two-, three-, fourstep, IRLS

All of them are working now and allow for comparison
Based on SFO, DFOP and HS fits to twelve test datasets, only
the combination of direct and threestep is needed to find
the lowest AIC
---
 R/mkinfit.R | 50 +++++++++++++++++++++++++++++++++++++++++---------
 1 file changed, 41 insertions(+), 9 deletions(-)

(limited to 'R')

diff --git a/R/mkinfit.R b/R/mkinfit.R
index 224dc7bd..e0a0e525 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -35,6 +35,7 @@ mkinfit <- function(mkinmod, observed,
   atol = 1e-8, rtol = 1e-10, n.outtimes = 100,
   error_model = c("const", "obs", "tc"),
   error_model_algorithm = c("direct", "twostep", "threestep", "fourstep", "IRLS"),
+  reweight.tol = 1e-8, reweight.max.iter = 10,
   trace_parms = FALSE,
   ...)
 {
@@ -285,7 +286,7 @@ mkinfit <- function(mkinmod, observed,
     # Trace parameter values if requested and if we are actually optimising
     if(trace_parms & update_data) cat(P, "\n")
 
-    if (fixed_degparms[1] != FALSE) {
+    if (is.numeric(fixed_degparms)) {
       degparms <- fixed_degparms
       errparms <- P # This version of errparms is local to the function
       degparms_fixed = TRUE
@@ -293,7 +294,7 @@ mkinfit <- function(mkinmod, observed,
       degparms_fixed = FALSE
     }
 
-    if (fixed_errparms[1] != FALSE) {
+    if (is.numeric(fixed_errparms)) {
       degparms <- P
       errparms <- fixed_errparms # Local to the function
       errparms_fixed = TRUE
@@ -364,7 +365,8 @@ mkinfit <- function(mkinmod, observed,
         sum(dnorm(x = value.observed, mean = value.predicted, sd = std, log = TRUE)))
     }
 
-    # We update the current likelihood and data during the optimisation, not during hessian calculations
+    # We update the current likelihood and data during the optimisation, not
+    # during hessian calculations
     if (update_data) {
 
       assign("out_predicted", out_long, inherits = TRUE)
@@ -434,11 +436,10 @@ mkinfit <- function(mkinmod, observed,
   # Show parameter names if tracing is requested
   if(trace_parms) cat(names_optim, "\n")
 
+  # browser()
+
   # Do the fit and take the time until the hessians are calculated
   fit_time <- system.time({
-    # In the inital run, we assume a constant standard deviation and do
-    # not optimise it, as this is unnecessary and leads to instability of the
-    # fit (most obvious when fitting the IORE model)
     degparms <- c(state.ini.optim, transparms.optim)
 
     if (err_mod == "const" | error_model_algorithm != "direct") {
@@ -463,7 +464,7 @@ mkinfit <- function(mkinmod, observed,
       fit <- nlminb(errparms, nlogLik, control = control,
         lower = lower[names(errparms)],
         upper = upper[names(errparms)],
-        degparms_fixed = degparms, ...)
+        fixed_degparms = degparms, ...)
       errparms <- fit$par
     }
     if (error_model_algorithm == "fourstep") {
@@ -471,10 +472,11 @@ mkinfit <- function(mkinmod, observed,
       fit <- nlminb(degparms, nlogLik, control = control,
         lower = lower[names(degparms)],
         upper = upper[names(degparms)],
-        errparms_fixed = errparms, ...)
+        fixed_errparms = errparms, ...)
       degparms <- fit$par
     }
-    if (error_model_algorithm %in% c("direct", "twostep", "threestep", "fourstep")) {
+    if (error_model_algorithm %in% c("direct", "twostep", "threestep", "fourstep") &
+        err_mod != "const") {
       if (!quiet) message("Optimising the complete model")
       parms.start <- c(degparms, errparms)
       fit <- nlminb(parms.start, nlogLik,
@@ -483,6 +485,36 @@ mkinfit <- function(mkinmod, observed,
         control = control, ...)
       fit$logLik <- - nlogLik.current
     }
+    if (err_mod != "const" & error_model_algorithm == "IRLS") {
+      reweight.diff <- 1
+      n.iter <- 0
+      errparms_last <- errparms
+
+      while (reweight.diff > reweight.tol &
+             n.iter < reweight.max.iter) {
+
+        if (!quiet) message("Optimising the error model")
+        fit <- nlminb(errparms, nlogLik, control = control,
+          lower = lower[names(errparms)],
+          upper = upper[names(errparms)],
+          fixed_degparms = degparms, ...)
+        errparms <- fit$par
+
+        if (!quiet) message("Optimising the degradation model")
+        fit <- nlminb(degparms, nlogLik, control = control,
+          lower = lower[names(degparms)],
+          upper = upper[names(degparms)],
+          fixed_errparms = errparms, ...)
+        degparms <- fit$par
+
+        reweight.diff <- dist(rbind(errparms, errparms_last))
+        errparms_last <- errparms
+
+        fit$par <- c(fit$par, errparms)
+        nlogLik.current <- nlogLik(c(degparms, errparms), OLS = FALSE)
+        fit$logLik <- - nlogLik.current
+      }
+    }
 
     # We include the error model in the parameter uncertainty analysis, also
     # for constant variance, to get a confidence interval for it
-- 
cgit v1.2.1


From 95178837d3f91e84837628446b5fd468179af2b9 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Tue, 4 Jun 2019 15:09:28 +0200
Subject: Additional algorithm "d_c", more tests, docs

The new algorithm tries direct optimization of the likelihood, as well
as a three step procedure. In this way, we consistently get the
model with the highest likelihood for SFO, DFOP and HS for all 12
new test datasets.
---
 R/mkinfit.R | 117 ++++++++++++++++++++++++++++++++++++------------------------
 1 file changed, 71 insertions(+), 46 deletions(-)

(limited to 'R')

diff --git a/R/mkinfit.R b/R/mkinfit.R
index e0a0e525..33e13d8e 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -34,7 +34,7 @@ mkinfit <- function(mkinmod, observed,
   quiet = FALSE,
   atol = 1e-8, rtol = 1e-10, n.outtimes = 100,
   error_model = c("const", "obs", "tc"),
-  error_model_algorithm = c("direct", "twostep", "threestep", "fourstep", "IRLS"),
+  error_model_algorithm = c("d_3", "direct", "twostep", "threestep", "fourstep", "IRLS"),
   reweight.tol = 1e-8, reweight.max.iter = 10,
   trace_parms = FALSE,
   ...)
@@ -442,77 +442,102 @@ mkinfit <- function(mkinmod, observed,
   fit_time <- system.time({
     degparms <- c(state.ini.optim, transparms.optim)
 
-    if (err_mod == "const" | error_model_algorithm != "direct") {
+    if (err_mod == "const") {
       if (!quiet) message("Ordinary least squares optimisation")
       fit <- nlminb(degparms, nlogLik, control = control,
         lower = lower[names(degparms)],
         upper = upper[names(degparms)], OLS = TRUE, ...)
       degparms <- fit$par
+
       # Get the maximum likelihood estimate for sigma at the optimum parameter values
       data_errmod$residual <- data_errmod$value.observed - data_errmod$value.predicted
       sigma_mle <- sqrt(sum(data_errmod$residual^2)/nrow(data_errmod))
 
-      if (err_mod == "const") {
-        errparms <- c(sigma = sigma_mle)
-      }
-
+      errparms <- c(sigma = sigma_mle)
       nlogLik.current <- nlogLik(c(degparms, errparms), OLS = FALSE)
       fit$logLik <- - nlogLik.current
-    }
-    if (error_model_algorithm %in% c("threestep", "fourstep")) {
-      if (!quiet) message("Optimising the error model")
-      fit <- nlminb(errparms, nlogLik, control = control,
-        lower = lower[names(errparms)],
-        upper = upper[names(errparms)],
-        fixed_degparms = degparms, ...)
-      errparms <- fit$par
-    }
-    if (error_model_algorithm == "fourstep") {
-      if (!quiet) message("Optimising the degradation model")
-      fit <- nlminb(degparms, nlogLik, control = control,
-        lower = lower[names(degparms)],
-        upper = upper[names(degparms)],
-        fixed_errparms = errparms, ...)
-      degparms <- fit$par
-    }
-    if (error_model_algorithm %in% c("direct", "twostep", "threestep", "fourstep") &
-        err_mod != "const") {
-      if (!quiet) message("Optimising the complete model")
-      parms.start <- c(degparms, errparms)
-      fit <- nlminb(parms.start, nlogLik,
-        lower = lower[names(parms.start)],
-        upper = upper[names(parms.start)],
-        control = control, ...)
-      fit$logLik <- - nlogLik.current
-    }
-    if (err_mod != "const" & error_model_algorithm == "IRLS") {
-      reweight.diff <- 1
-      n.iter <- 0
-      errparms_last <- errparms
-
-      while (reweight.diff > reweight.tol &
-             n.iter < reweight.max.iter) {
+    } else {
+      if (error_model_algorithm == "d_3") {
+        if (!quiet) message("Directly optimising the complete model")
+        parms.start <- c(degparms, errparms)
+        fit_direct <- nlminb(parms.start, nlogLik,
+          lower = lower[names(parms.start)],
+          upper = upper[names(parms.start)],
+          control = control, ...)
+        fit_direct$logLik <- - nlogLik.current
+        nlogLik.current <- Inf # reset to avoid conflict with the OLS step
+      }
+      if (error_model_algorithm != "direct") {
+        if (!quiet) message("Ordinary least squares optimisation")
+        fit <- nlminb(degparms, nlogLik, control = control,
+          lower = lower[names(degparms)],
+          upper = upper[names(degparms)], OLS = TRUE, ...)
+        degparms <- fit$par
+        # Get the maximum likelihood estimate for sigma at the optimum parameter values
+        data_errmod$residual <- data_errmod$value.observed - data_errmod$value.predicted
+        sigma_mle <- sqrt(sum(data_errmod$residual^2)/nrow(data_errmod))
 
+        nlogLik.current <- nlogLik(c(degparms, errparms), OLS = FALSE)
+        fit$logLik <- - nlogLik.current
+      }
+      if (error_model_algorithm %in% c("threestep", "fourstep", "d_3")) {
         if (!quiet) message("Optimising the error model")
         fit <- nlminb(errparms, nlogLik, control = control,
           lower = lower[names(errparms)],
           upper = upper[names(errparms)],
           fixed_degparms = degparms, ...)
         errparms <- fit$par
-
+      }
+      if (error_model_algorithm == "fourstep") {
         if (!quiet) message("Optimising the degradation model")
         fit <- nlminb(degparms, nlogLik, control = control,
           lower = lower[names(degparms)],
           upper = upper[names(degparms)],
           fixed_errparms = errparms, ...)
         degparms <- fit$par
-
-        reweight.diff <- dist(rbind(errparms, errparms_last))
+      }
+      if (error_model_algorithm %in% c("direct", "twostep", "threestep",
+                                       "fourstep", "d_3")) {
+        if (!quiet) message("Optimising the complete model")
+        parms.start <- c(degparms, errparms)
+        fit <- nlminb(parms.start, nlogLik,
+          lower = lower[names(parms.start)],
+          upper = upper[names(parms.start)],
+          control = control, ...)
+        fit$logLik <- - nlogLik.current
+        if (error_model_algorithm == "d_3" && fit_direct$logLik > fit$logLik) {
+          fit <- fit_direct
+        }
+      }
+      if (err_mod != "const" & error_model_algorithm == "IRLS") {
+        reweight.diff <- 1
+        n.iter <- 0
         errparms_last <- errparms
 
-        fit$par <- c(fit$par, errparms)
-        nlogLik.current <- nlogLik(c(degparms, errparms), OLS = FALSE)
-        fit$logLik <- - nlogLik.current
+        while (reweight.diff > reweight.tol &
+               n.iter < reweight.max.iter) {
+
+          if (!quiet) message("Optimising the error model")
+          fit <- nlminb(errparms, nlogLik, control = control,
+            lower = lower[names(errparms)],
+            upper = upper[names(errparms)],
+            fixed_degparms = degparms, ...)
+          errparms <- fit$par
+
+          if (!quiet) message("Optimising the degradation model")
+          fit <- nlminb(degparms, nlogLik, control = control,
+            lower = lower[names(degparms)],
+            upper = upper[names(degparms)],
+            fixed_errparms = errparms, ...)
+          degparms <- fit$par
+
+          reweight.diff <- dist(rbind(errparms, errparms_last))
+          errparms_last <- errparms
+
+          fit$par <- c(fit$par, errparms)
+          nlogLik.current <- nlogLik(c(degparms, errparms), OLS = FALSE)
+          fit$logLik <- - nlogLik.current
+        }
       }
     }
 
-- 
cgit v1.2.1


From 2bc4adb7080e5893ab423768fe2e24777b292f19 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Tue, 4 Jun 2019 19:46:55 +0200
Subject: For the d_3 algorithm, report which was better, if any

---
 R/mkinfit.R | 15 +++++++++++++--
 1 file changed, 13 insertions(+), 2 deletions(-)

(limited to 'R')

diff --git a/R/mkinfit.R b/R/mkinfit.R
index 33e13d8e..13c7aa90 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -505,8 +505,19 @@ mkinfit <- function(mkinmod, observed,
           upper = upper[names(parms.start)],
           control = control, ...)
         fit$logLik <- - nlogLik.current
-        if (error_model_algorithm == "d_3" && fit_direct$logLik > fit$logLik) {
-          fit <- fit_direct
+        if (error_model_algorithm == "d_3") {
+          if (abs((fit_direct$logLik - fit$logLik))/mean(c(fit_direct$logLik, fit$logLik)) < 0.001) {
+            if (!quiet) {
+              message("Direct fitting and three-step fitting yield approximately the same likelihood")
+            }
+          } else {
+            if (fit_direct$logLik < fit$logLik) {
+              if (!quiet) message("Three-step fitting yielded a higher likelihood than direct fitting")
+            } else {
+              if (!quiet) message("Direct fitting yielded a higher likelihood than three-step fitting")
+              fit <- fit_direct
+            }
+          }
         }
       }
       if (err_mod != "const" & error_model_algorithm == "IRLS") {
-- 
cgit v1.2.1


From 307a317666b8a1cdfe2293371ad8671403680a36 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Tue, 4 Jun 2019 21:10:58 +0200
Subject: Fix a bug introduced in the last commit

---
 R/mkinfit.R | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

(limited to 'R')

diff --git a/R/mkinfit.R b/R/mkinfit.R
index 13c7aa90..60697cb1 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -505,13 +505,15 @@ mkinfit <- function(mkinmod, observed,
           upper = upper[names(parms.start)],
           control = control, ...)
         fit$logLik <- - nlogLik.current
+
         if (error_model_algorithm == "d_3") {
-          if (abs((fit_direct$logLik - fit$logLik))/mean(c(fit_direct$logLik, fit$logLik)) < 0.001) {
+          rel_diff <- abs((fit_direct$logLik - fit$logLik))/-mean(c(fit_direct$logLik, fit$logLik))
+          if (rel_diff < 0.0001) {
             if (!quiet) {
               message("Direct fitting and three-step fitting yield approximately the same likelihood")
             }
           } else {
-            if (fit_direct$logLik < fit$logLik) {
+            if (fit$logLik > fit_direct$logLik) {
               if (!quiet) message("Three-step fitting yielded a higher likelihood than direct fitting")
             } else {
               if (!quiet) message("Direct fitting yielded a higher likelihood than three-step fitting")
-- 
cgit v1.2.1


From b6027bbd157734e1c7f8c3ba6373451f5c85fc38 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Wed, 5 Jun 2019 15:16:59 +0200
Subject: Add error model algorithm to output

---
 R/mkinfit.R | 27 +++++++++++++++++++--------
 1 file changed, 19 insertions(+), 8 deletions(-)

(limited to 'R')

diff --git a/R/mkinfit.R b/R/mkinfit.R
index 60697cb1..2af4e493 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -506,18 +506,23 @@ mkinfit <- function(mkinmod, observed,
           control = control, ...)
         fit$logLik <- - nlogLik.current
 
+        d_3_messages = c(
+           same = "Direct fitting and three-step fitting yield approximately the same likelihood",
+           threestep = "Three-step fitting yielded a higher likelihood than direct fitting",
+           direct = "Direct fitting yielded a higher likelihood than three-step fitting")
         if (error_model_algorithm == "d_3") {
           rel_diff <- abs((fit_direct$logLik - fit$logLik))/-mean(c(fit_direct$logLik, fit$logLik))
           if (rel_diff < 0.0001) {
-            if (!quiet) {
-              message("Direct fitting and three-step fitting yield approximately the same likelihood")
-            }
+            if (!quiet) message(d_3_messages["same"])
+            fit$d_3_message <- d_3_messages["same"]
           } else {
             if (fit$logLik > fit_direct$logLik) {
-              if (!quiet) message("Three-step fitting yielded a higher likelihood than direct fitting")
+              if (!quiet) message(d_3_messages["threestep"])
+              fit$d_3_message <- d_3_messages["threestep"]
             } else {
-              if (!quiet) message("Direct fitting yielded a higher likelihood than three-step fitting")
+              if (!quiet) message(d_3_messages["direct"])
               fit <- fit_direct
+              fit$d_3_message <- d_3_messages["direct"]
             }
           }
         }
@@ -553,6 +558,7 @@ mkinfit <- function(mkinmod, observed,
         }
       }
     }
+    fit$error_model_algorithm <- error_model_algorithm
 
     # We include the error model in the parameter uncertainty analysis, also
     # for constant variance, to get a confidence interval for it
@@ -725,6 +731,7 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05,
     solution_type = object$solution_type,
     warning = object$warning,
     use_of_ff = object$mkinmod$use_of_ff,
+    error_model_algorithm = object$error_model_algorithm,
     df = c(p, rdf),
     cov.unscaled = covar,
     err_mod = object$err_mod,
@@ -763,8 +770,9 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05,
   ep <- endpoints(object)
   if (length(ep$ff) != 0)
     ans$ff <- ep$ff
-  if(distimes) ans$distimes <- ep$distimes
-  if(length(ep$SFORB) != 0) ans$SFORB <- ep$SFORB
+  if (distimes) ans$distimes <- ep$distimes
+  if (length(ep$SFORB) != 0) ans$SFORB <- ep$SFORB
+  if (!is.null(object$d_3_message)) ans$d_3_message <- object$d_3_message
   class(ans) <- c("summary.mkinfit", "summary.modFit")
   return(ans)
 }
@@ -794,12 +802,15 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), .
 
   cat("\nFitted using", x$calls, "model solutions performed in", x$time[["elapsed"]],  "s\n")
 
-  cat("\nError model:\n")
+  cat("\nError model: ")
   cat(switch(x$err_mod,
              const = "Constant variance",
              obs = "Variance unique to each observed variable",
              tc = "Two-component variance function"), "\n")
 
+  cat("\nError model algorithm:", x$error_model_algorithm, "\n")
+  if (!is.null(x$d_3_message)) cat(x$d_3_message, "\n")
+
   cat("\nStarting values for parameters to be optimised:\n")
   print(x$start)
 
-- 
cgit v1.2.1