- 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

15 files changed, 443 insertions, 171 deletions

diff --git a/R/mkinfit.R b/R/mkinfit.R
index cf58a7d2..6e455e18 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 2b7e51d5..be43b0e4 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)))
diff --git a/TODO b/TODO
index a49957d2..a7503bda 100644
--- a/TODO
+++ b/TODO
@@ -1,10 +1,9 @@
 Must have:
 - Fix remaining problem with the schaefer complex test data (test.R)
-- Fix parent only SFORB model (see vignette)
-- Fix unit tests for mkinmod
 - Adapt mkinplot function
 
 Nice to have:
+- Make unit tests for mkinpredict, comparing different solution types
 - Calculate confidence intervals for parameters assuming normal distribution
 - Calculate confidence intervals for DT50 and DT90 values when only one parameter is involved
 - Add unit tests for mkinfit
diff --git a/data/FOCUS_2006_DFOP_ref_A_to_B.rda b/data/FOCUS_2006_DFOP_ref_A_to_B.rda
new file mode 100644
index 00000000..8b9dc9f8
--- /dev/null
+++ b/data/FOCUS_2006_DFOP_ref_A_to_B.rda
diff --git a/data/FOCUS_2006_FOMC_ref_A_to_F.rda b/data/FOCUS_2006_FOMC_ref_A_to_F.rda
new file mode 100644
index 00000000..246a53fa
--- /dev/null
+++ b/data/FOCUS_2006_FOMC_ref_A_to_F.rda
diff --git a/data/FOCUS_2006_HS_ref_A_to_F.rda b/data/FOCUS_2006_HS_ref_A_to_F.rda
new file mode 100644
index 00000000..1e3b3864
--- /dev/null
+++ b/data/FOCUS_2006_HS_ref_A_to_F.rda
diff --git a/data/FOCUS_2006_SFO_ref_A_to_F.rda b/data/FOCUS_2006_SFO_ref_A_to_F.rda
new file mode 100644
index 00000000..411c5368
--- /dev/null
+++ b/data/FOCUS_2006_SFO_ref_A_to_F.rda
diff --git a/inst/unitTests/runit.mkinfit.R b/inst/unitTests/runit.mkinfit.R
index 2a026ce0..a06e4ff1 100644
--- a/inst/unitTests/runit.mkinfit.R
+++ b/inst/unitTests/runit.mkinfit.R
@@ -18,36 +18,243 @@
 # You should have received a copy of the GNU General Public License along with

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

 

+# Test SFO model to a relative tolerance of 1% # {{{

+test.FOCUS_2006_SFO <- function()

+{

+  SFO.1 <- mkinmod(parent = list(type = "SFO"))

+  SFO.2 <- mkinmod(parent = list(type = "SFO"), use_of_ff = "max")

+

+  fit.A.SFO.1 <- mkinfit(SFO.1, FOCUS_2006_A, quiet=TRUE)

+  fit.A.SFO.2 <- mkinfit(SFO.2, FOCUS_2006_A, quiet=TRUE)

+

+  median.A.SFO <- as.numeric(lapply(subset(FOCUS_2006_SFO_ref_A_to_F, dataset == "A", 

+                   c(M0, k, DT50, DT90)), "median"))

+

+  fit.A.SFO.1.r <- as.numeric(c(fit.A.SFO.1$parms.all, fit.A.SFO.1$distimes))

+  dev.A.SFO.1 <- abs(round(100 * ((median.A.SFO - fit.A.SFO.1.r)/median.A.SFO), digits=1))

+  checkIdentical(dev.A.SFO.1 < 1, rep(TRUE, length(dev.A.SFO.1)))

+

+  fit.A.SFO.2.r <- as.numeric(c(fit.A.SFO.2$parms.all, fit.A.SFO.2$distimes))

+  dev.A.SFO.2 <- abs(round(100 * ((median.A.SFO - fit.A.SFO.2.r)/median.A.SFO), digits=1))

+  checkIdentical(dev.A.SFO.2 < 1, rep(TRUE, length(dev.A.SFO.2)))

+

+  fit.C.SFO.1 <- mkinfit(SFO.1, FOCUS_2006_C, quiet=TRUE)

+  fit.C.SFO.2 <- mkinfit(SFO.2, FOCUS_2006_C, quiet=TRUE)

+

+  median.C.SFO <- as.numeric(lapply(subset(FOCUS_2006_SFO_ref_A_to_F, dataset == "C", 

+                   c(M0, k, DT50, DT90)), "median"))

+

+  fit.C.SFO.1.r <- as.numeric(c(fit.C.SFO.1$parms.all, fit.C.SFO.1$distimes))

+  dev.C.SFO.1 <- abs(round(100 * ((median.C.SFO - fit.C.SFO.1.r)/median.C.SFO), digits=1))

+  checkIdentical(dev.C.SFO.1 < 1, rep(TRUE, length(dev.C.SFO.1)))

+

+  fit.C.SFO.2.r <- as.numeric(c(fit.C.SFO.2$parms.all, fit.C.SFO.2$distimes))

+  dev.C.SFO.2 <- abs(round(100 * ((median.C.SFO - fit.C.SFO.2.r)/median.C.SFO), digits=1))

+  checkIdentical(dev.C.SFO.2 < 1, rep(TRUE, length(dev.C.SFO.2)))

+} # }}}

