From 1ef7008be2a72a0847064ad9c2ddcfa16b055482 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Fri, 3 May 2019 19:14:15 +0200
Subject: Improve error model fitting

Now we have a three stage fitting process for
nonconstant error models:
- Unweighted least squares
- Only optimize the error model
- Optimize both

Static documentation rebuilt by pkgdown
---
 R/mkinds.R  |  1 -
 R/mkinfit.R | 50 ++++++++++++++++++++++++++++++++++++--------------
 2 files changed, 36 insertions(+), 15 deletions(-)

(limited to 'R')

diff --git a/R/mkinds.R b/R/mkinds.R
index 257a17c4..5333ca26 100644
--- a/R/mkinds.R
+++ b/R/mkinds.R
@@ -42,7 +42,6 @@ mkinds <- R6Class("mkinds",
   )
 )
 
-#' @export
 print.mkinds <- function(x, ...) {
   cat("<mkinds> with $title: ",  x$title, "\n")
   cat("Observed compounds $observed: ", paste(x$observed, collapse = ", "), "\n")
diff --git a/R/mkinfit.R b/R/mkinfit.R
index dca71ecf..cca75690 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -263,7 +263,7 @@ mkinfit <- function(mkinmod, observed,
       errparms = rep(3, length(obs_vars))
     }
     if (err_mod == "tc") {
-      errparms <- c(sigma_low = 3, rsd_high = 0.01)
+      errparms <- c(sigma_low = 0.1, rsd_high = 0.1)
     }
     names(errparms) <- errparm_names
   }
@@ -275,13 +275,20 @@ mkinfit <- function(mkinmod, observed,
                                               length.out = n.outtimes))))
 
   # Define log-likelihood function for optimisation, including (back)transformations
-  nlogLik <- function(P, trans = TRUE, OLS = FALSE, update_data = TRUE, ...)
+  nlogLik <- function(P, trans = TRUE, OLS = FALSE, local = FALSE, update_data = TRUE, ...)
   {
     assign("calls", calls + 1, inherits = TRUE) # Increase the model solution counter
+    P.orig <- P
 
     # Trace parameter values if requested
     if(trace_parms) 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
+    }
+
     # Initial states for t0
     if(length(state.ini.optim) > 0) {
       odeini <- c(P[1:length(state.ini.optim)], state.ini.fixed)
@@ -291,7 +298,7 @@ mkinfit <- function(mkinmod, observed,
       names(odeini) <- state.ini.fixed.boxnames
     }
 
-    if (OLS) {
+    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))]
@@ -314,6 +321,9 @@ 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") {
@@ -411,16 +421,17 @@ mkinfit <- function(mkinmod, observed,
 
   # Do the fit and take the time until the hessians are calculated
   fit_time <- system.time({
-    # For constant variance, we do not include sigma in the optimisation, as it
-    # is unnecessary and leads to instability of the fit which are most obvious
-    # when fitting the IORE model
-    if (err_mod == "const") {
-      fit.ols <- nlminb(c(state.ini.optim, transparms.optim),
-                    nlogLik, control = control,
-                    lower = lower, upper = upper, OLS = TRUE, ...)
+    # 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") {
       # Get the maximum likelihood estimate for sigma at the optimum parameter values
-      # This is equivalent to including sigma in the optimisation, but more stable
       data_errmod$residual <- data_errmod$value.observed - data_errmod$value.predicted
       sigma_mle = sqrt(sum(data_errmod$residual^2)/nrow(data_errmod))
 
@@ -428,9 +439,20 @@ mkinfit <- function(mkinmod, observed,
       fit <- fit.ols
       fit$logLik <- - nlogLik(c(fit$par, sigma = sigma_mle), OLS = FALSE)
     } else {
-      fit <- nlminb(c(state.ini.optim, transparms.optim, errparms),
-                    nlogLik, control = control,
-                    lower = lower, upper = upper, ...)
+      # After the OLS stop we fit the error model with fixed degradation model
+      # parameters
+      if (!quiet) message("Optimising the error model")
+      fit.err <- nlminb(errparms, nlogLik, control = control,
+        lower = lower[names(errparms)], 
+        upper = upper[names(errparms)],
+        local = fit.ols$par, ...)
+      errparms.tmp <- fit.err$par
+      if (!quiet) message("Optimising the complete model")
+      parms.start <- c(state.ini.optim, transparms.optim, errparms.tmp)
+      fit <- nlminb(parms.start, nlogLik,
+        lower = lower[names(parms.start)], 
+        upper = upper[names(parms.start)],
+        control = control, ...)
       fit$logLik <- - nlogLik.current
     }
 
-- 
cgit v1.2.1