Bug fix and unit tests for mkinerrmin

See NEWS.md for details

9 files changed, 146 insertions, 98 deletions

diff --git a/NEWS.md b/NEWS.md
index d05c2095..1ed94c97 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -4,12 +4,19 @@
 
 - The initial value (state.ini) for the observed variable with the highest observed residue is set to 100 in case it has no time zero observation and `state.ini = "auto"`
 
+- A basic unit test for `mkinerrmin()` was written
+
 ## BUG FIXES
 
 - `mkinfit()`: The internally fitted parameter for `g` was named `g_ilr` even when `transform_fractions=FALSE`
 
 - `mkinfit()`: The initial value (state.ini) for the parent compound was not set when the parent was not the (only) variable with the highest value in the observed data.
 
+- `mkinerrmin()`: When checking for degrees of freedom for metabolites, check
+  if their time zero value is fixed instead of checking if the observed value
+  is zero. This ensures correct calculation of degrees of freedom also in cases
+  where the metabolite residue at time zero is greater zero.
+
 ## MINOR CHANGES
 
 - The formatting of differential equations in the summary was improved by wrapping overly long lines
diff --git a/R/mkinerrmin.R b/R/mkinerrmin.R
index 09724730..2697d0a0 100644
--- a/R/mkinerrmin.R
+++ b/R/mkinerrmin.R
@@ -36,10 +36,11 @@ mkinerrmin <- function(fit, alpha = 0.05)
     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))

+  # Remove values at time zero for variables whose value for state.ini is fixed,

+  # as these will not have any effect in the optimization and should therefore not 

+  # be counted as degrees of freedom.

+  fixed_initials = gsub("_0$", "", rownames(subset(fit$fixed, type = "state")))

+  errdata <- subset(errdata, !(time == 0 & name %in% fixed_initials))

 

   n.optim.overall <- length(parms.optim)

 

diff --git a/TODO b/TODO
index 92a91069..f979d13a 100644
--- a/TODO
+++ b/TODO
@@ -2,11 +2,6 @@ TODO for version 1.0
 - Think about what a user would expect from version 1.0
 - Complete the main package vignette named mkin to include a method description
 - Improve order of parameters in output
-- Write unit tests for mkinerrmin
-- When checking for degrees of freedom for metabolites, check if their time
-  zero value (state.ini) is fixed instead of checking if the observed value is
-  zero (usually in regulatory kinetics it is set to zero anyway, but in the
-  case of known impurities this may not be the case).
 
 Nice to have:
 - Get starting values for formation fractions from data