+

+# Test FOMC model to a relative tolerance of 1% {{{

+# See kinfit vignette for a discussion of FOMC fits to FOCUS_2006_A

+# In this case, only M0, DT50 and DT90 are checked

+test.FOCUS_2006_FOMC <- function()

+{

+  FOMC <- mkinmod(parent = list(type = "FOMC"))

+

+  # FOCUS_2006_A (compare kinfit vignette for discussion) 

+  fit.A.FOMC <- mkinfit(FOMC, FOCUS_2006_A, quiet=TRUE)

+

+  median.A.FOMC <- as.numeric(lapply(subset(FOCUS_2006_FOMC_ref_A_to_F, dataset == "A", 

+                   c(M0, alpha, beta, DT50, DT90)), "median"))

+

+  fit.A.FOMC.r <- as.numeric(c(fit.A.FOMC$parms.all, fit.A.FOMC$distimes))

+  dev.A.FOMC <- abs(round(100 * ((median.A.FOMC - fit.A.FOMC.r)/median.A.FOMC), digits=1))

+  dev.A.FOMC <- dev.A.FOMC[c(1, 4, 5)]

+  checkIdentical(dev.A.FOMC < 1, rep(TRUE, length(dev.A.FOMC)))

+

+  # FOCUS_2006_B

+  fit.B.FOMC <- mkinfit(FOMC, FOCUS_2006_B, quiet=TRUE)

+

+  median.B.FOMC <- as.numeric(lapply(subset(FOCUS_2006_FOMC_ref_A_to_F, dataset == "B", 

+                   c(M0, alpha, beta, DT50, DT90)), "median"))

+

+  fit.B.FOMC.r <- as.numeric(c(fit.B.FOMC$parms.all, fit.B.FOMC$distimes))

+  dev.B.FOMC <- abs(round(100 * ((median.B.FOMC - fit.B.FOMC.r)/median.B.FOMC), digits=1))

+  dev.B.FOMC <- dev.B.FOMC[c(1, 4, 5)]

+  checkIdentical(dev.B.FOMC < 1, rep(TRUE, length(dev.B.FOMC)))

+

+  # FOCUS_2006_C

+  fit.C.FOMC <- mkinfit(FOMC, FOCUS_2006_C, quiet=TRUE)

+

+  median.C.FOMC <- as.numeric(lapply(subset(FOCUS_2006_FOMC_ref_A_to_F, dataset == "C", 

+                   c(M0, alpha, beta, DT50, DT90)), "median"))

+

+  fit.C.FOMC.r <- as.numeric(c(fit.C.FOMC$parms.all, fit.C.FOMC$distimes))

+  dev.C.FOMC <- abs(round(100 * ((median.C.FOMC - fit.C.FOMC.r)/median.C.FOMC), digits=1))

+  dev.C.FOMC <- dev.C.FOMC[c(1, 4, 5)]

+  checkIdentical(dev.C.FOMC < 1, rep(TRUE, length(dev.C.FOMC)))

+} # }}}

+

+# Test DFOP model, tolerance of 1% with the exception of f parameter for A {{{

+test.FOCUS_2006_DFOP <- function()

+{

+  DFOP <- mkinmod(parent = list(type = "DFOP"))

+

+  # FOCUS_2006_A

+  fit.A.DFOP <- mkinfit(DFOP, FOCUS_2006_A, quiet=TRUE)

+  fit.A.DFOP <- mkinfit(DFOP, FOCUS_2006_A, quiet=TRUE, plot=TRUE)

+

+  median.A.DFOP <- as.numeric(lapply(subset(FOCUS_2006_DFOP_ref_A_to_B, dataset == "A", 

+                   c(M0, k1, k2, f, DT50, DT90)), "median"))

+

+  fit.A.DFOP.r <- as.numeric(c(fit.A.DFOP$parms.all, fit.A.DFOP$distimes))

+  dev.A.DFOP <- abs(round(100 * ((median.A.DFOP - fit.A.DFOP.r)/median.A.DFOP), digits=1))

+  # about 6.7% deviation for parameter f, the others are < 0.1%

+  checkIdentical(dev.A.DFOP < c(1, 1, 1, 10, 1, 1), rep(TRUE, length(dev.A.DFOP)))

+

+  # FOCUS_2006_B

+  fit.B.DFOP <- mkinfit(DFOP, FOCUS_2006_B, quiet=TRUE)

+

+  median.B.DFOP <- as.numeric(lapply(subset(FOCUS_2006_DFOP_ref_A_to_B, dataset == "B", 

+                   c(M0, k1, k2, f, DT50, DT90)), "median"))

+

+  fit.B.DFOP.r <- as.numeric(c(fit.B.DFOP$parms.all, fit.B.DFOP$distimes))

+  dev.B.DFOP <- abs(round(100 * ((median.B.DFOP - fit.B.DFOP.r)/median.B.DFOP), digits=1))

+  # about 0.6% deviation for parameter f, the others are <= 0.1%

+  checkIdentical(dev.B.DFOP < 1, rep(TRUE, length(dev.B.DFOP)))

+} # }}}

