Nearly complete support for IORE, pending mkinerrmin

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 }