McCall 2,4,5-T dataset, lower bound for transform_rates=FALSE

7 files changed, 52 insertions, 6 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index f348248..345438e 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -17,7 +17,7 @@ Description: Calculation routines based on the FOCUS Kinetics Report (2006).
   implied for correctness of results or fitness for a particular purpose.
   An experimental graphical user interface is included which is based
   on the gWidgetsWWW2 package on GitHub (http://github.com/jverzani/gWidgetsWWW2/)
-Depends: minpack.lm
+Depends: minpack.lm, rootSolve
 Imports: FME, deSolve
 Suggests: knitr, RUnit, gWidgetsWWW2, RSVGTipsDevice
 License: GPL
diff --git a/NAMESPACE b/NAMESPACE
index 4128213..dfa9e9b 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -5,5 +5,6 @@ exportPattern(".")
 import(
   FME,
   deSolve,
-  minpack.lm
+  minpack.lm,
+  rootSolve
 )
diff --git a/R/mkinfit.R b/R/mkinfit.R
index a623146..16b4ff3 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -208,8 +208,17 @@ mkinfit <- function(mkinmod, observed,
     return(mC)

   }

 

+  lower <- rep(-Inf, length(c(state.ini.optim, parms.optim)))

+  names(lower) <- c(names(state.ini.optim), names(parms.optim))

+  if (!transform_rates) {

+    index_k <- grep("^k_", names(lower))

+    lower[index_k] <- 0

+    other_rate_parms <- intersect(c("alpha", "beta", "k1", "k2"), names(lower))

+    lower[other_rate_parms] <- 0

+  }

+

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

-                method = method.modFit, control = control.modFit, ...)

+                method = method.modFit, control = control.modFit, lower = lower, ...)

 

   # Reiterate the fit until convergence of the variance components (IRLS)

   # if requested by the user

@@ -234,7 +243,7 @@ mkinfit <- function(mkinmod, observed,
       sigma.old <- sqrt(fit$var_ms_unweighted)

       observed[err] <- sqrt(fit$var_ms_unweighted)[as.character(observed$name)]

       fit <- modFit(cost, fit$par, method = method.modFit,

-                    control = control.modFit, ...)

+                    control = control.modFit, lower = lower, ...)

       reweight.diff = sum((sqrt(fit$var_ms_unweighted) - sigma.old)^2)

       if (!quiet) {

         cat("Iteration", n.iter, "yields variance estimates:\n")

diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R
index 31200c7..24de074 100644
--- a/R/transform_odeparms.R
+++ b/R/transform_odeparms.R
@@ -44,8 +44,8 @@ transform_odeparms <- function(parms, mod_vars,
     }

   }

 

-  # Transform parameters also for FOMC, DFOP and HS models if requested

-  for (pname in c("alpha", "beta", "k1", "k2", "tb")) {

+  # Transform rate and fraction parameters also for FOMC, DFOP and HS models if requested

+  for (pname in c("alpha", "beta", "k1", "k2")) {

     if (!is.na(parms[pname])) {

       transparms[pname] <- ifelse(transform_rates, log(parms[pname]), parms[pname])

       transparms[pname] <- ifelse(transform_rates, log(parms[pname]), parms[pname])

diff --git a/data/mccall81_245T_soil_metabolism.RData b/data/mccall81_245T_soil_metabolism.RData
new file mode 100644
index 0000000..1e6efa7
--- /dev/null
+++ b/data/mccall81_245T_soil_metabolism.RData
diff --git a/man/mccall81_245T_soil_metabolism.Rd b/man/mccall81_245T_soil_metabolism.Rd
new file mode 100644
index 0000000..2ccc23d
--- /dev/null
+++ b/man/mccall81_245T_soil_metabolism.Rd
@@ -0,0 +1,35 @@
+\name{mccall81_245T_soil_metabolism}

+\Rdversion{1.1}

+\alias{commerce}

+\alias{catlin}

+\alias{fargo}

+\alias{keith}

+\alias{walla_walla}

+\alias{cecil}

+\docType{data}

+\title{

+  Datasets on aerobic soil metabolism of 2,4,5-T in six soils

+}

+\description{

+ Time course of 2,4,5-trichlorophenoxyacetic acid, and the corresponding 2,4,5-trichlorophenol and 

+ 2,4,5-trichloroanisole as recovered in diethylether extracts.

+}

+\usage{mccall81_245T_soil_metabolism}

+\format{

+  6 dataframes containing the following variables.

+  \describe{

+    \item{\code{name}}{the name of the compound observed. Note that T245 is used as an acronym for 2,4,5-T. T245 is 

+                       a legitimate object name in R, which is necessary for

+                       specifying models using \code{\link{mkinmod}}.}

+    \item{\code{time}}{a numeric vector containing sampling times in days after treatment}

+    \item{\code{value}}{a numeric vector containing concentrations in percent of applied radioactivity}

+  }

+}

+\source{

+  McCall P, Vrona SA, Kelley SS (1981) Fate of uniformly carbon-14 ring labeled 2,4,5-Trichlorophenoxyacetic acid and 2,4-dichlorophenoxyacetic acid. J Agric Chem 29, 100-107

+  \url{http://dx.doi.org/10.1021/jf00103a026}

+}

+\examples{

+  commerce

+}

+\keyword{datasets}

diff --git a/man/mkinfit.Rd b/man/mkinfit.Rd
index 839580a..823ceee 100644
--- a/man/mkinfit.Rd
+++ b/man/mkinfit.Rd
@@ -110,6 +110,7 @@ mkinfit(mkinmod, observed,
     assumption of normal distribution of the estimator. If TRUE, also 
     alpha and beta parameters of the FOMC model are log-transformed, as well 
     as k1 and k2 rate constants for the DFOP and HS models.
+    If TRUE, zero is used as a lower bound for the rates in the optimisation.
   }
   \item{transform_fractions}{
     Boolean specifying if formation fractions constants should be transformed in the