+

+# Test HS model to a relative tolerance of 1% excluding Mathematica values {{{

+# as they are unreliable

+test.FOCUS_2006_HS <- function()

+{

+  HS <- mkinmod(parent = list(type = "HS"))

+

+  # FOCUS_2006_A

+  fit.A.HS <- mkinfit(HS, FOCUS_2006_A, quiet=TRUE)

+

+  median.A.HS <- as.numeric(lapply(subset(FOCUS_2006_HS_ref_A_to_F, dataset == "A", 

+                   c(M0, k1, k2, tb, DT50, DT90)), "median"))

+

+  fit.A.HS.r <- as.numeric(c(fit.A.HS$parms.all, fit.A.HS$distimes))

+  dev.A.HS <- abs(round(100 * ((median.A.HS - fit.A.HS.r)/median.A.HS), digits=1))

+  # about 6.7% deviation for parameter f, the others are < 0.1%

+  checkIdentical(dev.A.HS < 1, rep(TRUE, length(dev.A.HS)))

+

+  # FOCUS_2006_B

+  fit.B.HS <- mkinfit(HS, FOCUS_2006_B, quiet=TRUE)

+

+  median.B.HS <- as.numeric(lapply(subset(FOCUS_2006_HS_ref_A_to_F, dataset == "B", 

+                   c(M0, k1, k2, tb, DT50, DT90)), "median"))

+

+  fit.B.HS.r <- as.numeric(c(fit.B.HS$parms.all, fit.B.HS$distimes))

+  dev.B.HS <- abs(round(100 * ((median.B.HS - fit.B.HS.r)/median.B.HS), digits=1))

+  # < 10% deviation for M0, k1, DT50 and DT90, others are problematic

+  dev.B.HS <- dev.B.HS[c(1, 2, 5, 6)]

+  checkIdentical(dev.B.HS < 10, rep(TRUE, length(dev.B.HS)))

+

+  # FOCUS_2006_C

+  fit.C.HS <- mkinfit(HS, FOCUS_2006_C, quiet=TRUE)

+

+  median.C.HS <- as.numeric(lapply(subset(FOCUS_2006_HS_ref_A_to_F, dataset == "C", 

+                   c(M0, k1, k2, tb, DT50, DT90)), "median"))

+

+  fit.A.HS.r <- as.numeric(c(fit.A.HS$parms.all, fit.A.HS$distimes))

+  dev.A.HS <- abs(round(100 * ((median.A.HS - fit.A.HS.r)/median.A.HS), digits=1))

+  # deviation <= 0.1%

+  checkIdentical(dev.A.HS < 1, rep(TRUE, length(dev.A.HS)))

+} # }}}

+

+# Test SFORB model against DFOP solutions to a relative tolerance of 1% # {{{

+test.FOCUS_2006_SFORB <- function()

