From a4421eba19eae98a0bd00adb4e8c6d72cc49f9fb Mon Sep 17 00:00:00 2001 From: jranke Date: Thu, 20 May 2010 18:02:42 +0000 Subject: 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 --- DESCRIPTION | 29 +++++++++++++++------------ R/mkin_wide_to_long.R | 2 +- R/mkinfit.R | 36 ++++++++++++++++++++++------------ data/schaefer07_complex_case.RData | Bin 0 -> 371 bytes man/mkinfit.Rd | 11 ++++++++++- man/schaefer07_complex_case.Rd | 39 +++++++++++++++++++++++++++++++++++++ 6 files changed, 91 insertions(+), 26 deletions(-) create mode 100644 data/schaefer07_complex_case.RData create mode 100644 man/schaefer07_complex_case.Rd 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 -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 +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 Binary files /dev/null and b/data/schaefer07_complex_case.RData 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} -- cgit v1.2.1