Overhaul of mkinmod in order to make it possible to choose between different model specifications -

a variant with minimal use of formation fractions, and a variant with formation fractions everywhere.



git-svn-id: svn+ssh://svn.r-forge.r-project.org/svnroot/kinfit/pkg/mkin@26 edb9625f-4e0d-4859-8d74-9fd3b1da38cb

4 files changed, 122 insertions, 55 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index a7f33e2..509c653 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -2,8 +2,8 @@ Package: mkin
 Type: Package
 Title: Routines for fitting kinetic models with one or more state
   variables to chemical degradation data
-Version: 0.9-01
-Date: 2012-04-10
+Version: 0.9-02
+Date: 2012-04-17
 Author: Johannes Ranke, Katrin Lindenberger, René Lehmann
 Maintainer: Johannes Ranke <jranke@uni-bremen.de>
 Description: Calculation routines based on the FOCUS Kinetics Report (2006).
diff --git a/R/mkinmod.R b/R/mkinmod.R
index 8095932..54b3b5f 100644
--- a/R/mkinmod.R
+++ b/R/mkinmod.R
@@ -18,10 +18,12 @@
 # You should have received a copy of the GNU General Public License along with

 # this program. If not, see <http://www.gnu.org/licenses/>#}}}

 

-mkinmod <- function(...)