+{

+  SFORB <- mkinmod(parent = list(type = "SFORB"))

+

+  # FOCUS_2006_A

+  fit.A.SFORB.1 <- mkinfit(SFORB, FOCUS_2006_A, quiet=TRUE)

+  fit.A.SFORB.2 <- mkinfit(SFORB, FOCUS_2006_A, solution_type = "deSolve", quiet=TRUE)

+

+  median.A.SFORB <- as.numeric(lapply(subset(FOCUS_2006_DFOP_ref_A_to_B, dataset == "A", 

+                   c(M0, k1, k2, DT50, DT90)), "median"))

+

+  fit.A.SFORB.1.r <- as.numeric(c(

+                      parent_0 = fit.A.SFORB.1$parms.all[[1]], 

+                      k1 = fit.A.SFORB.1$SFORB[[1]],

+                      k2 = fit.A.SFORB.1$SFORB[[2]],

+                      fit.A.SFORB.1$distimes))

+  dev.A.SFORB.1 <- abs(round(100 * ((median.A.SFORB - fit.A.SFORB.1.r)/median.A.SFORB), digits=1))

+  # The first Eigenvalue is a lot different from k1 in the DFOP fit

+  # The explanation is that the dataset is simply SFO

+  dev.A.SFORB.1 <- dev.A.SFORB.1[c(1, 3, 4, 5)]

+  checkIdentical(dev.A.SFORB.1 < 1, rep(TRUE, length(dev.A.SFORB.1)))

+

+  fit.A.SFORB.2.r <- as.numeric(c(

+                      parent_0 = fit.A.SFORB.2$parms.all[[1]], 

+                      k1 = fit.A.SFORB.2$SFORB[[1]],

+                      k2 = fit.A.SFORB.2$SFORB[[2]],

+                      fit.A.SFORB.2$distimes))

+  dev.A.SFORB.2 <- abs(round(100 * ((median.A.SFORB - fit.A.SFORB.2.r)/median.A.SFORB), digits=1))

+  # The first Eigenvalue is a lot different from k1 in the DFOP fit

+  # The explanation is that the dataset is simply SFO

+  dev.A.SFORB.2 <- dev.A.SFORB.2[c(1, 3, 4, 5)]

+  checkIdentical(dev.A.SFORB.2 < 1, rep(TRUE, length(dev.A.SFORB.2)))

+

+  # FOCUS_2006_B

+  fit.B.SFORB.1 <- mkinfit(SFORB, FOCUS_2006_B, quiet=TRUE)

+  fit.B.SFORB.2 <- mkinfit(SFORB, FOCUS_2006_B, solution_type = "deSolve", quiet=TRUE)

+

+  median.B.SFORB <- as.numeric(lapply(subset(FOCUS_2006_DFOP_ref_A_to_B, dataset == "B", 

+                   c(M0, k1, k2, DT50, DT90)), "median"))

+

+  fit.B.SFORB.1.r <- as.numeric(c(

+                      parent_0 = fit.B.SFORB.1$parms.all[[1]], 

+                      k1 = fit.B.SFORB.1$SFORB[[1]],

+                      k2 = fit.B.SFORB.1$SFORB[[2]],

+                      fit.B.SFORB.1$distimes))

+  dev.B.SFORB.1 <- abs(round(100 * ((median.B.SFORB - fit.B.SFORB.1.r)/median.B.SFORB), digits=1))

+  checkIdentical(dev.B.SFORB.1 < 1, rep(TRUE, length(dev.B.SFORB.1)))

+

+  fit.B.SFORB.2.r <- as.numeric(c(

+                      parent_0 = fit.B.SFORB.2$parms.all[[1]], 

+                      k1 = fit.B.SFORB.2$SFORB[[1]],

+                      k2 = fit.B.SFORB.2$SFORB[[2]],

+                      fit.B.SFORB.2$distimes))

+  dev.B.SFORB.2 <- abs(round(100 * ((median.B.SFORB - fit.B.SFORB.2.r)/median.B.SFORB), digits=1))

+  checkIdentical(dev.B.SFORB.2 < 1, rep(TRUE, length(dev.B.SFORB.2)))

+} # }}}

+

+# Test eigenvalue based fit to Schaefer 2007 data against solution from conference paper {{{

 test.mkinfit.schaefer07_complex_example <- function()

 {

   schaefer07_complex_model <- mkinmod(

-    parent = list(type = "SFO", to = c("A1", "B1", "C1")),

+    parent = list(type = "SFO", to = c("A1", "B1", "C1"), sink = FALSE),

     A1 = list(type = "SFO", to = "A2"),

     B1 = list(type = "SFO"),

     C1 = list(type = "SFO"),

     A2 = list(type = "SFO"))

   

-# Commented out because it takes too much time and is currently not used (see below)

-#  fit <- mkinfit(schaefer07_complex_model, 

-#    mkin_wide_to_long(schaefer07_complex_case, time = "time"))

-#  r <- schaefer07_complex_results

-#  r$mkin <- c(

-#    fit$parms.all["k_parent"],

-#    fit$distimes["parent", "DT50"],

-#    fit$parms.all["f_parent_to_A1"],

-#    fit$parms.all["k_A1"],

-#    fit$distimes["A1", "DT50"],

-#    fit$parms.all["f_parent_to_B1"],

-#    fit$parms.all["k_B1"],

-#    fit$distimes["B1", "DT50"],

-#    fit$parms.all["f_parent_to_C1"],

-#    fit$parms.all["k_C1"],

-#    fit$distimes["C1", "DT50"],

-#    fit$parms.all["f_A1_to_A2"],

-#    fit$parms.all["k_A2"],

-#    fit$distimes["A2", "DT50"])

-#  r$means <- (r$KinGUI + r$ModelMaker)/2

-#  r$mkin.deviation <- abs(round(100 * ((r$mkin - r$means)/r$means), digits=1))

-  # Commented out the check as mkin is fitting a different model

-  #checkIdentical(r$mkin.deviation < 10, rep(TRUE, length(r$mkin.deviation)))

-}

