aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorjranke <jranke@edb9625f-4e0d-4859-8d74-9fd3b1da38cb>2010-05-20 18:02:42 +0000
committerjranke <jranke@edb9625f-4e0d-4859-8d74-9fd3b1da38cb>2010-05-20 18:02:42 +0000
commita4421eba19eae98a0bd00adb4e8c6d72cc49f9fb (patch)
tree5692b056191b9197c900404410b17306da6526db
parent16f5b1d3c0136413e92b2be0f20d365e92e9cd1c (diff)
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
-rw-r--r--DESCRIPTION29
-rw-r--r--R/mkin_wide_to_long.R2
-rw-r--r--R/mkinfit.R36
-rw-r--r--data/schaefer07_complex_case.RDatabin0 -> 371 bytes
-rw-r--r--man/mkinfit.Rd11
-rw-r--r--man/schaefer07_complex_case.Rd39
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
Binary files differ
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}

Contact - Imprint