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

2 files changed, 36 insertions, 15 deletions

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

     }