+  fit <- mkinfit(schaefer07_complex_model, 

+    mkin_wide_to_long(schaefer07_complex_case, time = "time"))

+  s <- summary(fit)

+  r <- schaefer07_complex_results

+  attach(as.list(fit$parms.all))

+  k_parent <- sum(k_parent_A1, k_parent_B1, k_parent_C1)

+  r$mkin <- c(

+    k_parent,

+    s$distimes["parent", "DT50"],

+    s$ff["parent_A1"],

+    sum(k_A1_sink, k_A1_A2),

+    s$distimes["A1", "DT50"],

+    s$ff["parent_B1"],

+    k_B1_sink,

+    s$distimes["B1", "DT50"],

+    s$ff["parent_C1"],

+    k_C1_sink,

+    s$distimes["C1", "DT50"],

+    s$ff["A1_A2"],

+    k_A2_sink,

+    s$distimes["A2", "DT50"])

+  r$means <- (r$KinGUI + r$ModelMaker)/2

+  r$mkin.deviation <- abs(round(100 * ((r$mkin - r$means)/r$means), digits=1))

+  checkIdentical(r$mkin.deviation < 10, rep(TRUE, length(r$mkin.deviation)))

+} # }}}

+# vim: set foldmethod=marker ts=2 sw=2 expandtab:

diff --git a/inst/unitTests/runit.mkinmod.R b/inst/unitTests/runit.mkinmod.R
deleted file mode 100644
index 38b07cd5..00000000
--- a/inst/unitTests/runit.mkinmod.R
+++ /dev/null
@@ -1,134 +0,0 @@
-# $Id: runit.mkinmod.R 64 2010-09-01 13:33:51Z jranke $

-

-# 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/>

-

-test.mkinmod.SFO <- function()

-{  

-  SFO.diffs <- c(

-    parent = "d_parent = - k_parent * parent"

-  )

-  SFO.parms <- c("k_parent")

-  SFO.map <- list(parent = c(SFO = "parent"))

-  SFO.coefmat <- matrix("- k_parent", dimnames = list("parent", "parent"))

-  SFO <- list(diffs = SFO.diffs, parms = SFO.parms, map = SFO.map, 

-    coefmat = SFO.coefmat)

-  class(SFO) <- "mkinmod"

-  SFO.mkinmod <- mkinmod(

-    parent = list(type = "SFO")

-  )

-  checkIdentical(SFO, SFO.mkinmod)

-}

-

-test.mkinmod.SFORB <- function()

-{  

-  SFORB.diffs <- c(

-    parent_free = paste(

-      "d_parent_free = - k_parent_free * parent_free", 

-        "- k_parent_free_bound * parent_free",

-        "+ k_parent_bound_free * parent_bound"),

-    parent_bound = paste(

-      "d_parent_bound =",

-        "+ k_parent_free_bound * parent_free",

-        "- k_parent_bound_free * parent_bound")

-  )

-  SFORB.parms <- c("k_parent_free", "k_parent_free_bound", "k_parent_bound_free")

-  SFORB.map <- list(parent = c(SFORB = "parent_free", SFORB = "parent_bound"))

-  vars <- paste("parent", c("free", "bound"), sep="_")

-  SFORB.coefmat <- matrix(

-    c("- k_parent_free - k_parent_free_bound", "k_parent_bound_free",

-      "k_parent_free_bound", "- k_parent_bound_free"), nrow=2, byrow=TRUE, 

-    dimnames=list(vars, vars))

-  SFORB <- list(diffs = SFORB.diffs, parms = SFORB.parms, 

-    map = SFORB.map, coefmat = SFORB.coefmat)

-  class(SFORB) <- "mkinmod"

-  SFORB.mkinmod <- mkinmod(

-    parent = list(type = "SFORB")

-  )

-  #checkIdentical(SFORB, SFORB.mkinmod)

-}

-

-test.mkinmod.SFO_SFO <- function()

-{  

-  SFO_SFO.diffs <- c(

-    parent = "d_parent = - k_parent * parent",

-    m1 = "d_m1 = + f_parent_to_m1 * k_parent * parent - k_m1 * m1"

-  )

-  SFO_SFO.parms <- c("k_parent", "f_parent_to_m1", "k_m1")

-  SFO_SFO.map <- list(parent = c(SFO = "parent"), m1 = c(SFO = "m1"))

-  vars <- c("parent", "m1")

-  SFO_SFO.coefmat <- matrix(c("- k_parent", 

-          "0", "f_parent_to_m1 * k_parent", "- k_m1"), nrow=2, byrow=TRUE,

-      dimnames=list(vars, vars))

-  SFO_SFO <- list(diffs = SFO_SFO.diffs, parms = SFO_SFO.parms, 

-    map = SFO_SFO.map, coefmat = SFO_SFO.coefmat)

-  class(SFO_SFO) <- "mkinmod"

-  SFO_SFO.mkinmod <- mkinmod(

-    parent = list(type = "SFO", to = "m1"),

-    m1 = list(type = "SFO", sink=TRUE)

-  )

-  checkIdentical(SFO_SFO, SFO_SFO.mkinmod)

-}

