From a69bf39427ff4f93eebdc8bceacb8174ff13c085 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Mon, 14 Jul 2014 19:07:54 +0200 Subject: Nearly complete support for IORE, pending mkinerrmin --- NEWS.md | 6 ++++++ R/endpoints.R | 28 +++++++++++++++++++++++++--- R/mkinerrmin.R | 8 +++++--- R/mkinmod.R | 26 ++++++++++++++++---------- R/mkinparplot.R | 8 ++++++-- R/mkinpredict.R | 4 ++++ README.md | 3 +++ man/mkinmod.Rd | 26 +++++++++++++++++++++----- man/mkinpredict.Rd | 19 ++++++++++++------- 9 files changed, 98 insertions(+), 30 deletions(-) diff --git a/NEWS.md b/NEWS.md index 6c23f316..2281bcc2 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,9 @@ +# CHANGES in mkin VERSION 0.9-31 + +## NEW FEATURES + +- 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'). The numerical solution (method.ode = 'deSolve') of the IORE differential equations sometimes fails due to numerical problems. + # CHANGES in mkin VERSION 0.9-30 ## NEW FEATURES 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 671bcaab..d89a2f91 100644 --- a/R/mkinerrmin.R +++ b/R/mkinerrmin.R @@ -55,13 +55,15 @@ mkinerrmin <- function(fit, alpha = 0.05) # 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 + # 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.iore.optim <- length(grep(paste("^k.iore", obs_var, sep="_"), names(parms.optim))) + n.N.optim <- length(grep(paste("^N", obs_var, sep="_"), names(parms.optim))) # Formation fractions are attributed to the target variable n.ff.optim <- length(grep(paste("^f", ".*", obs_var, "$", sep=""), names(parms.optim))) - n.optim <- n.k.optim + n.initials.optim + n.ff.optim + 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 @@ -74,7 +76,7 @@ mkinerrmin <- function(fit, alpha = 0.05) if ("tb" %in% names(parms.optim)) n.optim <- n.optim + 1 } - # Calculate and add a line to the results + # 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 } diff --git a/R/mkinmod.R b/R/mkinmod.R index 321887fc..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,12 +33,10 @@ 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, IORE, DFOP or HS is used for the parent {{{ @@ -49,6 +47,8 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL) # 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 {{{ @@ -59,7 +59,8 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL) 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, @@ -88,7 +89,11 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL) 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 = "_") @@ -173,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 43b5d2f3..da013d50 100644 --- a/R/mkinpredict.R +++ b/R/mkinpredict.R @@ -109,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/README.md b/README.md index 4d394808..b968a72a 100644 --- a/README.md +++ b/README.md @@ -113,6 +113,9 @@ 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/mkinmod.Rd b/man/mkinmod.Rd index 63087d49..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,10 +19,11 @@ 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", "IORE", "DFOP" and "HS" can additionally be chosen for the first + "FOMC", "DFOP" and "HS" can additionally be chosen for the first variable which is assumed to be the source compartment. Additionally, each component of the list can include a character vector \code{to}, specifying names of variables to which a transfer is to be assumed in the @@ -49,6 +54,17 @@ mkinmod(..., use_of_ff = "min", speclist = NULL) \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 } -- cgit v1.2.1