+mkinmod <- function(..., use_of_ff = "min")

 {

   spec <- list(...)

   obs_vars <- names(spec)

+  if (!use_of_ff %in% c("min", "max"))

+    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

@@ -51,7 +53,7 @@ mkinmod <- function(...)
   {

     # Check the type component of the compartment specification#{{{

     if(is.null(spec[[varname]]$type)) stop(

-      "Every argument to mkinmod must be a list containing a type component")

+      "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("FOMC", "DFOP", "HS") & match(varname, obs_vars) != 1) {

@@ -77,46 +79,66 @@ mkinmod <- function(...)
   # Create content of differential equations and build parameter list#{{{

   for (varname in obs_vars)

   {

-    # Add first-order term to first (or only) box for SFO and SFORB#{{{

+    # Get the name of the box(es) we are working on for the decline term(s)

     box_1 = map[[varname]][[1]] # This is the only box unless type is SFORB

-    if(spec[[varname]]$type %in% c("SFO", "SFORB")) {

-      k_compound <- paste("k", box_1, sep="_")

-      origin_term <- paste(k_compound, "*", box_1)

-      parms <- c(parms, k_compound)

+    if(spec[[varname]]$type %in% c("SFO", "SFORB")) { # {{{ Add SFO or SFORB decline term

+      if (use_of_ff == "min") { # Minimum use of formation fractions

+        # 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]]$sink) {

+          # If sink is required, add first-order sink term

+          k_compound_sink <- paste("k", box_1, "sink", sep="_")

+          parms <- c(parms, k_compound_sink)

+          decline_term <- paste(k_compound_sink, "*", box_1)

+        } else { # otherwise no decline term needed here

+          decline_term = "" 

+        }

+      } else {

+        k_compound <- paste("k", box_1, sep="_")

+        parms <- c(parms, k_compound)

+        decline_term <- paste(k_compound, "*", box_1)

+      }

     }#}}}

-    # Construct and add FOMC term and add FOMC parameters if needed#{{{

-    if(spec[[varname]]$type == "FOMC") {

+    if(spec[[varname]]$type == "FOMC") { # {{{ Add FOMC decline term

       # From p. 53 of the FOCUS kinetics report

-      origin_term <- paste("(alpha/beta) * ((time/beta) + 1)^-1 *", box_1)

+      decline_term <- paste("(alpha/beta) * ((time/beta) + 1)^-1 *", box_1)

       parms <- c(parms, "alpha", "beta")

     } #}}}

-    # Construct and add DFOP term and add DFOP parameters if needed#{{{

-    if(spec[[varname]]$type == "DFOP") {

+    if(spec[[varname]]$type == "DFOP") { # {{{ Add DFOP decline term

       # From p. 57 of the FOCUS kinetics report

-      origin_term <- paste("((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) *", box_1)

+      decline_term <- paste("((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) *", box_1)

       parms <- c(parms, "k1", "k2", "g")

     } #}}}

-    # Construct and add HS term and add HS parameters if needed#{{{

-    if(spec[[varname]]$type == "HS") {

+    if(spec[[varname]]$type == "HS") { # {{{ Add HS decline term

       # From p. 55 of the FOCUS kinetics report

-      origin_term <- paste("ifelse(time <= tb, k1, k2)", "*", box_1)

+      decline_term <- paste("ifelse(time <= tb, k1, k2)", "*", box_1)

       parms <- c(parms, "k1", "k2", "tb")

     } #}}}

     # Add origin decline term to box 1 (usually the only box, unless type is SFORB)#{{{

-    diffs[[box_1]] <- paste(diffs[[box_1]], "-", origin_term)#}}}

-    # Add reversible binding for SFORB models#{{{

-    if(spec[[varname]]$type == "SFORB") {

+    diffs[[box_1]] <- paste(diffs[[box_1]], "-", decline_term)#}}}

+    if(spec[[varname]]$type == "SFORB") { # {{{ Add SFORB reversible binding terms

       box_2 = map[[varname]][[2]]

-      k_free_bound <- paste("k", varname, "free", "bound", sep="_")

-      k_bound_free <- paste("k", varname, "bound", "free", sep="_")

-      reversible_binding_term_1 <- paste("-", k_free_bound, "*", box_1, "+",

-        k_bound_free, "*", box_2)

+      if (use_of_ff == "min") { # Minimum use of formation fractions

+        k_free_bound <- paste("k", varname, "free", "bound", sep="_")

+        k_bound_free <- paste("k", varname, "bound", "free", sep="_")

+        parms <- c(parms, k_free_bound, k_bound_free)

+        reversible_binding_term_1 <- paste("-", k_free_bound, "*", box_1, "+",

+          k_bound_free, "*", box_2)

+        reversible_binding_term_2 <- paste("+", k_free_bound, "*", box_1, "-",

+          k_bound_free, "*", box_2)

+      } else { # Use formation fractions also for the free compartment

+        f_free_bound <- paste("f", varname, "free", "bound", sep="_")

+        k_bound_free <- paste("k", varname, "bound", "free", sep="_")

+        parms <- c(parms, f_free_bound, k_bound_free)

+        reversible_binding_term_1 <- paste("+", k_bound_free, "*", box_2)

+        reversible_binding_term_2 <- paste("+", f_free_bound, "*", k_compound, "*", box_1, "-",

+          k_bound_free, "*", box_2)

+      }

       diffs[[box_1]] <- paste(diffs[[box_1]], reversible_binding_term_1)

-      reversible_binding_term_2 <- paste("+", k_free_bound, "*", box_1, "-",

-        k_bound_free, "*", box_2)

       diffs[[box_2]] <- paste(diffs[[box_2]], reversible_binding_term_2)

-      parms <- c(parms, k_free_bound, k_bound_free)

     } #}}}

+

     # Transfer between compartments#{{{

     to <- spec[[varname]]$to

     if(!is.null(to)) {

@@ -128,45 +150,79 @@ mkinmod <- function(...)
         target_box <- switch(spec[[target]]$type,

           SFO = target,

           SFORB = paste(target, "free", sep="_"))

-        fraction_to_target = paste("f", origin_box, "to", target, sep="_")

-        diffs[[target_box]] <- paste(diffs[[target_box]], "+", 

-            fraction_to_target, "*", origin_term)

-        parms <- c(parms, fraction_to_target)

+        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]], "-", 

+            k_from_to, "*", origin_box)

+          diffs[[target_box]] <- paste(diffs[[target_box]], "+", 

+            k_from_to, "*", origin_box)

+        } else {

+          fraction_to_target = paste("f", origin_box, "to", target, sep="_")

+          parms <- c(parms, fraction_to_target)

+          diffs[[target_box]] <- paste(diffs[[target_box]], "+", 

+              fraction_to_target, "*", decline_term)

+        }

       }

-    }#}}}

+    } #}}}

   }#}}}

 

-  model <- list(diffs = diffs, parms = parms, map = map)

