aboutsummaryrefslogtreecommitdiff
path: root/R
diff options
context:
space:
mode:
Diffstat (limited to 'R')
-rw-r--r--R/IORE.solution.R4
-rw-r--r--R/endpoints.R28
-rw-r--r--R/mkinerrmin.R198
-rw-r--r--R/mkinfit.R4
-rw-r--r--R/mkinmod.R54
-rw-r--r--R/mkinparplot.R8
-rw-r--r--R/mkinpredict.R10
-rw-r--r--R/transform_odeparms.R14
8 files changed, 201 insertions, 119 deletions
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..a6889b0b 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -116,13 +116,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
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/transform_odeparms.R b/R/transform_odeparms.R
index e64bac59..778f56cd 100644
--- a/R/transform_odeparms.R
+++ b/R/transform_odeparms.R
@@ -33,6 +33,8 @@ 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)
@@ -40,6 +42,10 @@ transform_odeparms <- function(parms, mkinmod,
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,8 +114,10 @@ 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 {
@@ -117,6 +125,10 @@ backtransform_odeparms <- function(transparms, mkinmod,
parms[names(trans_k)] <- trans_k
}
+ # 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) {

Contact - Imprint