diff --git a/inst/unitTests/runit.mkinerrmin.R b/inst/unitTests/runit.mkinerrmin.R
new file mode 100644
index 00000000..56a33ff9
--- /dev/null
+++ b/inst/unitTests/runit.mkinerrmin.R
@@ -0,0 +1,62 @@
+# Test SFO_SFO model with FOCUS_2006_D against Schaefer 2007 paper, tolerance = 1% # {{{
+# and check chi2 error values against values obtained with mkin 0.33
+test.FOCUS_2006_D_SFO_SFO <- function()
+{
+  SFO_SFO.1 <- mkinmod(parent = list(type = "SFO", to = "m1"),
+         m1 = list(type = "SFO"), use_of_ff = "min")
+  SFO_SFO.2 <- mkinmod(parent = list(type = "SFO", to = "m1"),
+         m1 = list(type = "SFO"), use_of_ff = "max")
+
+  fit.1.e <- mkinfit(SFO_SFO.1, FOCUS_2006_D)
+  fit.1.d <- mkinfit(SFO_SFO.1, solution_type = "deSolve", FOCUS_2006_D)
+  fit.2.e <- mkinfit(SFO_SFO.2, FOCUS_2006_D)
+  fit.2.d <- mkinfit(SFO_SFO.2, solution_type = "deSolve", FOCUS_2006_D)
+
+  FOCUS_2006_D_results_schaefer07_means <- c(
+    parent_0 = 99.65, DT50_parent = 7.04, DT50_m1 = 131.34)
+
+  r.1.e <- c(fit.1.e$bparms.optim[[1]], endpoints(fit.1.e)$distimes[[1]])
+  r.1.d <- c(fit.1.d$bparms.optim[[1]], endpoints(fit.1.d)$distimes[[1]])
+  r.2.e <- c(fit.2.e$bparms.optim[[1]], endpoints(fit.2.e)$distimes[[1]])
+  r.2.d <- c(fit.2.d$bparms.optim[[1]], endpoints(fit.2.d)$distimes[[1]])
+
+  dev.1.e <- 100 * (r.1.e - FOCUS_2006_D_results_schaefer07_means)/r.1.e 
+  checkIdentical(as.numeric(abs(dev.1.e)) < 1, rep(TRUE, 3))
+  dev.1.d <- 100 * (r.1.d - FOCUS_2006_D_results_schaefer07_means)/r.1.d 
+  checkIdentical(as.numeric(abs(dev.1.d)) < 1, rep(TRUE, 3))
+  dev.2.e <- 100 * (r.2.e - FOCUS_2006_D_results_schaefer07_means)/r.2.e 
+  checkIdentical(as.numeric(abs(dev.2.e)) < 1, rep(TRUE, 3))
+  dev.2.d <- 100 * (r.2.d - FOCUS_2006_D_results_schaefer07_means)/r.2.d 
+  checkIdentical(as.numeric(abs(dev.2.d)) < 1, rep(TRUE, 3))
+
+  round(mkinerrmin(fit.2.e), 4)
+  round(mkinerrmin(fit.2.d), 4)
+
+  errmin.FOCUS_2006_D_rounded = data.frame(
+    err.min = c(0.0640, 0.0646, 0.0469),
+    n.optim = c(4, 2, 2),
+    df = c(15, 7, 8), 
+    row.names = c("All data", "parent", "m1"))
+  checkEqualsNumeric(round(mkinerrmin(fit.2.e), 4),
+                     errmin.FOCUS_2006_D_rounded)
+} # }}}
+
+# Test SFO_SFO model with FOCUS_2006_E against values obtained with mkin 0.33 {{{
+test.FOCUS_2006_E_SFO_SFO <- function()
+{
+  SFO_SFO.2 <- mkinmod(parent = list(type = "SFO", to = "m1"),
+         m1 = list(type = "SFO"), use_of_ff = "max")
+
+  fit.2.e <- mkinfit(SFO_SFO.2, FOCUS_2006_E)
+
+  round(mkinerrmin(fit.2.e), 4)
+  errmin.FOCUS_2006_E_rounded = data.frame(
+    err.min = c(0.1544, 0.1659, 0.1095),
+    n.optim = c(4, 2, 2),
+    df = c(13, 7, 6),
+    row.names = c("All data", "parent", "m1"))
+  checkEqualsNumeric(round(mkinerrmin(fit.2.e), 4),
+                     errmin.FOCUS_2006_E_rounded)
+} # }}}
+
+
diff --git a/inst/unitTests/runit.mkinfit.R b/inst/unitTests/runit.mkinfit.R
index fdbc86e0..8eefb995 100644
--- a/inst/unitTests/runit.mkinfit.R
+++ b/inst/unitTests/runit.mkinfit.R
@@ -1,6 +1,4 @@
-# $Id: runit.mkinfit.R 68 2010-09-09 22:40:04Z jranke $

-

-# Copyright (C) 2010-2013 Johannes Ranke

+# Copyright (C) 2010-2014 Johannes Ranke

 # Contact: jranke@uni-bremen.de

 

 # This file is part of the R package mkin

@@ -189,40 +187,6 @@ test.FOCUS_2006_SFORB <- function()
   checkIdentical(dev.B.SFORB.2 < 1, rep(TRUE, length(dev.B.SFORB.2)))

 } # }}}

 

-# Test SFO_SFO model with FOCUS_2006_D against Schaefer 2007 paper, tolerance = 1% # {{{

-test.FOCUS_2006_D_SFO_SFO <- function()

