Various improvements, the most prominent being the addition of the

test dataset from the original KinGUI Piacenza paper.



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

6 files changed, 91 insertions, 26 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index fb137a76..e518d6df 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,12 +1,17 @@
-Package: mkin

-Type: Package

-Title: Routines for fitting kinetic models with one or more state variables to chemical degradation data

-Version: 0.7-1

-Date: 2010-05-18

-Author: Johannes Ranke

-Maintainer: Johannes Ranke <jranke@harlan.com>

-Description: Calculation routines based on the FOCUS Kinetics Report (2006)

-Depends: FME

-License: GPL

-LazyLoad: yes

-LazyData: yes

+Package: mkin
+Type: Package
+Title: Routines for fitting kinetic models with one or more state
+        variables to chemical degradation data
+Version: 0.7-3
+Date: 2010-05-20
+Author: Johannes Ranke
+Maintainer: Johannes Ranke <jranke@harlan.com>
+Description: Calculation routines based on the FOCUS Kinetics Report
+        (2006)
+Depends: FME
+License: GPL
+LazyLoad: yes
+LazyData: yes
+Packaged: 2010-05-18 12:52:36 UTC; ranke
+Repository: CRAN
+Date/Publication: 2010-05-18 13:18:05
diff --git a/R/mkin_wide_to_long.R b/R/mkin_wide_to_long.R
index 6db2f337..21ab77b9 100644
--- a/R/mkin_wide_to_long.R
+++ b/R/mkin_wide_to_long.R
@@ -1,9 +1,9 @@
 mkin_wide_to_long <- function(wide_data, time = "t")

 {

   colnames <- names(wide_data)

+  if (!(time %in% colnames)) stop("The data in wide format have to contain a variable named ", time, ".")

   vars <- subset(colnames, colnames != time)

   n <- length(colnames) - 1

-  if (!(time %in% colnames)) stop("The data in wide format have to contain a variable named ", time, ".")

   long_data <- data.frame(

     name = rep(vars, each = length(wide_data[[time]])),

     time = rep(wide_data[[time]], n),

diff --git a/R/mkinfit.R b/R/mkinfit.R
index c1490276..9e872fcb 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -1,6 +1,7 @@
 mkinfit <- function(mkinmod, observed,

   parms.ini = rep(0.1, length(mkinmod$parms)),

   state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)), 

+  lower = 0, upper = Inf,

   fixed_parms = NULL,

   fixed_initials = names(mkinmod$diffs)[-1],

   plot = FALSE, quiet = FALSE,

@@ -58,11 +59,7 @@ mkinfit <- function(mkinmod, observed,
 

     odeparms <- c(P[(length(state.ini.optim) + 1):length(P)], parms.fixed)

 

-    # Get more timepoints if plotting is desired

-    if(plot) {

-      outtimes = unique(sort(c(observed$time, 

-        seq(min(observed$time), max(observed$time), length.out=100))))

-    } else outtimes = unique(observed$time)

+    outtimes = unique(observed$time)

 

     # Solve the ode

     out <- ode(

@@ -70,7 +67,7 @@ mkinfit <- function(mkinmod, observed,
       times = outtimes,

       func = mkindiff, 

       parms = odeparms)

-     

+

   

     # Output transformation for models with unobserved compartments like SFORB

     out_transformed <- data.frame(time = out[,"time"])

@@ -81,7 +78,7 @@ mkinfit <- function(mkinmod, observed,
         out_transformed[var] <- rowSums(out[, mkinmod$map[[var]]])

       }

     }    

-    assign("out_predicted", subset(out_transformed, time %in% observed$time), inherits=TRUE)

+    assign("out_predicted", out_transformed, inherits=TRUE)

 

     mC <- modCost(out_transformed, observed, y = "value",

       err = err, weight = weight, scaleVar = scaleVar)

@@ -90,8 +87,23 @@ mkinfit <- function(mkinmod, observed,
     if (mC$model < cost.old) {

       if(!quiet) cat("Model cost at call ", calls, ": ", mC$model, "\n")

 

-      # Plot the data and the current model output if requested

+      # Plot the data and current model output if requested

       if(plot) {

+        outtimes_plot = seq(min(observed$time), max(observed$time), length.out=100)

+        out_plot <- ode(

+          y = odeini,

+          times = outtimes_plot,

+          func = mkindiff, 

+          parms = odeparms)

+        out_transformed_plot <- data.frame(time = out_plot[,"time"])

+        for (var in names(mkinmod$map)) {

+          if(length(mkinmod$map[[var]]) == 1) {

+            out_transformed_plot[var] <- out_plot[, var]

+          } else {

+            out_transformed_plot[var] <- rowSums(out_plot[, mkinmod$map[[var]]])

+          }

+        }    

+

         plot(0, type="n", 

           xlim = range(observed$time), ylim = range(observed$value, na.rm=TRUE),

           xlab = "Time", ylab = "Observed")

@@ -101,7 +113,7 @@ mkinfit <- function(mkinmod, observed,
           points(subset(observed, name == obs_var, c(time, value)), 

             pch = pch_obs[obs_var], col = col_obs[obs_var])

         }

-        matlines(out_transformed$time, out_transformed[-1])

+        matlines(out_transformed_plot$time, out_transformed_plot[-1])

         legend("topright", inset=c(0.05, 0.05), legend=obs_vars, 

           col=col_obs, pch=pch_obs, lty=1:length(pch_obs))

       }

@@ -110,7 +122,7 @@ mkinfit <- function(mkinmod, observed,
     }

     return(mC)

   }

-  fit <- modFit(cost, c(state.ini.optim, parms.optim), ...)

+  fit <- modFit(cost, c(state.ini.optim, parms.optim), lower = lower, upper = upper, ...)

 

   # We need the function for plotting

   fit$mkindiff <- mkindiff

@@ -125,8 +137,8 @@ mkinfit <- function(mkinmod, observed,
   # Collect initial parameter values in two dataframes

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

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

-  if (exists("lower")) fit$start$lower <- lower

-  if (exists("upper")) fit$start$upper <- upper

+  fit$start$lower <- lower

+  fit$start$upper <- upper

 

   fit$fixed <- data.frame(

     value = c(state.ini.fixed, parms.fixed))

diff --git a/data/schaefer07_complex_case.RData b/data/schaefer07_complex_case.RData
new file mode 100644
index 00000000..0e3b3af2
--- /dev/null
+++ b/data/schaefer07_complex_case.RData
diff --git a/man/mkinfit.Rd b/man/mkinfit.Rd
index 640207ce..0e381afd 100644
--- a/man/mkinfit.Rd
+++ b/man/mkinfit.Rd
@@ -10,7 +10,7 @@
   values.
 }
 \usage{
-mkinfit(mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$parms)), state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)), fixed_parms = NULL, fixed_initials = names(mkinmod$diffs)[-1], plot = FALSE, quiet = FALSE, err = NULL, weight = "none", scaleVar = FALSE, ...)
+mkinfit(mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$parms)), state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)), lower = 0, upper = Inf, fixed_parms = NULL, fixed_initials = names(mkinmod$diffs)[-1], plot = FALSE, quiet = FALSE, err = NULL, weight = "none", scaleVar = FALSE, ...)
 }
 \arguments{
   \item{mkinmod}{
@@ -37,6 +37,15 @@ mkinfit(mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$parms)), state.in
     \code{\link{mkinmod}}). The default is to set the initial value of the first model
     variable to 100 and all others to 0.
   }