-

-test.mkinmod.SFO_SFO2 <- function()

-{  

-  SFO_SFO2.diffs <- c(

-    parent = "d_parent = - k_parent * parent",

-    m1 = "d_m1 = + f_parent_to_m1 * k_parent * parent - k_m1 * m1",

-    m2 = "d_m2 = + f_parent_to_m2 * k_parent * parent - k_m2 * m2"

-  )

-  SFO_SFO2.parms <- c("k_parent", "f_parent_to_m1", "f_parent_to_m2", "k_m1", "k_m2")

-  SFO_SFO2.map <- list(parent = c(SFO = "parent"), m1 = c(SFO = "m1"), m2 = c(SFO = "m2"))

-  vars <- c("parent", "m1", "m2")

-  SFO_SFO2.coefmat <- matrix(

-      c("- k_parent", "0", "0",

-          "f_parent_to_m1 * k_parent", "- k_m1", "0",

-          "f_parent_to_m2 * k_parent", "0", "- k_m2"), nrow=3, byrow=TRUE,

-      dimnames=list(vars, vars))

-  SFO_SFO2 <- list(diffs = SFO_SFO2.diffs, parms = SFO_SFO2.parms, 

-      map = SFO_SFO2.map, coefmat = SFO_SFO2.coefmat)

-  class(SFO_SFO2) <- "mkinmod"

-  SFO_SFO2.mkinmod <- mkinmod(

-    parent = list(type = "SFO", to = c("m1", "m2"), sink=TRUE),

-    m1 = list(type = "SFO", sink=TRUE),

-    m2 = list(type = "SFO", sink=TRUE)

-  )

-  checkIdentical(SFO_SFO2, SFO_SFO2.mkinmod)

-}

-

-test.mkinmod.FOMC_SFO2 <- function()

-{  

-  FOMC_SFO2.diffs <- c(

-    parent = "d_parent = - (alpha/beta) * ((time/beta) + 1)^-1 * parent",

-    m1 = "d_m1 = + f_parent_to_m1 * (alpha/beta) * ((time/beta) + 1)^-1 * parent - k_m1 * m1",

-    m2 = "d_m2 = + f_parent_to_m2 * (alpha/beta) * ((time/beta) + 1)^-1 * parent - k_m2 * m2"

-  )

-  FOMC_SFO2.parms <- c("alpha", "beta", "f_parent_to_m1", "f_parent_to_m2", "k_m1", "k_m2")

-  FOMC_SFO2.map <- list(parent = c(FOMC = "parent"), 

-    m1 = c(SFO = "m1"), 

-    m2 = c(SFO = "m2"))

-  FOMC_SFO2 <- list(diffs = FOMC_SFO2.diffs, parms = FOMC_SFO2.parms, 

-    map = FOMC_SFO2.map)

-  class(FOMC_SFO2) <- "mkinmod"

-  FOMC_SFO2.mkinmod <- mkinmod(

-    parent = list(type = "FOMC", to = c("m1", "m2"), sink=TRUE),

-    m1 = list(type = "SFO"),

-    m2 = list(type = "SFO")

-  )

-  checkIdentical(FOMC_SFO2, FOMC_SFO2.mkinmod)

-}

diff --git a/man/FOCUS_2006_DFOP_ref_A_to_B.Rd b/man/FOCUS_2006_DFOP_ref_A_to_B.Rd
new file mode 100644
index 00000000..88bd4ac6
--- /dev/null
+++ b/man/FOCUS_2006_DFOP_ref_A_to_B.Rd
@@ -0,0 +1,39 @@
+\name{FOCUS_2006_DFOP_ref_A_to_B}

+\Rdversion{1.1}

+\alias{FOCUS_2006_DFOP_ref_A_to_B}

+\docType{data}

+\title{

+Results of fitting the DFOP model to Datasets A to B of FOCUS (2006)

+}

+\description{

+A table with the fitted parameters and the resulting DT50 and DT90 values

+generated with different software packages. Taken directly from FOCUS (2006).

+The results from fitting the data with the Topfit software was removed, as

+the initial concentration of the parent compound was fixed to a value of 100

+in this fit.

+}

+\usage{data(FOCUS_2006_DFOP_ref_A_to_B)}