-{

-  SFO_SFO.1 <- mkinmod(parent = list(type = "SFO", to = "m1"),

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

-  SFO_SFO.2 <- mkinmod(parent = list(type = "SFO", to = "m1"),

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

-

-  fit.1.e <- mkinfit(SFO_SFO.1, FOCUS_2006_D)

-  fit.1.d <- mkinfit(SFO_SFO.1, solution_type = "deSolve", FOCUS_2006_D)

-  fit.2.e <- mkinfit(SFO_SFO.2, FOCUS_2006_D)

-  SFO <- mkinmod(parent = list(type = "SFO"))

-  f.SFO <- mkinfit(SFO, FOCUS_2006_D)

-  fit.2.d <- mkinfit(SFO_SFO.2, solution_type = "deSolve", FOCUS_2006_D)

-  fit.2.e <- mkinfit(SFO_SFO.2, FOCUS_2006_D)

-

-  FOCUS_2006_D_results_schaefer07_means <- c(

-    parent_0 = 99.65, DT50_parent = 7.04, DT50_m1 = 131.34)

-

-  r.1.e <- c(fit.1.e$bparms.optim[[1]], endpoints(fit.1.e)$distimes[[1]])

-  r.1.d <- c(fit.1.d$bparms.optim[[1]], endpoints(fit.1.d)$distimes[[1]])

-  r.2.e <- c(fit.2.e$bparms.optim[[1]], endpoints(fit.2.e)$distimes[[1]])

-  r.2.d <- c(fit.2.d$bparms.optim[[1]], endpoints(fit.2.d)$distimes[[1]])

-

-  dev.1.e <- 100 * (r.1.e - FOCUS_2006_D_results_schaefer07_means)/r.1.e 

-  checkIdentical(as.numeric(abs(dev.1.e)) < 1, rep(TRUE, 3))

-  dev.1.d <- 100 * (r.1.d - FOCUS_2006_D_results_schaefer07_means)/r.1.d 

-  checkIdentical(as.numeric(abs(dev.1.d)) < 1, rep(TRUE, 3))

-  dev.2.e <- 100 * (r.2.e - FOCUS_2006_D_results_schaefer07_means)/r.2.e 

-  checkIdentical(as.numeric(abs(dev.2.e)) < 1, rep(TRUE, 3))

-  dev.2.d <- 100 * (r.2.d - FOCUS_2006_D_results_schaefer07_means)/r.2.d 

-  checkIdentical(as.numeric(abs(dev.2.d)) < 1, rep(TRUE, 3))

-} # }}}

-

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

 test.mkinfit.schaefer07_complex_example <- function()

 {

diff --git a/man/mkinerrmin.Rd b/man/mkinerrmin.Rd
index c43d87a1..78ab414e 100644
--- a/man/mkinerrmin.Rd
+++ b/man/mkinerrmin.Rd
@@ -34,6 +34,16 @@ mkinerrmin(fit, alpha = 0.05)
 \details{

     This function is used internally by \code{\link{summary.mkinfit}}.

 }

+\examples{

+SFO_SFO = mkinmod(parent = list(type = "SFO", to = "m1"),

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

+                  use_of_ff = "max")

+

+fit_FOCUS_D = mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE)

+round(mkinerrmin(fit_FOCUS_D), 4)

+fit_FOCUS_E = mkinfit(SFO_SFO, FOCUS_2006_E, quiet = TRUE)

+round(mkinerrmin(fit_FOCUS_E), 4)

+}

 \references{ 

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

   Degradation Kinetics from Environmental Fate Studies on Pesticides in EU

diff --git a/tests/doRUnit.R b/tests/doRUnit.R
index f0f82812..9faee940 100644
--- a/tests/doRUnit.R
+++ b/tests/doRUnit.R
@@ -1,4 +1,3 @@
-# $Id: doRUnit.R 96 2011-04-29 11:10:40Z jranke $
 # Adapted from a version around 2.9 of the rcdk package by Rajarshi Guha
 if(require("RUnit", quietly=TRUE)) {
  
diff --git a/vignettes/FOCUS_L.html b/vignettes/FOCUS_L.html
index ab7ccaee..2dd186de 100644
--- a/vignettes/FOCUS_L.html
+++ b/vignettes/FOCUS_L.html
@@ -193,7 +193,13 @@ hr {
 report, p. 284:</p>
 
 <pre><code class="r">library(&quot;mkin&quot;)
-FOCUS_2006_L1 = data.frame(
+</code></pre>
+
+<pre><code>## Loading required package: minpack.lm
+## Loading required package: rootSolve
+</code></pre>
+
+<pre><code class="r">FOCUS_2006_L1 = data.frame(
   t = rep(c(0, 1, 2, 3, 5, 7, 14, 21, 30), each = 2),
   parent = c(88.3, 91.4, 85.6, 84.5, 78.9, 77.6, 
              72.0, 71.9, 50.3, 59.4, 47.0, 45.1,
@@ -215,17 +221,17 @@ given in the FOCUS report. </p>
 summary(m.L1.SFO)
 </code></pre>
 
-<pre><code>## mkin version:    0.9.32 
+<pre><code>## mkin version:    0.9.33 
 ## R version:       3.1.1 
-## Date of fit:     Thu Jul 24 10:32:09 2014 
-## Date of summary: Thu Jul 24 10:32:09 2014 
+## Date of fit:     Mon Aug 25 10:34:14 2014 
+## Date of summary: Mon Aug 25 10:34:14 2014 
 ## 
 ## Equations:
-## [1] d_parent = - k_parent_sink * parent
+## d_parent = - k_parent_sink * parent
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted with method Marq using 14 model solutions performed in 0.081 s
+## Fitted with method Marq using 14 model solutions performed in 0.083 s
 ## 
 ## Weighting: none
 ## 
@@ -318,17 +324,17 @@ is checked.</p>
 summary(m.L1.FOMC, data = FALSE)
 </code></pre>
 
-<pre><code>## mkin version:    0.9.32 
+<pre><code>## mkin version:    0.9.33 
 ## R version:       3.1.1 
-## Date of fit:     Thu Jul 24 10:32:10 2014 
-## Date of summary: Thu Jul 24 10:32:11 2014 
+## Date of fit:     Mon Aug 25 10:34:17 2014 
+## Date of summary: Mon Aug 25 10:34:17 2014 
 ## 
 ## Equations:
-## [1] d_parent = - (alpha/beta) * ((time/beta) + 1)^-1 * parent
+## d_parent = - (alpha/beta) * ((time/beta) + 1)^-1 * parent
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted with method Marq using 53 model solutions performed in 0.321 s
+## Fitted with method Marq using 53 model solutions performed in 0.3 s
 ## 
 ## Weighting: none
 ## 
@@ -412,17 +418,17 @@ FOCUS_2006_L2_mkin &lt;- mkin_wide_to_long(FOCUS_2006_L2)
 summary(m.L2.SFO)
 </code></pre>
 
-<pre><code>## mkin version:    0.9.32 
+<pre><code>## mkin version:    0.9.33 
 ## R version:       3.1.1 
-## Date of fit:     Thu Jul 24 10:32:11 2014 
-## Date of summary: Thu Jul 24 10:32:11 2014 
+## Date of fit:     Mon Aug 25 10:34:17 2014 
+## Date of summary: Mon Aug 25 10:34:17 2014 
 ## 
 ## Equations:
-## [1] d_parent = - k_parent_sink * parent
+## d_parent = - k_parent_sink * parent
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted with method Marq using 29 model solutions performed in 0.196 s
+## Fitted with method Marq using 29 model solutions performed in 0.184 s
 ## 
 ## Weighting: none
 ## 
@@ -522,17 +528,17 @@ mkinresplot(m.L2.FOMC)
 <pre><code class="r">summary(m.L2.FOMC, data = FALSE)
 </code></pre>
 
-<pre><code>## mkin version:    0.9.32 
+<pre><code>## mkin version:    0.9.33 
 ## R version:       3.1.1 
-## Date of fit:     Thu Jul 24 10:32:11 2014 
-## Date of summary: Thu Jul 24 10:32:11 2014 
+## Date of fit:     Mon Aug 25 10:34:17 2014 
+## Date of summary: Mon Aug 25 10:34:17 2014 
 ## 
 ## Equations:
-## [1] d_parent = - (alpha/beta) * ((time/beta) + 1)^-1 * parent
+## d_parent = - (alpha/beta) * ((time/beta) + 1)^-1 * parent
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted with method Marq using 35 model solutions performed in 0.223 s
+## Fitted with method Marq using 35 model solutions performed in 0.2 s
 ## 
 ## Weighting: none
 ## 
@@ -608,17 +614,19 @@ plot(m.L2.DFOP)
 <pre><code class="r">summary(m.L2.DFOP, data = FALSE)
 </code></pre>
 
-<pre><code>## mkin version:    0.9.32 
+<pre><code>## mkin version:    0.9.33 
 ## R version:       3.1.1 
-## Date of fit:     Thu Jul 24 10:32:12 2014 
-## Date of summary: Thu Jul 24 10:32:12 2014 
+## Date of fit:     Mon Aug 25 10:34:18 2014 
+## Date of summary: Mon Aug 25 10:34:18 2014 
 ## 
 ## Equations:
-## [1] d_parent = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) * parent
+## d_parent = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * 
+##            time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * 
+##            time))) * parent
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted with method Marq using 43 model solutions performed in 0.271 s
+## Fitted with method Marq using 43 model solutions performed in 0.26 s
 ## 
 ## Weighting: none
 ## 
@@ -695,17 +703,17 @@ plot(m.L3.SFO)
 <pre><code class="r">summary(m.L3.SFO)
 </code></pre>
 
-<pre><code>## mkin version:    0.9.32 
+<pre><code>## mkin version:    0.9.33 
 ## R version:       3.1.1 
-## Date of fit:     Thu Jul 24 10:32:14 2014 
-## Date of summary: Thu Jul 24 10:32:14 2014 
+## Date of fit:     Mon Aug 25 10:34:18 2014 
+## Date of summary: Mon Aug 25 10:34:18 2014 
 ## 
 ## Equations:
-## [1] d_parent = - k_parent_sink * parent
+## d_parent = - k_parent_sink * parent
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted with method Marq using 44 model solutions performed in 0.251 s
+## Fitted with method Marq using 44 model solutions performed in 0.252 s
 ## 
 ## Weighting: none
 ## 
@@ -781,17 +789,17 @@ plot(m.L3.FOMC)
 <pre><code class="r">summary(m.L3.FOMC, data = FALSE)
 </code></pre>
 
-<pre><code>## mkin version:    0.9.32 
+<pre><code>## mkin version:    0.9.33 
 ## R version:       3.1.1 
-## Date of fit:     Thu Jul 24 10:32:14 2014 
-## Date of summary: Thu Jul 24 10:32:14 2014 
+## Date of fit:     Mon Aug 25 10:34:19 2014 
+## Date of summary: Mon Aug 25 10:34:19 2014 
 ## 
 ## Equations:
-## [1] d_parent = - (alpha/beta) * ((time/beta) + 1)^-1 * parent
+## d_parent = - (alpha/beta) * ((time/beta) + 1)^-1 * parent
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted with method Marq using 26 model solutions performed in 0.154 s
+## Fitted with method Marq using 26 model solutions performed in 0.148 s
 ## 
 ## Weighting: none
 ## 
@@ -854,17 +862,19 @@ plot(m.L3.DFOP)
 <pre><code class="r">summary(m.L3.DFOP, data = FALSE)
 </code></pre>
 
-<pre><code>## mkin version:    0.9.32 
+<pre><code>## mkin version:    0.9.33 
 ## R version:       3.1.1 
-## Date of fit:     Thu Jul 24 10:32:14 2014 
-## Date of summary: Thu Jul 24 10:32:14 2014 
+## Date of fit:     Mon Aug 25 10:34:19 2014 
+## Date of summary: Mon Aug 25 10:34:19 2014 
 ## 
 ## Equations:
-## [1] d_parent = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) * parent
+## d_parent = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * 
+##            time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * 
+##            time))) * parent
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted with method Marq using 37 model solutions performed in 0.228 s
+## Fitted with method Marq using 37 model solutions performed in 0.236 s
 ## 
 ## Weighting: none
 ## 
@@ -950,17 +960,17 @@ plot(m.L4.SFO)
 <pre><code class="r">summary(m.L4.SFO, data = FALSE)
 </code></pre>
 
-<pre><code>## mkin version:    0.9.32 
+<pre><code>## mkin version:    0.9.33 
 ## R version:       3.1.1 
-## Date of fit:     Thu Jul 24 10:32:15 2014 
-## Date of summary: Thu Jul 24 10:32:15 2014 
+## Date of fit:     Mon Aug 25 10:34:19 2014 
+## Date of summary: Mon Aug 25 10:34:19 2014 
 ## 
 ## Equations:
-## [1] d_parent = - k_parent_sink * parent
+## d_parent = - k_parent_sink * parent
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted with method Marq using 20 model solutions performed in 0.141 s
+## Fitted with method Marq using 20 model solutions performed in 0.123 s
 ## 
 ## Weighting: none
 ## 
@@ -1025,17 +1035,17 @@ plot(m.L4.FOMC)
 <pre><code class="r">summary(m.L4.FOMC, data = FALSE)
 </code></pre>
 
-<pre><code>## mkin version:    0.9.32 
+<pre><code>## mkin version:    0.9.33 
 ## R version:       3.1.1 
-## Date of fit:     Thu Jul 24 10:32:15 2014 
-## Date of summary: Thu Jul 24 10:32:15 2014 
+## Date of fit:     Mon Aug 25 10:34:20 2014 
+## Date of summary: Mon Aug 25 10:34:20 2014 
 ## 
 ## Equations:
-## [1] d_parent = - (alpha/beta) * ((time/beta) + 1)^-1 * parent
+## d_parent = - (alpha/beta) * ((time/beta) + 1)^-1 * parent
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted with method Marq using 48 model solutions performed in 0.296 s
+## Fitted with method Marq using 48 model solutions performed in 0.281 s
 ## 
 ## Weighting: none
 ## 
diff --git a/vignettes/FOCUS_Z.pdf b/vignettes/FOCUS_Z.pdf
index 210ce099..ca6d2506 100644
--- a/vignettes/FOCUS_Z.pdf
+++ b/vignettes/FOCUS_Z.pdf