aboutsummaryrefslogtreecommitdiff
path: root/R
diff options
context:
space:
mode:
authorJohannes Ranke <jranke@uni-bremen.de>2014-07-14 20:18:53 +0200
committerJohannes Ranke <jranke@uni-bremen.de>2014-07-14 20:18:53 +0200
commite0a3413892c7330d496b448a561e87d2bdd67aa9 (patch)
tree20227b52969496a03daf1fe01588783e4a112501 /R
parenta69bf39427ff4f93eebdc8bceacb8174ff13c085 (diff)
parent759e693e9af8e794bbfa62b001117fabbdbc8bfa (diff)
Merge bugfix branch 'master' into iore
Add IORE support to mkinerrmin
Diffstat (limited to 'R')
-rw-r--r--R/mkinerrmin.R185
-rw-r--r--R/mkinfit.R36
2 files changed, 117 insertions, 104 deletions
diff --git a/R/mkinerrmin.R b/R/mkinerrmin.R
index d89a2f91..4137d33a 100644
--- a/R/mkinerrmin.R
+++ b/R/mkinerrmin.R
@@ -1,85 +1,100 @@
-# 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), ]
-
- # Any value that is set to exactly zero is not really an observed value
- # Remove those at time 0 - those are caused by the FOCUS recommendation
- # to set metabolites occurring at time 0 to 0
- errdata <- subset(errdata, !(time == 0 & value_mean == 0))
-
- 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"
-
- 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.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 <- 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]]) {
- if ("alpha" %in% names(parms.optim)) n.optim <- n.optim + 1
- if ("beta" %in% names(parms.optim)) n.optim <- n.optim + 1
- if ("k1" %in% names(parms.optim)) n.optim <- n.optim + 1
- if ("k2" %in% names(parms.optim)) n.optim <- n.optim + 1
- if ("g" %in% names(parms.optim)) n.optim <- n.optim + 1
- if ("tb" %in% names(parms.optim)) n.optim <- n.optim + 1
- }
-
- # 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)
-}
+# 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), ]
+
+ # Any value that is set to exactly zero is not really an observed value
+ # Remove those at time 0 - those are caused by the FOCUS recommendation
+ # to set metabolites occurring at time 0 to 0
+ errdata <- subset(errdata, !(time == 0 & value_mean == 0))
+
+ 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) {
+ for (target_var in fit$mkinmod$spec[[source_var]]$to) {
+ if (obs_var == target_var) {
+ n.ff.optim <- n.ff.optim +
+ length(grep(paste("^f", source_var, sep = "_"),
+ names(parms.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
+ 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 26fc0e6b..46121c6d 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -107,26 +107,24 @@ mkinfit <- function(mkinmod, observed,
if (parmname == "tb") parms.ini[parmname] = 5
if (parmname == "g") parms.ini[parmname] = 0.5
}
- # Default values for formation fractions in case they are used
- if (mkinmod$use_of_ff == "max") {
- for (box in mod_vars) {
- f_names <- mkinmod$parms[grep(paste0("^f_", box), mkinmod$parms)]
- if (length(f_names) > 0) {
- # We need to differentiate between default and specified fractions
- # and set the unspecified to 1 - sum(specified)/n_unspecified
- f_default_names <- intersect(f_names, defaultpar.names)
- f_specified_names <- setdiff(f_names, defaultpar.names)
- sum_f_specified = sum(parms.ini[f_specified_names])
- if (sum_f_specified > 1) {
- stop("Starting values for the formation fractions originating from ",
- box, " sum up to more than 1.")
- }
- if (mkinmod$spec[[box]]$sink) n_unspecified = length(f_default_names) + 1
- else {
- n_unspecified = length(f_default_names)
- }
- parms.ini[f_default_names] <- (1 - sum_f_specified) / n_unspecified
+ # Default values for formation fractions in case they are present
+ for (box in mod_vars) {
+ f_names <- mkinmod$parms[grep(paste0("^f_", box), mkinmod$parms)]
+ if (length(f_names) > 0) {
+ # We need to differentiate between default and specified fractions
+ # and set the unspecified to 1 - sum(specified)/n_unspecified
+ f_default_names <- intersect(f_names, defaultpar.names)
+ f_specified_names <- setdiff(f_names, defaultpar.names)
+ sum_f_specified = sum(parms.ini[f_specified_names])
+ if (sum_f_specified > 1) {
+ stop("Starting values for the formation fractions originating from ",
+ box, " sum up to more than 1.")
+ }
+ if (mkinmod$spec[[box]]$sink) n_unspecified = length(f_default_names) + 1
+ else {
+ n_unspecified = length(f_default_names)
}
+ parms.ini[f_default_names] <- (1 - sum_f_specified) / n_unspecified
}
}

Contact - Imprint