diff options
| -rw-r--r-- | NEWS.md | 6 | ||||
| -rw-r--r-- | R/IORE.solution.R | 4 | ||||
| -rw-r--r-- | R/endpoints.R | 28 | ||||
| -rw-r--r-- | R/mkinerrmin.R | 198 | ||||
| -rw-r--r-- | R/mkinfit.R | 13 | ||||
| -rw-r--r-- | R/mkinmod.R | 54 | ||||
| -rw-r--r-- | R/mkinparplot.R | 8 | ||||
| -rw-r--r-- | R/mkinpredict.R | 10 | ||||
| -rw-r--r-- | R/mkinsub.R | 23 | ||||
| -rw-r--r-- | R/transform_odeparms.R | 19 | ||||
| -rw-r--r-- | README.md | 3 | ||||
| -rw-r--r-- | man/FOMC.solution.Rd | 2 | ||||
| -rw-r--r-- | man/IORE.solution.Rd | 34 | ||||
| -rw-r--r-- | man/mkinfit.Rd | 6 | ||||
| -rw-r--r-- | man/mkinmod.Rd | 27 | ||||
| -rw-r--r-- | man/mkinpredict.Rd | 19 | ||||
| -rw-r--r-- | man/mkinsub.Rd | 44 | ||||
| -rw-r--r-- | vignettes/FOCUS_Z.pdf | bin | 220189 -> 398825 bytes | |||
| -rw-r--r-- | vignettes/mkin.pdf | bin | 160333 -> 295336 bytes |
19 files changed, 360 insertions, 138 deletions
@@ -2,8 +2,11 @@ ## NEW FEATURES -- Switch to using the Port algorithm (using a model/trust region approach) per default. While needing more iterations than the Levenberg-Marquardt algorithm previously used per default, it is less sensitive to starting parameters. +- Add the convenience function `mkinsub()` for creating the lists used in `mkinmod()` + +- Add the possibility to fit indeterminate order rate equation (IORE) models using an analytical solution (parent only) or a numeric solution. Paths from IORE compounds to metabolites are supported when using of formation fractions (use_of_ff = 'max'). Note that the numerical solution (method.ode = 'deSolve') of the IORE differential equations sometimes fails due to numerical problems. +- Switch to using the Port algorithm (using a model/trust region approach) per default. While needing more iterations than the Levenberg-Marquardt algorithm previously used per default, it is less sensitive to starting parameters. ## MINOR CHANGES - The formatting of differential equations in the summary was further improved @@ -79,7 +82,6 @@ ## 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. - # CHANGES in mkin VERSION 0.9-30 ## NEW FEATURES diff --git a/R/IORE.solution.R b/R/IORE.solution.R new file mode 100644 index 00000000..9546ce56 --- /dev/null +++ b/R/IORE.solution.R @@ -0,0 +1,4 @@ +IORE.solution <- function(t, parent.0, k.iore, N)
+{
+ parent = (parent.0^(1 - N) - (1 - N) * k.iore * t)^(1/(1 - N))
+}
diff --git a/R/endpoints.R b/R/endpoints.R index 676831a0..ae9aa3ec 100644 --- a/R/endpoints.R +++ b/R/endpoints.R @@ -1,11 +1,11 @@ endpoints <- function(fit) {
# Calculate dissipation times DT50 and DT90 and formation
# fractions as well as SFORB eigenvalues from optimised parameters
- # Additional DT50 values are calculated from the FOMC DT90 and k1 and k2 from HS and DFOP,
- # as well as from Eigenvalues b1 and b2 of any SFORB models
+ # Additional DT50 values are calculated from the FOMC DT90 and k1 and k2 from
+ # HS and DFOP, as well as from Eigenvalues b1 and b2 of any SFORB models
ep <- list()
obs_vars <- fit$obs_vars
- parms.all <- fit$bparms.ode
+ parms.all <- c(fit$bparms.optim, fit$bparms.fixed)
ep$ff <- vector()
ep$SFORB <- vector()
ep$distimes <- data.frame(DT50 = rep(NA, length(obs_vars)),
@@ -49,6 +49,28 @@ endpoints <- function(fit) { DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011
ep$distimes[obs_var, c("DT50back")] = DT50_back
}
+ if (type == "IORE") {
+ k_names = grep(paste("k.iore", obs_var, sep="_"), names(parms.all), value=TRUE)
+ k_tot = sum(parms.all[k_names])
+ # From the NAFTA kinetics guidance, p. 5
+ n = parms.all[paste("N", obs_var, sep = "_")]
+ k = k_tot
+ # Use the initial concentration of the parent compound
+ source_name = fit$mkinmod$map[[1]][[1]]
+ c0 = parms.all[paste(source_name, "0", sep = "_")]
+ alpha = 1 / (n - 1)
+ beta = (c0^(1 - n))/(k * (n - 1))
+ DT50 = beta * (2^(1/alpha) - 1)
+ DT90 = beta * (10^(1/alpha) - 1)
+ DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011
+ ep$distimes[obs_var, c("DT50back")] = DT50_back
+ if (fit$mkinmod$use_of_ff == "min") {
+ for (k_name in k_names)
+ {
+ ep$ff[[sub("k_", "", k_name)]] = parms.all[[k_name]] / k_tot
+ }
+ }
+ }
if (type == "DFOP") {
k1 = parms.all["k1"]
k2 = parms.all["k2"]
diff --git a/R/mkinerrmin.R b/R/mkinerrmin.R index 2697d0a0..b28235a6 100644 --- a/R/mkinerrmin.R +++ b/R/mkinerrmin.R @@ -1,96 +1,102 @@ -# Copyright (C) 2010-2014 Johannes Ranke
-# Contact: jranke@uni-bremen.de
-
-# This file is part of the R package mkin
-
-# mkin is free software: you can redistribute it and/or modify it under the
-# terms of the GNU General Public License as published by the Free Software
-# Foundation, either version 3 of the License, or (at your option) any later
-# version.
-
-# This program is distributed in the hope that it will be useful, but WITHOUT
-# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
-# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-# details.
-
-# You should have received a copy of the GNU General Public License along with
-# this program. If not, see <http://www.gnu.org/licenses/>
-if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "value_mean"))
-
-mkinerrmin <- function(fit, alpha = 0.05)
-{
- parms.optim <- fit$par
-
- kinerrmin <- function(errdata, n.parms) {
- means.mean <- mean(errdata$value_mean, na.rm = TRUE)
- df = length(errdata$value_mean) - n.parms
-
- err.min <- sqrt((1 / qchisq(1 - alpha, df)) *
- sum((errdata$value_mean - errdata$value_pred)^2)/(means.mean^2))
-
- return(list(err.min = err.min, n.optim = n.parms, df = df))
- }
-
- means <- aggregate(value ~ time + name, data = fit$observed, mean, na.rm=TRUE)
- errdata <- merge(means, fit$predicted, by = c("time", "name"),
- suffixes = c("_mean", "_pred"))
- errdata <- errdata[order(errdata$time, errdata$name), ]
-
- # Remove values at time zero for variables whose value for state.ini is fixed,
- # as these will not have any effect in the optimization and should therefore not
- # be counted as degrees of freedom.
- fixed_initials = gsub("_0$", "", rownames(subset(fit$fixed, type = "state")))
- errdata <- subset(errdata, !(time == 0 & name %in% fixed_initials))
-
- n.optim.overall <- length(parms.optim)
-
- errmin.overall <- kinerrmin(errdata, n.optim.overall)
- errmin <- data.frame(err.min = errmin.overall$err.min,
- n.optim = errmin.overall$n.optim, df = errmin.overall$df)
- rownames(errmin) <- "All data"
-
- # The degrees of freedom are counted according to FOCUS kinetics (2011, p. 164)
- for (obs_var in fit$obs_vars)
- {
- errdata.var <- subset(errdata, name == obs_var)
-
- # Check if initial value is optimised
- n.initials.optim <- length(grep(paste(obs_var, ".*", "_0", sep=""), names(parms.optim)))
-
- # Rate constants are attributed to the source variable
- n.k.optim <- length(grep(paste("^k", obs_var, sep="_"), names(parms.optim)))
- n.k.optim <- n.k.optim + length(grep(paste("^log_k", obs_var, sep="_"),
- names(parms.optim)))
-
- n.ff.optim <- 0
- # Formation fractions are attributed to the target variable, so look
- # for source compartments with formation fractions
- for (source_var in fit$obs_vars) {
- n.ff.source = length(grep(paste("^f", source_var, sep = "_"),
- names(parms.optim)))
- n.paths.source = length(fit$mkinmod$spec[[source_var]]$to)
- for (target_var in fit$mkinmod$spec[[source_var]]$to) {
- if (obs_var == target_var) {
- n.ff.optim <- n.ff.optim + n.ff.source/n.paths.source
- }
- }
- }
-
- n.optim <- n.k.optim + n.initials.optim + n.ff.optim
-
- # FOMC, DFOP and HS parameters are only counted if we are looking at the
- # first variable in the model which is always the source variable
- if (obs_var == fit$obs_vars[[1]]) {
- special_parms = c("alpha", "log_alpha", "beta", "log_beta",
- "k1", "log_k1", "k2", "log_k2",
- "g", "g_ilr", "tb", "log_tb")
- n.optim <- n.optim + length(intersect(special_parms, names(parms.optim)))
- }
-
- # Calculate and add a line to the results
- errmin.tmp <- kinerrmin(errdata.var, n.optim)
- errmin[obs_var, c("err.min", "n.optim", "df")] <- errmin.tmp
- }
-
- return(errmin)
-}
+# Copyright (C) 2010-2014 Johannes Ranke +# Contact: jranke@uni-bremen.de + +# This file is part of the R package mkin + +# mkin is free software: you can redistribute it and/or modify it under the +# terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. + +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. + +# You should have received a copy of the GNU General Public License along with +# this program. If not, see <http://www.gnu.org/licenses/> +if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "value_mean")) + +mkinerrmin <- function(fit, alpha = 0.05) +{ + parms.optim <- fit$par + + kinerrmin <- function(errdata, n.parms) { + means.mean <- mean(errdata$value_mean, na.rm = TRUE) + df = length(errdata$value_mean) - n.parms + + err.min <- sqrt((1 / qchisq(1 - alpha, df)) * + sum((errdata$value_mean - errdata$value_pred)^2)/(means.mean^2)) + + return(list(err.min = err.min, n.optim = n.parms, df = df)) + } + + means <- aggregate(value ~ time + name, data = fit$observed, mean, na.rm=TRUE) + errdata <- merge(means, fit$predicted, by = c("time", "name"), + suffixes = c("_mean", "_pred")) + errdata <- errdata[order(errdata$time, errdata$name), ] + + # Remove values at time zero for variables whose value for state.ini is fixed, + # as these will not have any effect in the optimization and should therefore not + # be counted as degrees of freedom. + fixed_initials = gsub("_0$", "", rownames(subset(fit$fixed, type = "state"))) + errdata <- subset(errdata, !(time == 0 & name %in% fixed_initials)) + + n.optim.overall <- length(parms.optim) + + errmin.overall <- kinerrmin(errdata, n.optim.overall) + errmin <- data.frame(err.min = errmin.overall$err.min, + n.optim = errmin.overall$n.optim, df = errmin.overall$df) + rownames(errmin) <- "All data" + + # The degrees of freedom are counted according to FOCUS kinetics (2011, p. 164) + for (obs_var in fit$obs_vars) + { + errdata.var <- subset(errdata, name == obs_var) + + # Check if initial value is optimised + n.initials.optim <- length(grep(paste(obs_var, ".*", "_0", sep=""), names(parms.optim))) + + # Rate constants and IORE exponents are attributed to the source variable + n.k.optim <- length(grep(paste("^k", obs_var, sep="_"), names(parms.optim))) + n.k.optim <- n.k.optim + length(grep(paste("^log_k", obs_var, sep="_"), + names(parms.optim))) + n.k.iore.optim <- length(grep(paste("^k.iore", obs_var, sep="_"), names(parms.optim))) + n.k.iore.optim <- n.k.iore.optim + length(grep(paste("^log_k.iore", obs_var, + sep = "_"), + names(parms.optim))) + + n.N.optim <- length(grep(paste("^N", obs_var, sep="_"), names(parms.optim))) + + n.ff.optim <- 0 + # Formation fractions are attributed to the target variable, so look + # for source compartments with formation fractions + for (source_var in fit$obs_vars) { + n.ff.source = length(grep(paste("^f", source_var, sep = "_"), + names(parms.optim))) + n.paths.source = length(fit$mkinmod$spec[[source_var]]$to) + for (target_var in fit$mkinmod$spec[[source_var]]$to) { + if (obs_var == target_var) { + n.ff.optim <- n.ff.optim + n.ff.source/n.paths.source + } + } + } + + n.optim <- sum(n.initials.optim, n.k.optim, n.k.iore.optim, n.N.optim, n.ff.optim) + + # FOMC, DFOP and HS parameters are only counted if we are looking at the + # first variable in the model which is always the source variable + if (obs_var == fit$obs_vars[[1]]) { + special_parms = c("alpha", "log_alpha", "beta", "log_beta", + "k1", "log_k1", "k2", "log_k2", + "g", "g_ilr", "tb", "log_tb") + n.optim <- n.optim + length(intersect(special_parms, names(parms.optim))) + } + + # Calculate and add a line to the dataframe holding the results + errmin.tmp <- kinerrmin(errdata.var, n.optim) + errmin[obs_var, c("err.min", "n.optim", "df")] <- errmin.tmp + } + + return(errmin) +} diff --git a/R/mkinfit.R b/R/mkinfit.R index 59598631..7a465c2b 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -43,7 +43,7 @@ mkinfit <- function(mkinmod, observed, {
# Check mkinmod and generate a model for the variable whith the highest value
# if a suitable string is given
- parent_models_available = c("SFO", "FOMC", "DFOP", "HS", "SFORB")
+ parent_models_available = c("SFO", "FOMC", "DFOP", "HS", "SFORB", "IORE")
if (class(mkinmod) != "mkinmod") {
presumed_parent_name = observed[which.max(observed$value), "name"]
if (mkinmod[[1]] %in% parent_models_available) {
@@ -60,7 +60,10 @@ mkinfit <- function(mkinmod, observed, 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 == "Port") {
+ control.modFit$iter.max = maxit.modFit
+ control.modFit$eval.max = maxit.modFit
+ }
if (method.modFit %in% c("SANN", "Nelder-Mead", "BFGS", "CG", "L-BFGS-B")) {
control.modFit$maxit = maxit.modFit
}
@@ -116,13 +119,15 @@ mkinfit <- function(mkinmod, observed, defaultpar.names <- setdiff(mkinmod$parms, names(parms.ini))
for (parmname in defaultpar.names) {
# Default values for rate constants, depending on the parameterisation
- if (substr(parmname, 1, 2) == "k_") {
+ if (grepl("^k", parmname)) {
parms.ini[parmname] = 0.1 + k_salt
k_salt = k_salt + 1e-4
}
# Default values for rate constants for reversible binding
if (grepl("free_bound$", parmname)) parms.ini[parmname] = 0.1
if (grepl("bound_free$", parmname)) parms.ini[parmname] = 0.02
+ # Default values for IORE exponents
+ if (grepl("^N", parmname)) parms.ini[parmname] = 1
# Default values for the FOMC, DFOP and HS models
if (parmname == "alpha") parms.ini[parmname] = 1
if (parmname == "beta") parms.ini[parmname] = 10
@@ -306,6 +311,8 @@ mkinfit <- function(mkinmod, observed, if (!transform_rates) {
index_k <- grep("^k_", names(lower))
lower[index_k] <- 0
+ index_k.iore <- grep("^k.iore_", names(lower))
+ lower[index_k.iore] <- 0
other_rate_parms <- intersect(c("alpha", "beta", "k1", "k2", "tb"), names(lower))
lower[other_rate_parms] <- 0
}
diff --git a/R/mkinmod.R b/R/mkinmod.R index eb290719..8b86303f 100644 --- a/R/mkinmod.R +++ b/R/mkinmod.R @@ -1,4 +1,4 @@ -# Copyright (C) 2010-2013 Johannes Ranke {{{
+# Copyright (C) 2010-2014 Johannes Ranke {{{
# Contact: jranke@uni-bremen.de
# This file is part of the R package mkin
@@ -33,38 +33,39 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL) stop("The use of formation fractions 'use_of_ff' can only be 'min' or 'max'")
# The returned model will be a list of character vectors, containing {{{
- # differential equations, parameter names and a mapping from model variables
- # to observed variables. If possible, a matrix representation of the
- # differential equations is included
+ # differential equations (if supported), parameter names and a mapping from
+ # model variables to observed variables. If possible, a matrix representation
+ # of the differential equations is included
parms <- vector()
- diffs <- vector()
- map <- list()
# }}}
- # Do not return a coefficient matrix mat when FOMC, DFOP or HS are used for
- # the parent compound {{{
- if(spec[[1]]$type %in% c("FOMC", "DFOP", "HS")) {
+ # Do not return a coefficient matrix mat when FOMC, IORE, DFOP or HS is used for the parent {{{
+ if(spec[[1]]$type %in% c("FOMC", "IORE", "DFOP", "HS")) {
mat = FALSE
} else mat = TRUE
#}}}
# Establish a list of differential equations as well as a map from observed {{{
# compartments to differential equations
+ diffs <- vector()
+ map <- list()
for (varname in obs_vars)
{
# Check the type component of the compartment specification {{{
if(is.null(spec[[varname]]$type)) stop(
"Every part of the model specification must be a list containing a type component")
- if(!spec[[varname]]$type %in% c("SFO", "FOMC", "DFOP", "HS", "SFORB")) stop(
- "Available types are SFO, FOMC, DFOP, HS and SFORB only")
+ if(!spec[[varname]]$type %in% c("SFO", "FOMC", "IORE", "DFOP", "HS", "SFORB")) stop(
+ "Available types are SFO, FOMC, IORE, DFOP, HS and SFORB only")
if(spec[[varname]]$type %in% c("FOMC", "DFOP", "HS") & match(varname, obs_vars) != 1) {
stop(paste("Types FOMC, DFOP and HS are only implemented for the first compartment,",
"which is assumed to be the source compartment"))
- } #}}}
+ }
+ #}}}
# New (sub)compartments (boxes) needed for the model type {{{
new_boxes <- switch(spec[[varname]]$type,
SFO = varname,
FOMC = varname,
+ IORE = varname,
DFOP = varname,
HS = varname,
SFORB = paste(varname, c("free", "bound"), sep = "_")
@@ -85,20 +86,40 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL) # Turn on sink if this is not explicitly excluded by the user by
# specifying sink=FALSE
if(is.null(spec[[varname]]$sink)) spec[[varname]]$sink <- TRUE
- if(spec[[varname]]$type %in% c("SFO", "SFORB")) { # {{{ Add SFO or SFORB decline term
+ if(spec[[varname]]$type %in% c("SFO", "IORE", "SFORB")) { # {{{ Add decline term
if (use_of_ff == "min") { # Minimum use of formation fractions
if(spec[[varname]]$sink) {
- # If sink is required, add first-order sink term
+ if(spec[[varname]]$type == "IORE" && length(spec[[varname]]$to) > 0) {
+ stop("Transformation reactions from compounds modelled with IORE\n",
+ "are only supported with formation fractions (use_of_ff = 'max')")
+ }
+ # If sink is required, add first-order/IORE sink term
k_compound_sink <- paste("k", box_1, "sink", sep = "_")
+ if(spec[[varname]]$type == "IORE") {
+ k_compound_sink <- paste("k.iore", box_1, "sink", sep = "_")
+ }
parms <- c(parms, k_compound_sink)
decline_term <- paste(k_compound_sink, "*", box_1)
+ if(spec[[varname]]$type == "IORE") {
+ N <- paste("N", box_1, sep = "_")
+ parms <- c(parms, N)
+ decline_term <- paste0(decline_term, "^", N)
+ }
} else { # otherwise no decline term needed here
decline_term = "0"
}
} else {
k_compound <- paste("k", box_1, sep = "_")
+ if(spec[[varname]]$type == "IORE") {
+ k_compound <- paste("k.iore", box_1, sep = "_")
+ }
parms <- c(parms, k_compound)
decline_term <- paste(k_compound, "*", box_1)
+ if(spec[[varname]]$type == "IORE") {
+ N <- paste("N", box_1, sep = "_")
+ parms <- c(parms, N)
+ decline_term <- paste0(decline_term, "^", N)
+ }
}
} #}}}
if(spec[[varname]]$type == "FOMC") { # {{{ Add FOMC decline term
@@ -157,9 +178,10 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL) for (target in to) {
target_box <- switch(spec[[target]]$type,
SFO = target,
+ IORE = target,
SFORB = paste(target, "free", sep = "_"))
- if (use_of_ff == "min" && spec[[varname]]$type %in% c("SFO",
- "SFORB")) {
+ if (use_of_ff == "min" && spec[[varname]]$type %in% c("SFO", "SFORB"))
+ {
k_from_to <- paste("k", origin_box, target_box, sep = "_")
parms <- c(parms, k_from_to)
diffs[[origin_box]] <- paste(diffs[[origin_box]], "-",
diff --git a/R/mkinparplot.R b/R/mkinparplot.R index 683a8b24..4abecab9 100644 --- a/R/mkinparplot.R +++ b/R/mkinparplot.R @@ -19,11 +19,13 @@ mkinparplot <- function(object) { state.optim = rownames(subset(object$start, type == "state")) deparms.optim = rownames(subset(object$start, type == "deparm")) fractions.optim = grep("^f_", deparms.optim, value = TRUE) + N.optim = grep("^N_", deparms.optim, value = TRUE) if ("g" %in% deparms.optim) fractions.optim <- c("g", fractions.optim) - rates.optim.unsorted = setdiff(deparms.optim, fractions.optim) + rates.optim.unsorted = setdiff(deparms.optim, union(fractions.optim, N.optim)) rates.optim <- rownames(object$start[rates.optim.unsorted, ]) n.plot <- c(state.optim = length(state.optim), rates.optim = length(rates.optim), + N.optim = length(N.optim), fractions.optim = length(fractions.optim)) n.plot <- n.plot[n.plot > 0] @@ -42,6 +44,8 @@ mkinparplot <- function(object) { na.rm = TRUE, finite = TRUE), rates.optim = range(c(0, unlist(values)), na.rm = TRUE, finite = TRUE), + N.optim = range(c(0, 1, unlist(values)), + na.rm = TRUE, finite = TRUE), fractions.optim = range(c(0, 1, unlist(values)), na.rm = TRUE, finite = TRUE)) stripchart(values["Estimate", ][length(parnames):1], @@ -49,7 +53,7 @@ mkinparplot <- function(object) { ylim = c(0.5, length(get(type)) + 0.5), yaxt = "n") if (type %in% c("rates.optim", "fractions.optim")) abline(v = 0, lty = 2) - if (type %in% c("fractions.optim")) abline(v = 1, lty = 2) + if (type %in% c("N.optim", "fractions.optim")) abline(v = 1, lty = 2) position <- ifelse(values["Estimate", ] < mean(xlim), "right", "left") text(ifelse(position == "left", min(xlim), max(xlim)), length(parnames):1, parnames, diff --git a/R/mkinpredict.R b/R/mkinpredict.R index 14efc568..da013d50 100644 --- a/R/mkinpredict.R +++ b/R/mkinpredict.R @@ -50,6 +50,12 @@ mkinpredict <- function(mkinmod, odeparms, odeini, FOMC = FOMC.solution(outtimes, evalparse(parent.name), evalparse("alpha"), evalparse("beta")), + IORE = IORE.solution(outtimes, + evalparse(parent.name), + ifelse(mkinmod$use_of_ff == "min", + evalparse(paste("k.iore", parent.name, "sink", sep="_")), + evalparse(paste("k.iore", parent.name, sep="_"))), + evalparse("N_parent")), DFOP = DFOP.solution(outtimes, evalparse(parent.name), evalparse("k1"), evalparse("k2"), @@ -103,6 +109,10 @@ mkinpredict <- function(mkinmod, odeparms, odeini, rtol = rtol, ... ) + if (sum(is.na(out)) > 0) { + stop("Differential equations were not integrated for all output times because\n", + "NaN values occurred in output from ode()") + } } if (map_output) { # Output transformation for models with unobserved compartments like SFORB diff --git a/R/mkinsub.R b/R/mkinsub.R new file mode 100644 index 00000000..f92af54b --- /dev/null +++ b/R/mkinsub.R @@ -0,0 +1,23 @@ +# Copyright (C) 2014 Johannes Ranke +# Portions of this code are copyright (C) 2013 Eurofins Regulatory AG +# Contact: jranke@uni-bremen.de + +# This file is part of the R package mkin + +# mkin is free software: you can redistribute it and/or modify it under the +# terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. + +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. + +# You should have received a copy of the GNU General Public License along with +# this program. If not, see <http://www.gnu.org/licenses/> +mkinsub <- function(submodel, to = NULL, sink = TRUE) +{ + return(list(type = submodel, to = to, sink = sink)) +} +# vim: set ts=2 sw=2 expandtab: diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R index e64bac59..aa70a0b3 100644 --- a/R/transform_odeparms.R +++ b/R/transform_odeparms.R @@ -33,13 +33,18 @@ transform_odeparms <- function(parms, mkinmod, # Log transformation for rate constants if requested
k <- parms[grep("^k_", names(parms))]
+ k.iore <- parms[grep("^k.iore_", names(parms))]
+ k <- c(k, k.iore)
if (length(k) > 0) {
if(transform_rates) {
transparms[paste0("log_", names(k))] <- log(k)
- }
- else transparms[names(k)] <- k
+ } else transparms[names(k)] <- k
}
+ # Do not transform exponents in IORE models
+ N <- parms[grep("^N", names(parms))]
+ transparms[names(N)] <- N
+
# Go through state variables and apply isometric logratio transformation to
# formation fractions if requested
mod_vars = names(spec)
@@ -108,15 +113,23 @@ backtransform_odeparms <- function(transparms, mkinmod, # Exponential transformation for rate constants
if(transform_rates) {
trans_k <- transparms[grep("^log_k_", names(transparms))]
+ trans_k.iore <- transparms[grep("^log_k.iore_", names(transparms))]
+ trans_k = c(trans_k, trans_k.iore)
if (length(trans_k) > 0) {
- k_names <- gsub("^log_k_", "k_", names(trans_k))
+ k_names <- gsub("^log_k", "k", names(trans_k))
parms[k_names] <- exp(trans_k)
}
} else {
trans_k <- transparms[grep("^k_", names(transparms))]
parms[names(trans_k)] <- trans_k
+ trans_k.iore <- transparms[grep("^k.iore_", names(transparms))]
+ parms[names(trans_k.iore)] <- trans_k.iore
}
+ # Do not transform exponents in IORE models
+ N <- transparms[grep("^N", names(transparms))]
+ parms[names(N)] <- N
+
# Go through state variables and apply inverse isometric logratio transformation
mod_vars = names(spec)
for (box in mod_vars) {
@@ -118,6 +118,9 @@ documentation or the package vignettes referenced from the `reweight = "obs"` to your call to `mkinfit` and a separate variance componenent for each of the observed variables will be optimised in a second stage after the primary optimisation algorithm has converged. +* When a metabolite decline phase is not described well by SFO kinetics, + either IORE kinetics or SFORB kinetics can be used for the metabolite, + adding one respectively two parameters to the system. ## GUI diff --git a/man/FOMC.solution.Rd b/man/FOMC.solution.Rd index d04d34e1..f4a26e41 100644 --- a/man/FOMC.solution.Rd +++ b/man/FOMC.solution.Rd @@ -44,6 +44,6 @@ FOMC.solution(t, parent.0, alpha, beta) Technology} \bold{24}, 1032-1038
}
\examples{
- \dontrun{plot(function(x) FOMC.solution(x, 100, 10, 2), 0, 2)}
+ \dontrun{plot(function(x) FOMC.solution(x, 100, 10, 2), 0, 2, ylim = c(0, 100))}
}
\keyword{ manip }
diff --git a/man/IORE.solution.Rd b/man/IORE.solution.Rd new file mode 100644 index 00000000..65dac9a7 --- /dev/null +++ b/man/IORE.solution.Rd @@ -0,0 +1,34 @@ +\name{IORE.solution}
+\Rdversion{1.1}
+\alias{IORE.solution}
+\title{ Indeterminate order rate equation kinetics }
+\description{
+ Function describing exponential decline from a defined starting value, with
+ a concentration dependent rate constant.
+}
+\usage{
+ IORE.solution(t, parent.0, k.iore, N)
+}
+\arguments{
+ \item{t}{ Time. }
+ \item{parent.0}{ Starting value for the response variable at time zero. }
+ \item{k.iore}{ Rate constant. Note that this depends on the concentration units used. }
+ \item{N}{ Exponent describing the nonlinearity of the rate equation }
+}
+\note{
+ The solution of the IORE kinetic model reduces to the
+ \code{\link{SFO.solution}} if N = 1.
+}
+\value{
+ The value of the response variable at time \code{t}.
+}
+\references{
+ NAFTA Technical Working Group on Pesticides (not dated) Guidance for
+ Evaluating and Calculating Degradation Kinetics in Environmental
+ Media
+}
+\examples{
+ \dontrun{plot(function(x) IORE.solution(x, 100, 0.2, 1.3), 0, 2,
+ ylim = c(0, 100))}
+}
+\keyword{ manip }
diff --git a/man/mkinfit.Rd b/man/mkinfit.Rd index f4174c29..9c82f5ff 100644 --- a/man/mkinfit.Rd +++ b/man/mkinfit.Rd @@ -216,6 +216,12 @@ mkinfit(mkinmod, observed, Gao). A similar implemention can also be found in CAKE 2.0, which is the other GUI derivative of mkin, sponsored by Syngenta. } +\note{ + When using the "IORE" submodel for metabolites, fitting with + "transform_rates = TRUE" (the default) often leads to failures of the + numerical ODE solver. In this situation it may help to switch off the + internal rate transformation. +} \author{ Johannes Ranke } diff --git a/man/mkinmod.Rd b/man/mkinmod.Rd index 76127c58..e6749ded 100644 --- a/man/mkinmod.Rd +++ b/man/mkinmod.Rd @@ -4,9 +4,13 @@ Function to set up a kinetic model with one or more state variables.
}
\description{
- The function takes a specification, consisting of a list of the observed variables
- in the data. Each observed variable is again represented by a list, specifying the
- kinetic model type and reaction or transfer to other observed compartments.
+ The function usually takes several expressions, each assigning a compound name to
+ a list, specifying the kinetic model type and reaction or transfer to other
+ observed compartments. Instead of specifying several expressions, a list
+ of lists can be given in the speclist argument.
+
+ For the definition of model types and their parameters, the equations given
+ in the FOCUS and NAFTA guidance documents are used.
}
\usage{
mkinmod(..., use_of_ff = "min", speclist = NULL)
@@ -15,7 +19,8 @@ mkinmod(..., use_of_ff = "min", speclist = NULL) \item{...}{
For each observed variable, a list has to be specified as an argument, containing
at least a component \code{type}, specifying the type of kinetics to use
- for the variable. Currently, single first order kinetics "SFO" or
+ for the variable. Currently, single first order kinetics "SFO",
+ indeterminate order rate equation kinetics "IORE", or
single first order with reversible binding "SFORB" are implemented for all
variables, while
"FOMC", "DFOP" and "HS" can additionally be chosen for the first
@@ -46,7 +51,19 @@ mkinmod(..., use_of_ff = "min", speclist = NULL) \item{map}{ A list containing named character vectors for each observed variable, specifying
the modelling variables by which it is represented. }
\item{use_of_ff}{ The content of \code{use_of_ff} is passed on in this list component. }
- \item{coefmat}{ The coefficient matrix, if the system of differential equations can be represented by one. }
+ \item{coefmat}{ The coefficient matrix, if the system of differential equations can be
+ represented by one. }
+}
+\references{
+ FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence and
+ Degradation Kinetics from Environmental Fate Studies on Pesticides in EU
+ Registration} Report of the FOCUS Work Group on Degradation Kinetics,
+ EC Document Reference Sanco/10058/2005 version 2.0, 434 pp,
+ \url{http://focus.jrc.ec.europa.eu/dk}
+
+ NAFTA Technical Working Group on Pesticides (not dated) Guidance for
+ Evaluating and Calculating Degradation Kinetics in Environmental
+ Media
}
\author{
Johannes Ranke
diff --git a/man/mkinpredict.Rd b/man/mkinpredict.Rd index 97db90e3..7d8979e4 100644 --- a/man/mkinpredict.Rd +++ b/man/mkinpredict.Rd @@ -65,21 +65,25 @@ } \examples{ SFO <- mkinmod(degradinol = list(type = "SFO")) + # Compare solution types mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, solution_type = "analytical") mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, + solution_type = "deSolve") + mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, solution_type = "eigen") - mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 1:20, - solution_type = "analytical")[20,] - mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, - atol = 1e-20)[20,] - # The integration method does not make a lot of difference + # Compare integration methods to analytical solution mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, - atol = 1e-20, method = "ode45")[20,] + solution_type = "analytical")[21,] mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, - atol = 1e-20, method = "rk4")[20,] + method = "lsoda")[21,] + mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, + method = "ode45")[21,] + mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, + method = "rk4")[21,] + # rk4 is not as precise here # The number of output times used to make a lot of difference until the # default for atol was adjusted @@ -87,5 +91,6 @@ seq(0, 20, by = 0.1))[201,] mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), seq(0, 20, by = 0.01))[2001,] + } \keyword{ manip } diff --git a/man/mkinsub.Rd b/man/mkinsub.Rd new file mode 100644 index 00000000..637a671c --- /dev/null +++ b/man/mkinsub.Rd @@ -0,0 +1,44 @@ +\name{mkinsub} +\alias{mkinsub} +\title{ + Function to set up a kinetic submodel for one state variable +} +\description{ + This is a convenience function to set up the lists used as arguments for + \code{\link{mkinmod}}. +} +\usage{ +mkinsub(type, to = NULL, sink = TRUE) +} +\arguments{ + \item{type}{ + Character vector of length one to specify the submodel type. See + \code{\link{mkinmod}} for the list of allowed submodel names. + } + \item{to}{ + Vector of the names of the state variable to which a transformation + shall be included in the model. + } + \item{sink}{ + Should a pathway to sink be included in the model in addition to the + pathways to other state variables? + } +} +\value{ + A list for use with \code{\link{mkinmod}}. +} +\author{ + Johannes Ranke +} +\examples{ +# One parent compound, one metabolite, both single first order. +SFO_SFO <- mkinmod( + parent = list(type = "SFO", to = "m1"), + m1 = list(type = "SFO")) + +# The same model using mkinsub +SFO_SFO.2 <- mkinmod( + parent = mkinsub("SFO", "m1"), + m1 = mkinsub("SFO")) +} + diff --git a/vignettes/FOCUS_Z.pdf b/vignettes/FOCUS_Z.pdf Binary files differindex 0013cd5e..a22a3a2e 100644 --- a/vignettes/FOCUS_Z.pdf +++ b/vignettes/FOCUS_Z.pdf diff --git a/vignettes/mkin.pdf b/vignettes/mkin.pdf Binary files differindex 83182e65..d14b86bf 100644 --- a/vignettes/mkin.pdf +++ b/vignettes/mkin.pdf |