+\format{

+  A data frame containing the following variables.

+  \describe{

+    \item{\code{package}}{a factor giving the name of the software package}

+    \item{\code{M0}}{The fitted initial concentration of the parent compound}

+    \item{\code{f}}{The fitted f parameter}

+    \item{\code{k1}}{The fitted k1 parameter}

+    \item{\code{k2}}{The fitted k2 parameter}

+    \item{\code{DT50}}{The resulting half-life of the parent compound}

+    \item{\code{DT90}}{The resulting DT90 of the parent compound}

+    \item{\code{dataset}}{The FOCUS dataset that was used}

+  }

+}

+\source{

+  FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence and

+  Degradation Kinetics from Environmental Fate Studies on Pesticides in EU

+  Registration} Report of the FOCUS Work Group on Degradation Kinetics,

+  EC Document Reference Sanco/10058/2005 version 2.0, 434 pp,

+  \url{http://focus.jrc.ec.europa.eu/dk} 

+}

+\examples{

+data(FOCUS_2006_DFOP_ref_A_to_B)

+}

+\keyword{datasets}

diff --git a/man/FOCUS_2006_FOMC_ref_A_to_F.Rd b/man/FOCUS_2006_FOMC_ref_A_to_F.Rd
new file mode 100644
index 00000000..2fcc2db7
--- /dev/null
+++ b/man/FOCUS_2006_FOMC_ref_A_to_F.Rd
@@ -0,0 +1,38 @@
+\name{FOCUS_2006_FOMC_ref_A_to_F}

+\Rdversion{1.1}

+\alias{FOCUS_2006_FOMC_ref_A_to_F}

+\docType{data}

+\title{

+Results of fitting the FOMC model to Datasets A to F of FOCUS (2006)

+}

+\description{

+A table with the fitted parameters and the resulting DT50 and DT90 values

+generated with different software packages. Taken directly from FOCUS (2006).

+The results from fitting the data with the Topfit software was removed, as

+the initial concentration of the parent compound was fixed to a value of 100

+in this fit.

+}

+\usage{data(FOCUS_2006_FOMC_ref_A_to_F)}

+\format{

+  A data frame containing the following variables.

+  \describe{

+    \item{\code{package}}{a factor giving the name of the software package}

+    \item{\code{M0}}{The fitted initial concentration of the parent compound}

+    \item{\code{alpha}}{The fitted alpha parameter}

+    \item{\code{beta}}{The fitted beta parameter}

+    \item{\code{DT50}}{The resulting half-life of the parent compound}

+    \item{\code{DT90}}{The resulting DT90 of the parent compound}

+    \item{\code{dataset}}{The FOCUS dataset that was used}

+  }

+}

+\source{

+  FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence and

+  Degradation Kinetics from Environmental Fate Studies on Pesticides in EU

+  Registration} Report of the FOCUS Work Group on Degradation Kinetics,

+  EC Document Reference Sanco/10058/2005 version 2.0, 434 pp,

+  \url{http://focus.jrc.ec.europa.eu/dk} 

+}

+\examples{

+data(FOCUS_2006_FOMC_ref_A_to_F)

+}

+\keyword{datasets}

diff --git a/man/FOCUS_2006_HS_ref_A_to_F.Rd b/man/FOCUS_2006_HS_ref_A_to_F.Rd
new file mode 100644
index 00000000..6fc99933
--- /dev/null
+++ b/man/FOCUS_2006_HS_ref_A_to_F.Rd
@@ -0,0 +1,39 @@
+\name{FOCUS_2006_HS_ref_A_to_F}

+\Rdversion{1.1}

+\alias{FOCUS_2006_HS_ref_A_to_F}

+\docType{data}

+\title{

+Results of fitting the HS model to Datasets A to F of FOCUS (2006)

+}

+\description{

+A table with the fitted parameters and the resulting DT50 and DT90 values

+generated with different software packages. Taken directly from FOCUS (2006).

+The results from fitting the data with the Topfit software was removed, as

+the initial concentration of the parent compound was fixed to a value of 100

+in this fit.

+}

+\usage{data(FOCUS_2006_HS_ref_A_to_F)}

+\format{

+  A data frame containing the following variables.

+  \describe{

+    \item{\code{package}}{a factor giving the name of the software package}

+    \item{\code{M0}}{The fitted initial concentration of the parent compound}

+    \item{\code{tb}}{The fitted tb parameter}

+    \item{\code{k1}}{The fitted k1 parameter}

+    \item{\code{k2}}{The fitted k2 parameter}

+    \item{\code{DT50}}{The resulting half-life of the parent compound}

+    \item{\code{DT90}}{The resulting DT90 of the parent compound}

+    \item{\code{dataset}}{The FOCUS dataset that was used}

+  }

+}

