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

2 files changed, 86 insertions, 14 deletions

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

diff --git a/tests/testthat/test_error_models.R b/tests/testthat/test_error_models.R
index 94703e7f..c656f7d2 100644
--- a/tests/testthat/test_error_models.R
+++ b/tests/testthat/test_error_models.R
@@ -180,7 +180,6 @@ test_that("Reweighting method 'tc' produces reasonable variance estimates", {
 
 test_that("The two-component error model finds the best known AIC values for parent models", {
   skip_on_cran()
-  experimental_data_for_UBA_2019
   library(parallel)
   source("~/git/mkin/R/mkinfit.R")
   source("~/git/mkin/R/mmkin.R")
@@ -195,27 +194,68 @@ test_that("The two-component error model finds the best known AIC values for par
                  error_model = "tc", error_model_algorithm = "fourstep")
   f_9 <- mkinfit("SFO", experimental_data_for_UBA_2019[[9]]$data,
                  error_model = "tc", error_model_algorithm = "IRLS")
+  f_9 <- mkinfit("SFO", experimental_data_for_UBA_2019[[9]]$data,
+                 error_model = "tc", error_model_algorithm = "d_3")
   AIC(f_9)
-  f_tc_exp <- mmkin(c("SFO"),
+  f_10 <- mkinfit("DFOP", experimental_data_for_UBA_2019[[10]]$data,
+                 error_model = "tc", error_model_algorithm = "IRLS")
+  f_tc_exp_direct <- mmkin(c("SFO", "DFOP", "HS"),
     lapply(experimental_data_for_UBA_2019, function(x) x$data),
     error_model = "tc",
     error_model_algorithm = "direct",
     quiet = TRUE)
-  f_tc_exp <- mmkin(c("SFO"),
+  f_tc_exp_twostep <- mmkin(c("SFO", "DFOP", "HS"),
     lapply(experimental_data_for_UBA_2019, function(x) x$data),
     error_model = "tc",
     error_model_algorithm = "twostep",
     quiet = TRUE)
-  f_tc_exp <- mmkin(c("SFO"),
+  f_tc_exp_threestep <- mmkin(c("SFO", "DFOP", "HS"),
     lapply(experimental_data_for_UBA_2019, function(x) x$data),
     error_model = "tc",
     error_model_algorithm = "threestep",
     quiet = TRUE)
+  f_tc_exp_fourstep <- mmkin(c("SFO", "DFOP", "HS"),
+    lapply(experimental_data_for_UBA_2019, function(x) x$data),
+    error_model = "tc",
+    error_model_algorithm = "fourstep",
+    quiet = TRUE)
+  f_tc_exp_IRLS <- mmkin(c("SFO", "DFOP", "HS"),
+    lapply(experimental_data_for_UBA_2019, function(x) x$data),
+    error_model = "tc",
+    error_model_algorithm = "IRLS",
+    quiet = TRUE)
+
+  AIC_exp_direct <- lapply(f_tc_exp_direct, AIC)
+  AIC_exp_direct <- lapply(AIC_exp_direct, round, 1)
+  dim(AIC_exp_direct) <- dim(f_tc_exp_direct)
+  dimnames(AIC_exp_direct) <- dimnames(f_tc_exp_direct)
+
+  AIC_exp_twostep <- lapply(f_tc_exp_twostep, AIC)
+  AIC_exp_twostep <- lapply(AIC_exp_twostep, round, 1)
+  dim(AIC_exp_twostep) <- dim(f_tc_exp_twostep)
+  dimnames(AIC_exp_twostep) <- dimnames(f_tc_exp_twostep)
+
+  AIC_exp_threestep <- lapply(f_tc_exp_threestep, AIC)
+  AIC_exp_threestep <- lapply(AIC_exp_threestep, round, 1)
+  dim(AIC_exp_threestep) <- dim(f_tc_exp_threestep)
+  dimnames(AIC_exp_threestep) <- dimnames(f_tc_exp_threestep)
+
+  AIC_exp_fourstep <- lapply(f_tc_exp_fourstep, AIC)
+  AIC_exp_fourstep <- lapply(AIC_exp_fourstep, round, 1)
+  dim(AIC_exp_fourstep) <- dim(f_tc_exp_fourstep)
+  dimnames(AIC_exp_fourstep) <- dimnames(f_tc_exp_fourstep)
+
+  AIC_exp_IRLS <- lapply(f_tc_exp_IRLS, AIC)
+  AIC_exp_IRLS <- lapply(AIC_exp_IRLS, round, 1)
+  dim(AIC_exp_IRLS) <- dim(f_tc_exp_IRLS)
+  dimnames(AIC_exp_IRLS) <- dimnames(f_tc_exp_IRLS)
 
   AIC_exp <- lapply(f_tc_exp, AIC)
   dim(AIC_exp) <- dim(f_tc_exp)
   dimnames(AIC_exp) <- dimnames(f_tc_exp)
-  expect_equivalent(round(AIC_exp["SFO", c(9, 11, 12)], 1), c(134.9, 125.5, 82.0))
+  unlist(AIC_exp["SFO", c(9, 11, 12)])
+  expect_equivalent(round(unlist(AIC_exp["SFO", c(9, 11, 12)]), 1),
+                    c(134.9, 125.5, 82.0))
 })