+  \item{lower}{
+    Lower bounds for the parameters, passed to \code{\link{modFit}}. Defaults to 0 because
+    negative values to not make sense for the models currently created by \code{\link{mkinmod}}.
+  }
+  \item{upper}{
+    Upper bounds for the parameters, passed to \code{\link{modFit}}. Defaults to \code{Inf}. Setting
+    non-infinite upper bounds has a strong impact on performance, and using a method like "L-BFGS-B" (by 
+    specifying an additional \code{method} argument) is recommended.
+  }
   \item{fixed_parms}{
     The names of parameters that should not be optimised but rather kept at the values
     specified in \code{parms.ini}.
diff --git a/man/schaefer07_complex_case.Rd b/man/schaefer07_complex_case.Rd
new file mode 100644
index 00000000..2978e25d
--- /dev/null
+++ b/man/schaefer07_complex_case.Rd
@@ -0,0 +1,39 @@
+\name{schaefer07_complex_case}

+\alias{schaefer07_complex_case}

+\encoding{latin1}

+\docType{data}

+\title{

+  Metabolism data set used for checking the software quality of KinGUI

+}

+\description{

+  This dataset was used for a comparison of KinGUI and ModelMaker to check the

+  software quality of KinGUI in the original publication (Schäfer et al., 2007).

+}

+\usage{data(schaefer07_complex_case)}

+\format{

+  A data frame with 8 observations on the following 6 variables.

+  \describe{

+    \item{\code{time}}{a numeric vector}

+    \item{\code{parent}}{a numeric vector}

+    \item{\code{A1}}{a numeric vector}

+    \item{\code{B1}}{a numeric vector}

+    \item{\code{C1}}{a numeric vector}

+    \item{\code{A2}}{a numeric vector}

+  }

+}

+\source{

+  Schäfer D, Mikolasch M, Rainbird P and Harvey B (2007). KinGUI: a new kinetic

+  software tool for evaluations according to FOCUS degradation kinetics. In: Del

+  Re AAM, Capri E, Fragoulis G and Trevisan M (Eds.). Proceedings of the XIII

+  Symposium Pesticide Chemistry, Piacenza, 2007, p. 916-923.  }

+\examples{

+data <- mkin_wide_to_long(schaefer07_complex_case, time = "time")

+model <- mkinmod(

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

+\dontrun{mkinfit(model, data, plot=TRUE)}

+}

+\keyword{datasets}