+  model <- list(diffs = diffs, parms = parms, map = map, use_of_ff = use_of_ff)

 

   # Create coefficient matrix if appropriate#{{{

   if (mat) {

     boxes <- names(diffs)

     n <- length(boxes)

     m <- matrix(nrow=n, ncol=n, dimnames=list(boxes, boxes))

-    for (from in boxes) {

-      for (to in boxes) {

-        if (from == to) { # diagonal elements

-          k.candidate = paste("k", from, sep="_")

-          m[from,to] = ifelse(k.candidate %in% model$parms,

-              paste("-", k.candidate), "0")

-          if(grepl("_free", from)) { # add transfer to bound compartment for SFORB

-            m[from,to] = paste(m[from,to], "-", paste("k", from, "bound", sep="_"))

+

+    if (use_of_ff == "min") { # {{{ Minimum use of formation fractions

+      for (from in boxes) {

+        for (to in boxes) {

+          if (from == to) { # diagonal elements

+            k.candidate = paste("k", from, c(boxes, "sink"), sep="_")

+            k.candidate = sub("free.*bound", "free_bound", k.candidate)

+            k.candidate = sub("bound.*free", "bound_free", k.candidate)

+            k.effective = intersect(model$parms, k.candidate)

+            m[from,to] = ifelse(length(k.effective) > 0,

+                paste("-", k.effective, collapse = " "), "0")

+

+          } else {          # off-diagonal elements

+            k.candidate = paste("k", from, to, sep="_")

+            k.candidate = sub("free.*bound", "free_bound", k.candidate)

+            k.candidate = sub("bound.*free", "bound_free", k.candidate)

+            k.effective = intersect(model$parms, k.candidate)

+            m[to, from] = ifelse(length(k.effective) > 0,

+                k.effective, "0")

           }

-          if(grepl("_bound", from)) { # add backtransfer to free compartment for SFORB

-            m[from,to] = paste("- k", from, "free", sep="_")

+        }

+      } # }}}

+    } else { # {{{ Use formation fractions where possible

+      for (from in boxes) {

+        for (to in boxes) {

+          if (from == to) { # diagonal elements

+            k.candidate = paste("k", from, sep="_")

+            m[from,to] = ifelse(k.candidate %in% model$parms,

+                paste("-", k.candidate), "0")

+            if(grepl("_free", from)) { # add transfer to bound compartment for SFORB

+              m[from,to] = paste(m[from,to], "-", paste("k", from, "bound", sep="_"))

+            }

+            if(grepl("_bound", from)) { # add backtransfer to free compartment for SFORB

+              m[from,to] = paste("- k", from, "free", sep="_")

+            }

+            m[from,to] = m[from,to]

+          } else {          # off-diagonal elements

+            f.candidate = paste("f", from, "to", to, sep="_")

+            k.candidate = paste("k", from, to, sep="_")

+            k.candidate = sub("free.*bound", "free_bound", k.candidate)

+            k.candidate = sub("bound.*free", "bound_free", k.candidate)

+            m[to, from] = ifelse(f.candidate %in% model$parms,

+              paste(f.candidate, " * k_", from, sep=""), 

+              ifelse(k.candidate %in% model$parms, k.candidate, "0"))

           }

-          m[from,to] = m[from,to]

-        } else {          # off-diagonal elements

-          f.candidate = paste("f", from, "to", to, sep="_")

-          k.candidate = paste("k", from, to, sep="_")

-          k.candidate = sub("free.*bound", "free_bound", k.candidate)

-          k.candidate = sub("bound.*free", "bound_free", k.candidate)

-          m[to, from] = ifelse(f.candidate %in% model$parms,

-            paste(f.candidate, " * k_", from, sep=""), 

-            ifelse(k.candidate %in% model$parms, k.candidate, "0"))

         }

       }

-    }

+    } # }}}

     model$coefmat <- m

   }#}}}

 

diff --git a/TODO b/TODO
index da6c5f5..d6e567e 100644
--- a/TODO
+++ b/TODO
@@ -1,3 +1,4 @@
+- Fix coefmat in mkinmod after transition to both possibilities of model formation regarding ff
 - Adapt mkinplot function
 - Calculate confidence intervals for parameters assuming normal distribution
 - Calculate confidence intervals for DT50 and DT90 values when only one parameter is involved
diff --git a/man/mkinmod.Rd b/man/mkinmod.Rd
index 16ad38a..b528d83 100644
--- a/man/mkinmod.Rd
+++ b/man/mkinmod.Rd
@@ -9,7 +9,7 @@
   kinetic model type and reaction or transfer to other observed compartments.

 }

 \usage{

-mkinmod(...)

+mkinmod(..., use_of_ff = "min")

 }

 \arguments{

   \item{...}{

@@ -23,6 +23,15 @@ mkinmod(...)
     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

     model.

+    If the argument \code{use_of_ff} is set to "min" (default) and the model for 

+    the compartment is "SFO" or "SFORB", an additional component of the list

+    can be "sink=FALSE" effectively fixing the flux to sink to zero.

+  }

+  \item{use_of_ff}{

+    Specification of the use of formation fractions in the model equations and, if 

+    applicable, the coefficient matrix. If "min", a minimum use of formation 

+    fractions is made in order to avoid fitting the product of formation fractions 

+    and rate constants. If "max", formation fractions are always used.

   }

 }

 \value{

@@ -32,6 +41,7 @@ mkinmod(...)
     \item{parms}{ A vector of parameter names occurring in the differential equations. }

     \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. }

 }

 \author{