+\source{

+  FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence and

+  Degradation Kinetics from Environmental Fate Studies on Pesticides in EU

+  Registration} Report of the FOCUS Work Group on Degradation Kinetics,

+  EC Document Reference Sanco/10058/2005 version 2.0, 434 pp,

+  \url{http://focus.jrc.ec.europa.eu/dk} 

+}

+\examples{

+data(FOCUS_2006_HS_ref_A_to_F)

+}

+\keyword{datasets}

diff --git a/man/FOCUS_2006_SFO_ref_A_to_F.Rd b/man/FOCUS_2006_SFO_ref_A_to_F.Rd
new file mode 100644
index 00000000..19650ed4
--- /dev/null
+++ b/man/FOCUS_2006_SFO_ref_A_to_F.Rd
@@ -0,0 +1,37 @@
+\name{FOCUS_2006_SFO_ref_A_to_F}

+\Rdversion{1.1}

+\alias{FOCUS_2006_SFO_ref_A_to_F}

+\docType{data}

+\title{

+Results of fitting the SFO model to Datasets A to F of FOCUS (2006)

+}

+\description{

+A table with the fitted parameters and the resulting DT50 and DT90 values

+generated with different software packages. Taken directly from FOCUS (2006).

+The results from fitting the data with the Topfit software was removed, as

+the initial concentration of the parent compound was fixed to a value of 100

+in this fit.

+}

+\usage{data(FOCUS_2006_SFO_ref_A_to_F)}

+\format{

+  A data frame containing the following variables.

+  \describe{

+    \item{\code{package}}{a factor giving the name of the software package}

+    \item{\code{M0}}{The fitted initial concentration of the parent compound}

+    \item{\code{k}}{The fitted first-order degradation rate constant}

+    \item{\code{DT50}}{The resulting half-life of the parent compound}

+    \item{\code{DT90}}{The resulting DT90 of the parent compound}

+    \item{\code{dataset}}{The FOCUS dataset that was used}

+  }

+}

+\source{

+  FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence and

+  Degradation Kinetics from Environmental Fate Studies on Pesticides in EU

+  Registration} Report of the FOCUS Work Group on Degradation Kinetics,

+  EC Document Reference Sanco/10058/2005 version 2.0, 434 pp,

+  \url{http://focus.jrc.ec.europa.eu/dk} 

+}

+\examples{

+data(FOCUS_2006_SFO_ref_A_to_F)

+}

+\keyword{datasets}

diff --git a/man/mkinfit.Rd b/man/mkinfit.Rd
index 0c8e48f1..5f70dae9 100644
--- a/man/mkinfit.Rd
+++ b/man/mkinfit.Rd
@@ -17,7 +17,8 @@ mkinfit(mkinmod, observed,
   solution_type = "auto",
   plot = FALSE, quiet = FALSE, err = NULL, weight = "none", 
   scaleVar = FALSE, 
-  atol = 1e-6, n.outtimes, ...)
+  atol = 1e-6, n.outtimes, 
+  trace_parms, ...)
 }
 \arguments{
   \item{mkinmod}{
@@ -92,6 +93,9 @@ mkinfit(mkinmod, observed,
     the numerical solver if that is used (see \code{solution} argument. 
     The default value is 100.
   }
+  \item{trace_parms}{
+    Should a trace of the parameter values be listed?
+  }
   \item{\dots}{
     Further arguments that will be passed to \code{\link{modFit}}. 
   }
@@ -101,7 +105,7 @@ mkinfit(mkinmod, observed,
     A summary can be obtained by \code{\link{summary.mkinfit}}. 
 }
 \author{
-  Johannes Ranke <jranke@{harlan.com,uni-bremen.de}>
+  Johannes Ranke <jranke@uni-bremen.de>
 }
 \examples{
 # One parent compound, one metabolite, both single first order.
diff --git a/vignettes/mkin.Rnw b/vignettes/mkin.Rnw
index ee472ba5..f674db7e 100644
--- a/vignettes/mkin.Rnw
+++ b/vignettes/mkin.Rnw
@@ -147,8 +147,8 @@ on the fundamental system of the coefficient matrix, as proposed by
 <<model_fitting, echo=TRUE>>=

 SFO.fit <- mkinfit(SFO, FOCUS_2006_C)

 summary(SFO.fit)

-#SFORB.fit <- mkinfit(SFORB, FOCUS_2006_C)

-#summary(SFORB.fit)

+SFORB.fit <- mkinfit(SFORB, FOCUS_2006_C)

+summary(SFORB.fit)

 SFO_SFO.fit <- mkinfit(SFO_SFO, FOCUS_2006_D, plot=TRUE)

 summary(SFO_SFO.fit, data=FALSE)

 SFORB_SFO.fit <- mkinfit(SFORB_SFO, FOCUS_2006_D, plot=TRUE)