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