diff options
| -rw-r--r-- | DESCRIPTION | 4 | ||||
| -rw-r--r-- | R/mkinmod.R | 160 | ||||
| -rw-r--r-- | TODO | 1 | ||||
| -rw-r--r-- | man/mkinmod.Rd | 12 |
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
}#}}}
@@ -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{
|