- Added the reference fit data for FOCUS 2006 datasets from the kinfit package

- Used these data in unit tests for parent only models
- Fixed SFORB data and calculation of formation fractions along the way
- Reintroduced the test for the Schaefer 2007 data
- Got rid of the mkinmod unit tests - they are too hard to maintain and the mkinfit
  tests test the model definitions as well

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

2 files changed, 49 insertions, 6 deletions

diff --git a/R/mkinfit.R b/R/mkinfit.R
index cf58a7d..6e455e1 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -30,6 +30,7 @@ mkinfit <- function(mkinmod, observed,
   plot = FALSE, quiet = FALSE,

   err = NULL, weight = "none", scaleVar = FALSE,

   atol = 1e-6, n.outtimes = 100,

+  trace_parms = FALSE,

   ...)

 {

   # Get the names of the state variables in the model

@@ -112,7 +113,7 @@ mkinfit <- function(mkinmod, observed,
     assign("calls", calls+1, inherits=TRUE) # Increase the model solution counter

 

     # Trace parameter values if quiet is off

-    if(!quiet) cat(P, "\n")

+    if(trace_parms) cat(P, "\n")

 

     # Time points at which observed data are available

     # Make sure we include time 0, so initial values for state variables are for time 0

@@ -216,10 +217,12 @@ mkinfit <- function(mkinmod, observed,
   }

   fit$errmin <- errmin

 

-  # Calculate dissipation times DT50 and DT90 from parameters

+  # 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)

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

     row.names = obs_vars)

+  fit$ff <- vector()

   fit$SFORB <- vector()

   for (obs_var in obs_vars) {

     type = names(mkinmod$map[[obs_var]])[1]  

@@ -230,7 +233,7 @@ mkinfit <- function(mkinmod, observed,
       DT90 = log(10)/k_tot

       for (k_name in k_names)

       {

-        fit$ff[[sub("^k_", "", k_name)]] = parms.all[[k_name]] / k_tot

+        fit$ff[[sub("k_", "", k_name)]] = parms.all[[k_name]] / k_tot

       }

     }

     if (type == "FOMC") {

@@ -290,10 +293,12 @@ mkinfit <- function(mkinmod, observed,
       f_90 <- function(t) (SFORB_fraction(t) - 0.1)^2

       DT90.o <- optimize(f_90, c(0.01, max_DT))$minimum

       if (abs(DT90.o - max_DT) < 0.01) DT90 = NA else DT90 = DT90.o

+

       for (k_out_name in k_out_names)

       {

-        fit$ff[[sub("^k_", "", k_out_name)]] = parms.all[[k_out_name]] / k_1output

+        fit$ff[[sub("k_", "", k_out_name)]] = parms.all[[k_out_name]] / k_1output

       }

+

       # Return the eigenvalues for comparison with DFOP rate constants

       fit$SFORB[[paste(obs_var, "b1", sep="_")]] = b1

       fit$SFORB[[paste(obs_var, "b2", sep="_")]] = b2

@@ -310,6 +315,7 @@ mkinfit <- function(mkinmod, observed,
   fit$atol <- atol

   fit$parms.all <- parms.all # Return all backtransformed parameters for summary

   fit$odeparms.final <- parms.all[mkinmod$parms] # Return ode parameters for further fitting

+  fit$date <- date()

 

   class(fit) <- c("mkinfit", "modFit") 

   return(fit)

@@ -337,7 +343,12 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, ...) {
   param <- cbind(param, se)

   dimnames(param) <- list(pnames, c("Estimate", "Std. Error"))

 

-  ans <- list(residuals = object$residuals,

+  ans <- list(

+          version = as.character(packageVersion("mkin")),

+          Rversion = paste(R.version$major, R.version$minor, sep="."),

+	  date.fit = object$date,

+	  date.summary = date(),

+          residuals = object$residuals,

           residualVariance = resvar,

           sigma = sqrt(resvar),

           modVariance = modVariance,

@@ -354,6 +365,7 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, ...) {
   ans$fixed <- object$fixed

   ans$errmin <- object$errmin 

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

+  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") 

@@ -362,6 +374,11 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, ...) {
 

 # Expanded from print.summary.modFit

 print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), ...) {

+  cat("mkin version:   ", x$version, "\n")

+  cat("R version:      ", x$Rversion, "\n")

+  cat("Date of fit:    ", x$date.fit, "\n")

+  cat("Date of summary:", x$date.summary, "\n")

+

   cat("\nEquations:\n")

   print(noquote(as.character(x[["diffs"]])))

   df  <- x$df

@@ -393,6 +410,12 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), .
     print(x$distimes, digits=digits,...)

   }    

 

+  printff <- !is.null(x$ff)

+  if(printff){

+    cat("\nEstimated formation fractions:\n")

+    print(data.frame(ff = x$ff), digits=digits,...)

+  }    

+

   printSFORB <- !is.null(x$SFORB)

   if(printSFORB){

     cat("\nEstimated Eigenvalues of SFORB model(s):\n")

diff --git a/R/mkinpredict.R b/R/mkinpredict.R
index 2b7e51d..be43b0e 100644
--- a/R/mkinpredict.R
+++ b/R/mkinpredict.R
@@ -1,3 +1,23 @@
+# $Id$
+
+# Copyright (C) 2010-2012 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/>
+
 mkinpredict <- function(mkinmod, odeparms, odeini, outtimes, solution_type = "deSolve", map_output = TRUE, atol = 1e-6, ...) {
 
   # Get the names of the state variables in the model
@@ -35,7 +55,7 @@ mkinpredict <- function(mkinmod, odeparms, odeini, outtimes, solution_type = "de
           evalparse(parent.name),
           evalparse(paste("k", parent.name, "bound", sep="_")),
           evalparse(paste("k", sub("free", "bound", parent.name), "free", sep="_")),
-          evalparse(paste("k", parent.name, sep="_")))
+          evalparse(paste("k", parent.name, "sink", sep="_")))
     )
     out <- cbind(outtimes, o)
     dimnames(out) <- list(outtimes, c("time", sub("_free", "", parent.name)))