diff options
| -rw-r--r-- | DESCRIPTION | 2 | ||||
| -rw-r--r-- | NEWS.md | 30 | ||||
| -rw-r--r-- | R/mkinfit.R | 117 | ||||
| -rw-r--r-- | man/mkinfit.Rd | 18 | ||||
| -rw-r--r-- | vignettes/FOCUS_L.html | 118 | ||||
| -rw-r--r-- | vignettes/FOCUS_Z.pdf | bin | 212996 -> 214130 bytes | |||
| -rw-r--r-- | vignettes/mkin.pdf | bin | 160326 -> 160326 bytes |
7 files changed, 190 insertions, 95 deletions
diff --git a/DESCRIPTION b/DESCRIPTION index 833acd8e..fb32b807 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -3,7 +3,7 @@ Type: Package Title: Routines for fitting kinetic models with one or more state variables to chemical degradation data Version: 0.9-32 -Date: 2014-07-14 +Date: 2014-07-17 Authors@R: c(person("Johannes", "Ranke", role = c("aut", "cre", "cph"), email = "jranke@uni-bremen.de"), person("Katrin", "Lindenberger", role = "ctb"), @@ -1,11 +1,33 @@ +# CHANGES in mkin VERSION 0.9-32 + +## NEW FEATURES + +- Optimisation method, number of model evaluations and time elapsed during optimisation are given in the summary of mkinfit objects. + +- The maximum number of iterations in the optimisation algorithm can be specified using the argument `maxit.modFit` to the mkinfit function. + +- mkinfit gives a warning when the fit does not converge (does not apply to SANN method). This warning is included in the summary. + +## BUG FIXES + +- Initial values for formation fractions were not set in all cases. + +- No warning was given when the fit did not converge when a method other than the default Levenberg-Marquardt method `Marq` was used. + +## MINOR CHANGES + +- Vignettes were rebuilt to reflect the changes in the summary method. + +- Algorithm `Pseudo` was excluded because it needs user-defined parameter limits which are not supported. + +- Algorithm `Newton` was excluded because of its different way to specify the maximum number of iterations and because it does not appear to provide additional benefits. + # CHANGES in mkin VERSION 0.9-31 ## BUG FIXES - The internal renaming of optimised parameters in Version 0.9-30 led to errors in the determination of the degrees of freedom for the chi2 error level calulations in `mkinerrmin()` used by the summary function. -- Initial values for formation fractions were not set in all cases - # CHANGES in mkin VERSION 0.9-30 ## NEW FEATURES @@ -16,13 +38,13 @@ - The original and the transformed parameters now have different names (e.g. `k_parent` and `log_k_parent`. They also differ in how many they are when we have formation fractions but no pathway to sink. -- The order of some of the information blocks in `print.summary.mkinfit.R()` has been ordered in a more logical way +- The order of some of the information blocks in `print.summary.mkinfit.R()` has been ordered in a more logical way. ## MINOR CHANGES - The vignette FOCUS_Z has been simplified to use formation fractions with turning off the sink, and slightly amended to use the new versions of DT50 values calculated since mkin 0.9-29. -- All vignettes have been rebuilt so they reflect all changes +- All vignettes have been rebuilt so they reflect all changes. - The ChangeLog was renamed to NEWS.md and the entries were converted to markdown syntax compatible with the `tools::news()` function built into R. diff --git a/R/mkinfit.R b/R/mkinfit.R index c6e13b97..d591c42a 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -28,7 +28,8 @@ mkinfit <- function(mkinmod, observed, fixed_initials = names(mkinmod$diffs)[-1],
solution_type = "auto",
method.ode = "lsoda",
- method.modFit = "Marq",
+ method.modFit = c("Marq", "Port", "SANN", "Nelder-Mead", "BFSG", "CG", "L-BFGS-B"),
+ maxit.modFit = "auto",
control.modFit = list(),
transform_rates = TRUE,
transform_fractions = TRUE,
@@ -40,6 +41,16 @@ mkinfit <- function(mkinmod, observed, trace_parms = FALSE,
...)
{
+ # Check optimisation method and set maximum number of iterations if specified
+ method.modFit = match.arg(method.modFit)
+ if (maxit.modFit != "auto") {
+ if (method.modFit == "Marq") control.modFit$maxiter = maxit.modFit
+ if (method.modFit == "Port") control.modFit$iter.max = maxit.modFit
+ if (method.modFit %in% c("SANN", "Nelder-Mead", "BFGS", "CG", "L-BFGS-B")) {
+ control.modFit$maxit = maxit.modFit
+ }
+ }
+
# Get the names of the state variables in the model
mod_vars <- names(mkinmod$diffs)
@@ -281,48 +292,76 @@ mkinfit <- function(mkinmod, observed, upper[other_fraction_parms] <- 1
}
- fit <- modFit(cost, c(state.ini.optim, transparms.optim),
- method = method.modFit, control = control.modFit,
- lower = lower, upper = upper, ...)
-
- # Reiterate the fit until convergence of the variance components (IRLS)
- # if requested by the user
- weight.ini = weight
- if (!is.null(err)) weight.ini = "manual"
-
- if (!is.null(reweight.method)) {
- if (reweight.method != "obs") stop("Only reweighting method 'obs' is implemented")
- if(!quiet) {
- cat("IRLS based on variance estimates for each observed variable\n")
- }
- if (!quiet) {
- cat("Initial variance estimates are:\n")
- print(signif(fit$var_ms_unweighted, 8))
- }
- reweight.diff = 1
- n.iter <- 0
- if (!is.null(err)) observed$err.ini <- observed[[err]]
- err = "err.irls"
- while (reweight.diff > reweight.tol & n.iter < reweight.max.iter) {
- n.iter <- n.iter + 1
- sigma.old <- sqrt(fit$var_ms_unweighted)
- observed[err] <- sqrt(fit$var_ms_unweighted)[as.character(observed$name)]
- fit <- modFit(cost, fit$par, method = method.modFit,
- control = control.modFit, lower = lower, upper = upper, ...)
- reweight.diff = sum((sqrt(fit$var_ms_unweighted) - sigma.old)^2)
+ # Do the fit and take the time
+ fit_time <- system.time({
+ fit <- modFit(cost, c(state.ini.optim, transparms.optim),
+ method = method.modFit, control = control.modFit,
+ lower = lower, upper = upper, ...)
+
+ # Reiterate the fit until convergence of the variance components (IRLS)
+ # if requested by the user
+ weight.ini = weight
+ if (!is.null(err)) weight.ini = "manual"
+
+ if (!is.null(reweight.method)) {
+ if (reweight.method != "obs") stop("Only reweighting method 'obs' is implemented")
+ if(!quiet) {
+ cat("IRLS based on variance estimates for each observed variable\n")
+ }
if (!quiet) {
- cat("Iteration", n.iter, "yields variance estimates:\n")
+ cat("Initial variance estimates are:\n")
print(signif(fit$var_ms_unweighted, 8))
- cat("Sum of squared differences to last variance estimates:",
- signif(reweight.diff, 2), "\n")
+ }
+ reweight.diff = 1
+ n.iter <- 0
+ if (!is.null(err)) observed$err.ini <- observed[[err]]
+ err = "err.irls"
+ while (reweight.diff > reweight.tol & n.iter < reweight.max.iter) {
+ n.iter <- n.iter + 1
+ sigma.old <- sqrt(fit$var_ms_unweighted)
+ observed[err] <- sqrt(fit$var_ms_unweighted)[as.character(observed$name)]
+ fit <- modFit(cost, fit$par, method = method.modFit,
+ control = control.modFit, lower = lower, upper = upper, ...)
+ reweight.diff = sum((sqrt(fit$var_ms_unweighted) - sigma.old)^2)
+ if (!quiet) {
+ cat("Iteration", n.iter, "yields variance estimates:\n")
+ print(signif(fit$var_ms_unweighted, 8))
+ cat("Sum of squared differences to last variance estimates:",
+ signif(reweight.diff, 2), "\n")
+ }
}
}
+ })
+
+ # Check for convergence
+ if (method.modFit == "Marq") {
+ if (!fit$info %in% c(1, 2, 3)) {
+ fit$warning = paste0("Optimisation by method ", method.modFit,
+ " did not converge.\n",
+ "The message returned by nls.lm is:\n",
+ fit$message)
+ warning(fit$warning)
+ }
+ }
+ if (method.modFit %in% c("Port", "SANN", "Nelder-Mead", "BFGS", "CG", "L-BFGS-B")) {
+ if (fit$convergence != 0) {
+ fit$warning = paste0("Optimisation by method ", method.modFit,
+ " did not converge.\n",
+ "Convergence code is ", fit$convergence,
+ ifelse(is.null(fit$message), "",
+ paste0("\nMessage is ", fit$message)))
+ warning(fit$warning)
+ }
}
# We need to return some more data for summary and plotting
fit$solution_type <- solution_type
fit$transform_rates <- transform_rates
fit$transform_fractions <- transform_fractions
+ fit$method.modFit <- method.modFit
+ fit$maxit.modFit <- maxit.modFit
+ fit$calls <- calls
+ fit$time <- fit_time
# We also need the model for summary and plotting
fit$mkinmod <- mkinmod
@@ -449,6 +488,8 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, date.fit = object$date,
date.summary = date(),
solution_type = object$solution_type,
+ method.modFit = object$method.modFit,
+ warning = object$warning,
use_of_ff = object$mkinmod$use_of_ff,
weight.ini = object$weight.ini,
reweight.method = object$reweight.method,
@@ -461,6 +502,8 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, cov.scaled = covar * resvar,
info = object$info,
niter = object$iterations,
+ calls = object$calls,
+ time = object$time,
stopmess = message,
par = param,
bpar = bparam)
@@ -491,13 +534,17 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), . cat("Date of fit: ", x$date.fit, "\n")
cat("Date of summary:", x$date.summary, "\n")
+ if (!is.null(x$warning)) cat("\n\nWarning:", x$warning, "\n\n")
+
cat("\nEquations:\n")
print(noquote(as.character(x[["diffs"]])))
df <- x$df
rdf <- df[2]
- cat("\nMethod used for solution of differential equation system:\n")
- cat(x$solution_type, "\n")
+ cat("\nModel predictions using solution type", x$solution_type, "\n")
+
+ cat("\nFitted with method", x$method.modFit,
+ "using", x$calls, "model solutions performed in", x$time[["elapsed"]], "s\n")
cat("\nWeighting:", x$weight.ini)
if(!is.null(x$reweight.method)) cat(" then iterative reweighting method",
diff --git a/man/mkinfit.Rd b/man/mkinfit.Rd index bd7f73b7..c99e146c 100644 --- a/man/mkinfit.Rd +++ b/man/mkinfit.Rd @@ -21,7 +21,8 @@ mkinfit(mkinmod, observed, fixed_parms = NULL, fixed_initials = names(mkinmod$diffs)[-1], solution_type = "auto", method.ode = "lsoda", - method.modFit = "Marq", + method.modFit = c("Marq", "Port", "SANN", "Nelder-Mead", "BFSG", "CG", "L-BFGS-B"), + maxit.modFit = "auto", control.modFit = list(), transform_rates = TRUE, transform_fractions = TRUE, @@ -102,12 +103,21 @@ mkinfit(mkinmod, observed, recommended as it is less prone to get trapped in local minima and depends less on starting values for parameters. However, it needs more iterations. - When using "Pseudo", "upper" and "lower" need to be specified for the - transformed parameters. + The "Pseudo" algorithm is not included because it needs finite parameter bounds + which are currently not supported. + + The "Newton" algorithm is not included because its number of iterations + can not be controlled by \code{control.modFit} and it does not appear + to provide advantages over the other algorithms. + } + \item{maxit.modFit}{ + Maximum number of iterations in the optimisation. If not "auto", this will + be passed to the method called by \code{\link{modFit}}, overriding + what may be specified in the next argument \code{control.modFit}. } \item{control.modFit}{ Additional arguments passed to the optimisation method used by - \code{\link{modFit}}. + \code{\link{modFit}}. } \item{transform_rates}{ Boolean specifying if kinetic rate constants should be transformed in the diff --git a/vignettes/FOCUS_L.html b/vignettes/FOCUS_L.html index bb02ec3e..85fadbfe 100644 --- a/vignettes/FOCUS_L.html +++ b/vignettes/FOCUS_L.html @@ -193,7 +193,13 @@ hr { report, p. 284</p> <pre><code class="r">library("mkin") -FOCUS_2006_L1 = data.frame( +</code></pre> + +<pre><code>## Loading required package: minpack.lm +## Loading required package: rootSolve +</code></pre> + +<pre><code class="r">FOCUS_2006_L1 = data.frame( t = rep(c(0, 1, 2, 3, 5, 7, 14, 21, 30), each = 2), parent = c(88.3, 91.4, 85.6, 84.5, 78.9, 77.6, 72.0, 71.9, 50.3, 59.4, 47.0, 45.1, @@ -223,16 +229,17 @@ FOCUS report.</p> summary(m.L1.SFO) </code></pre> -<pre><code>## mkin version: 0.9.31 +<pre><code>## mkin version: 0.9.32 ## R version: 3.1.1 -## Date of fit: Mon Jul 14 20:32:20 2014 -## Date of summary: Mon Jul 14 20:32:20 2014 +## Date of fit: Thu Jul 17 12:37:41 2014 +## Date of summary: Thu Jul 17 12:37:41 2014 ## ## Equations: ## [1] d_parent = - k_parent_sink * parent ## -## Method used for solution of differential equation system: -## analytical +## Model predictions using solution type analytical +## +## Fitted with method Marq using 14 model solutions performed in 0.087 s ## ## Weighting: none ## @@ -325,16 +332,17 @@ is checked.</p> summary(m.L1.FOMC, data = FALSE) </code></pre> -<pre><code>## mkin version: 0.9.31 +<pre><code>## mkin version: 0.9.32 ## R version: 3.1.1 -## Date of fit: Mon Jul 14 20:32:20 2014 -## Date of summary: Mon Jul 14 20:32:20 2014 +## Date of fit: Thu Jul 17 12:37:42 2014 +## Date of summary: Thu Jul 17 12:37:42 2014 ## ## Equations: ## [1] d_parent = - (alpha/beta) * ((time/beta) + 1)^-1 * parent ## -## Method used for solution of differential equation system: -## analytical +## Model predictions using solution type analytical +## +## Fitted with method Marq using 45 model solutions performed in 0.266 s ## ## Weighting: none ## @@ -417,16 +425,17 @@ FOCUS_2006_L2_mkin <- mkin_wide_to_long(FOCUS_2006_L2) summary(m.L2.SFO) </code></pre> -<pre><code>## mkin version: 0.9.31 +<pre><code>## mkin version: 0.9.32 ## R version: 3.1.1 -## Date of fit: Mon Jul 14 20:32:20 2014 -## Date of summary: Mon Jul 14 20:32:20 2014 +## Date of fit: Thu Jul 17 12:37:42 2014 +## Date of summary: Thu Jul 17 12:37:42 2014 ## ## Equations: ## [1] d_parent = - k_parent_sink * parent ## -## Method used for solution of differential equation system: -## analytical +## Model predictions using solution type analytical +## +## Fitted with method Marq using 32 model solutions performed in 0.357 s ## ## Weighting: none ## @@ -526,16 +535,17 @@ mkinresplot(m.L2.FOMC) <pre><code class="r">summary(m.L2.FOMC, data = FALSE) </code></pre> -<pre><code>## mkin version: 0.9.31 +<pre><code>## mkin version: 0.9.32 ## R version: 3.1.1 -## Date of fit: Mon Jul 14 20:32:21 2014 -## Date of summary: Mon Jul 14 20:32:21 2014 +## Date of fit: Thu Jul 17 12:37:43 2014 +## Date of summary: Thu Jul 17 12:37:43 2014 ## ## Equations: ## [1] d_parent = - (alpha/beta) * ((time/beta) + 1)^-1 * parent ## -## Method used for solution of differential equation system: -## analytical +## Model predictions using solution type analytical +## +## Fitted with method Marq using 39 model solutions performed in 0.235 s ## ## Weighting: none ## @@ -611,16 +621,17 @@ plot(m.L2.DFOP) <pre><code class="r">summary(m.L2.DFOP, data = FALSE) </code></pre> -<pre><code>## mkin version: 0.9.31 +<pre><code>## mkin version: 0.9.32 ## R version: 3.1.1 -## Date of fit: Mon Jul 14 20:32:23 2014 -## Date of summary: Mon Jul 14 20:32:23 2014 +## Date of fit: Thu Jul 17 12:37:44 2014 +## Date of summary: Thu Jul 17 12:37:44 2014 ## ## Equations: ## [1] d_parent = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) * parent ## -## Method used for solution of differential equation system: -## analytical +## Model predictions using solution type analytical +## +## Fitted with method Marq using 54 model solutions performed in 0.423 s ## ## Weighting: none ## @@ -697,16 +708,17 @@ plot(m.L3.SFO) <pre><code class="r">summary(m.L3.SFO) </code></pre> -<pre><code>## mkin version: 0.9.31 +<pre><code>## mkin version: 0.9.32 ## R version: 3.1.1 -## Date of fit: Mon Jul 14 20:32:23 2014 -## Date of summary: Mon Jul 14 20:32:23 2014 +## Date of fit: Thu Jul 17 12:37:45 2014 +## Date of summary: Thu Jul 17 12:37:45 2014 ## ## Equations: ## [1] d_parent = - k_parent_sink * parent ## -## Method used for solution of differential equation system: -## analytical +## Model predictions using solution type analytical +## +## Fitted with method Marq using 44 model solutions performed in 0.241 s ## ## Weighting: none ## @@ -782,16 +794,17 @@ plot(m.L3.FOMC) <pre><code class="r">summary(m.L3.FOMC, data = FALSE) </code></pre> -<pre><code>## mkin version: 0.9.31 +<pre><code>## mkin version: 0.9.32 ## R version: 3.1.1 -## Date of fit: Mon Jul 14 20:32:24 2014 -## Date of summary: Mon Jul 14 20:32:24 2014 +## Date of fit: Thu Jul 17 12:37:45 2014 +## Date of summary: Thu Jul 17 12:37:45 2014 ## ## Equations: ## [1] d_parent = - (alpha/beta) * ((time/beta) + 1)^-1 * parent ## -## Method used for solution of differential equation system: -## analytical +## Model predictions using solution type analytical +## +## Fitted with method Marq using 26 model solutions performed in 0.208 s ## ## Weighting: none ## @@ -854,16 +867,17 @@ plot(m.L3.DFOP) <pre><code class="r">summary(m.L3.DFOP, data = FALSE) </code></pre> -<pre><code>## mkin version: 0.9.31 +<pre><code>## mkin version: 0.9.32 ## R version: 3.1.1 -## Date of fit: Mon Jul 14 20:32:24 2014 -## Date of summary: Mon Jul 14 20:32:24 2014 +## Date of fit: Thu Jul 17 12:37:46 2014 +## Date of summary: Thu Jul 17 12:37:46 2014 ## ## Equations: ## [1] d_parent = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) * parent ## -## Method used for solution of differential equation system: -## analytical +## Model predictions using solution type analytical +## +## Fitted with method Marq using 37 model solutions performed in 0.338 s ## ## Weighting: none ## @@ -944,16 +958,17 @@ plot(m.L4.SFO) <pre><code class="r">summary(m.L4.SFO, data = FALSE) </code></pre> -<pre><code>## mkin version: 0.9.31 +<pre><code>## mkin version: 0.9.32 ## R version: 3.1.1 -## Date of fit: Mon Jul 14 20:32:27 2014 -## Date of summary: Mon Jul 14 20:32:27 2014 +## Date of fit: Thu Jul 17 12:37:46 2014 +## Date of summary: Thu Jul 17 12:37:46 2014 ## ## Equations: ## [1] d_parent = - k_parent_sink * parent ## -## Method used for solution of differential equation system: -## analytical +## Model predictions using solution type analytical +## +## Fitted with method Marq using 20 model solutions performed in 0.127 s ## ## Weighting: none ## @@ -1018,16 +1033,17 @@ plot(m.L4.FOMC) <pre><code class="r">summary(m.L4.FOMC, data = FALSE) </code></pre> -<pre><code>## mkin version: 0.9.31 +<pre><code>## mkin version: 0.9.32 ## R version: 3.1.1 -## Date of fit: Mon Jul 14 20:32:28 2014 -## Date of summary: Mon Jul 14 20:32:28 2014 +## Date of fit: Thu Jul 17 12:37:46 2014 +## Date of summary: Thu Jul 17 12:37:46 2014 ## ## Equations: ## [1] d_parent = - (alpha/beta) * ((time/beta) + 1)^-1 * parent ## -## Method used for solution of differential equation system: -## analytical +## Model predictions using solution type analytical +## +## Fitted with method Marq using 53 model solutions performed in 0.355 s ## ## Weighting: none ## diff --git a/vignettes/FOCUS_Z.pdf b/vignettes/FOCUS_Z.pdf Binary files differindex 559504cc..43f3e2e2 100644 --- a/vignettes/FOCUS_Z.pdf +++ b/vignettes/FOCUS_Z.pdf diff --git a/vignettes/mkin.pdf b/vignettes/mkin.pdf Binary files differindex f06d38f9..9cf1b3e5 100644 --- a/vignettes/mkin.pdf +++ b/vignettes/mkin.pdf |