- Moved the call to mkinerrmin to summary.mkinfit

- The argument to mkinerrmin is now an object of class mkinfit
- Fixed the allocation of parameters to observed variables in mkinerrmin



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

2 files changed, 71 insertions, 46 deletions

diff --git a/R/mkinerrmin.R b/R/mkinerrmin.R
index 648ba45c..5415b41d 100644
--- a/R/mkinerrmin.R
+++ b/R/mkinerrmin.R
@@ -1,7 +1,7 @@
 # $Id$

 

-# Copyright (C) 2010 Johannes Ranke

-# Contact: mkin-devel@lists.berlios.de

+# Copyright (C) 2010-2012 Johannes Ranke

+# Contact: jranke@uni-bremen.de

 

 # This file is part of the R package mkin

 

@@ -18,17 +18,63 @@
 # You should have received a copy of the GNU General Public License along with

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

 

-mkinerrmin <- function(errdata, n.parms, alpha = 0.05)

+mkinerrmin <- function(fit, alpha = 0.05)

 {

-  means.mean <- mean(errdata$value_mean, na.rm=TRUE)

-

-	df = length(errdata$value_mean) - n.parms

+  parms.optim <- fit$parms.all

+  kinerrmin <- function(errdata, n.parms) {

+    means.mean <- mean(errdata$value_mean, na.rm=TRUE)

+    df = length(errdata$value_mean) - n.parms

   

-	f <- function(err)

-	{

-		(sum((errdata$value_mean - errdata$value_pred)^2/((err * means.mean)^2)) - 

-		 qchisq(1 - alpha,df))^2

-	}

-	err.min <- optimize(f, c(0.01,0.9))$minimum

-	return(list(err.min = err.min, n.optim = n.parms, df = df))

+    f <- function(err)

+    {

+      (sum((errdata$value_mean - errdata$value_pred)^2/((err * means.mean)^2)) - 

+       qchisq(1 - alpha,df))^2

+    }

+    err.min <- optimize(f, c(0.01,0.9))$minimum

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

+  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 are attributed to the source variable

+    n.k.optim <- length(grep(paste("^k", 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 <- 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]]) {

+      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 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 cb0396f7..b2641e64 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -175,9 +175,9 @@ mkinfit <- function(mkinmod, observed,
   fit$mkinmod <- mkinmod

 

   # We need data and predictions for summary and plotting

-  fit$observed <- mkin_long_to_wide(observed)

-  predicted_long <- mkin_wide_to_long(out_predicted, time = "time")

-  fit$predicted <- out_predicted

+  fit$observed <- observed

+  fit$obs_vars <- obs_vars

+  fit$predicted <- mkin_wide_to_long(out_predicted, time = "time")

 

   # Collect initial parameter values in two dataframes

   fit$start <- data.frame(initial = c(state.ini.optim, 

@@ -189,35 +189,11 @@ mkinfit <- function(mkinmod, observed,
     value = c(state.ini.fixed, parms.fixed))

   fit$fixed$type = c(rep("state", length(state.ini.fixed)), rep("deparm", length(parms.fixed)))

 

-  # Calculate chi2 error levels according to FOCUS (2006)

-  means <- aggregate(value ~ time + name, data = observed, mean, na.rm=TRUE)

-  errdata <- merge(means, predicted_long, by = c("time", "name"), suffixes = c("_mean", "_pred"))

-  errdata <- errdata[order(errdata$time, errdata$name), ]

-  errmin.overall <- mkinerrmin(errdata, length(parms.optim) + length(state.ini.optim))

-  

-  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 obs_vars)

-  {

-    errdata.var <- subset(errdata, name == obs_var)

-    n.k.optim <- length(grep(paste("k", obs_var, sep="_"), names(parms.optim)))

-    n.initials.optim <- length(grep(paste(obs_var, ".*", "_0", sep=""), names(state.ini.optim)))

-    n.optim <- n.k.optim + n.initials.optim

-    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

-    errmin.tmp <- mkinerrmin(errdata.var, n.optim)

-    errmin[obs_var, c("err.min", "n.optim", "df")] <- errmin.tmp

-  }

-  fit$errmin <- errmin

 

-  # Calculate dissipation times DT50 and DT90 and, if necessary, formation fractions

-  # from optimised parameters

   parms.all = backtransform_odeparms(c(fit$par, parms.fixed), mod_vars)

+

+  # Calculate dissipation times DT50 and DT90 and, if necessary, formation

+  # fractions and SFORB eigenvalues from optimised parameters

   fit$distimes <- data.frame(DT50 = rep(NA, length(obs_vars)), DT90 = rep(NA, length(obs_vars)), 

     row.names = obs_vars)

   fit$ff <- vector()

@@ -305,7 +281,7 @@ mkinfit <- function(mkinmod, observed,
   }

 

   # Collect observed, predicted and residuals

-  data <- merge(observed, predicted_long, by = c("time", "name"))

+  data <- merge(fit$observed, fit$predicted, by = c("time", "name"))

   names(data) <- c("time", "variable", "observed", "predicted")

   data$residual <- data$observed - data$predicted

   data$variable <- ordered(data$variable, levels = obs_vars)

@@ -362,9 +338,12 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, ...) {
   ans$start <- object$start

 

   ans$fixed <- object$fixed

-  ans$errmin <- object$errmin 

+

+  ans$errmin <- mkinerrmin(object, alpha = 0.05)

+

   ans$parms.all <- object$parms.all

-  ans$ff <- object$ff

+  if (!is.null(object$ff))

+    ans$ff <- object$ff

   if(distimes) ans$distimes <- object$distimes

   if(length(object$SFORB) != 0) ans$SFORB <- object$SFORB

   class(ans) <- c("summary.mkinfit", "summary.modFit")