Merge bugfix branch 'master' into iore

Add IORE support to mkinerrmin

1 files changed, 100 insertions, 85 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)
+}