diff options
| author | Johannes Ranke <jranke@uni-bremen.de> | 2014-05-07 14:47:28 +0200 |
|---|---|---|
| committer | Johannes Ranke <jranke@uni-bremen.de> | 2014-05-07 14:47:28 +0200 |
| commit | e959fde98f95f3595e01490b67892678bbcd1b27 (patch) | |
| tree | 992c56223a31c6937091dd5f9eeef63c2dd9e579 | |
| parent | d846ac7691ab648afbb5a98bbca91911396a95bf (diff) | |
Fork the gmkin GUI from mkin. See ChangeLog for details
98 files changed, 43 insertions, 5292 deletions
@@ -1,3 +1,11 @@ +2014-05-07 Johannes Ranke <jranke@uni-bremen.de> for gmkin (0.5-01)
+
+ * Remove all content that does not belong to the newly established
+ GUI package gmkin
+ * Update DESCRIPTION, README.md, NAMESPACE, GNUmakefile
+ * Update R/gmkin.R and inst/GUI/gmkin.R
+ * Move inst/GUI/TODO to TODO
+
2014-04-28 Johannes Ranke <jranke@uni-bremen.de> for mkin (0.9-26)
* DESCRIPTION: Added copyright information and location of gWidgetsWWW2
diff --git a/DESCRIPTION b/DESCRIPTION index 594d68c..2631194 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,29 +1,18 @@ -Package: mkin +Package: gmkin Type: Package -Title: Routines for fitting kinetic models with one or more state - variables to chemical degradation data -Version: 0.9-26 -Date: 2014-05-02 +Title: GUI for fitting kinetic models to chemical degradation data with mkin +Version: 0.5-1 +Date: 2014-05-07 Authors@R: c(person("Johannes", "Ranke", role = c("aut", "cre", "cph"), email = "jranke@uni-bremen.de"), - person("Katrin", "Lindenberger", role = "ctb"), - person("René", "Lehmann", role = "ctb"), person("Eurofins Regulatory AG", role = "cph")) -Description: Calculation routines based on the FOCUS Kinetics Report (2006). - Includes a function for conveniently defining differential equation models, - model solution based on eigenvalues if possible or using numerical solvers - and a choice of the optimisation methods made available by the FME package - (default is a Levenberg-Marquardt variant). Please note that no warranty is - 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, rootSolve -Imports: FME, deSolve -Suggests: knitr, RUnit, gWidgetsWWW2 +Description: This package contains a browser based graphical user interface + for R package mkin, based on the gWidgetsWWW2 package on GitHub + (http://github.com/jverzani/gWidgetsWWW2/) +Depends: mkin (>= 0.9-26), gWidgetsWWW2 License: GPL LazyLoad: yes LazyData: yes Encoding: UTF-8 -VignetteBuilder: knitr -BugReports: http://github.com/jranke/mkin/issues -URL: http://cran.r-project.org/package=mkin, http://kinfit.r-forge.r-project.org, http://github.com/jranke/mkin +BugReports: http://github.com/jranke/gmkin/issues +URL: http://github.com/jranke/gmkin diff --git a/GNUmakefile b/GNUmakefile index 4511b07..7466d85 100644 --- a/GNUmakefile +++ b/GNUmakefile @@ -17,25 +17,10 @@ help: @echo "Development Tasks" @echo "-----------------" @echo " build Create the package" - @echo " build-no-vignettes Create the package without rebuilding vignettes" @echo " check Invoke build and then check the package" - @echo " check-no-vignettes Invoke build without rebuilding vignettes, and then check" @echo " install Invoke build and then install the result" - @echo " install-no-vignettes Invoke build without rebuilding vignettes and then install the result" - @echo " test Install a new copy of the package and run it " - @echo " through the testsuite" - @echo " test-no-vignettes Invoke build without rebuilding vignettes, and then run it" - @echo " through the testsuite" - @echo "" - @echo "Packaging Tasks" - @echo "---------------" - @echo " release Give some reminders" - @echo "" - @echo "Using R in: $(RBIN)" - @echo "Set the RBIN environment variable to change this." @echo "" - #------------------------------------------------------------------------------ # Development Tasks #------------------------------------------------------------------------------ @@ -44,45 +29,10 @@ build: cd ..;\ "$(RBIN)/R" CMD build $(PKGSRC) -build-no-vignettes: - cd ..;\ - "$(RBIN)/R" CMD build $(PKGSRC) --no-build-vignettes - install: build cd ..;\ "$(RBIN)/R" CMD INSTALL $(PKGNAME)_$(PKGVERS).tar.gz -install-no-vignettes: build-no-vignettes - cd ..;\ - "$(RBIN)/R" CMD INSTALL $(PKGNAME)_$(PKGVERS).tar.gz - check: build cd ..;\ - "$(RBIN)/R" CMD check --as-cran --no-tests $(PKGNAME)_$(PKGVERS).tar.gz - -check-no-vignettes: build-no-vignettes - cd ..;\ - "$(RBIN)/R" CMD check --as-cran --no-tests $(PKGNAME)_$(PKGVERS).tar.gz - -test: install - cd tests;\ - "$(RBIN)/Rscript" doRUnit.R - -test-no-vignettes: install-no-vignettes - cd tests;\ - "$(RBIN)/Rscript" doRUnit.R - -#------------------------------------------------------------------------------ -# Packaging Tasks -#------------------------------------------------------------------------------ -release: - @echo "\nHow about make test and make check?" - @echo "\nIs the DESCRIPTION file up to date?" - @echo "\nTo update the svn repository tied to the local r-forge branch with" - @echo "changes in the local master branch, run:" - @echo "'git checkout r-forge'" - @echo "'git merge --squash master'" - @echo "'git commit'" - @echo "'git svn dcommit'" - @echo "\nThen change back to the master branch:" - @echo "'git checkout master'" + "$(RBIN)/R" CMD check --as-cran $(PKGNAME)_$(PKGVERS).tar.gz @@ -2,9 +2,4 @@ exportPattern(".") # Import all packages listed as Imports or Depends -import( - FME, - deSolve, - minpack.lm, - rootSolve -) +import(mkin, gWidgetsWWW2) diff --git a/R/DFOP.solution.R b/R/DFOP.solution.R deleted file mode 100644 index 8531cfe..0000000 --- a/R/DFOP.solution.R +++ /dev/null @@ -1,5 +0,0 @@ -DFOP.solution <- function(t, parent.0, k1, k2, g)
-{
- parent = g * parent.0 * exp(-k1 * t) +
- (1 - g) * parent.0 * exp(-k2 * t)
-}
diff --git a/R/FOMC.solution.R b/R/FOMC.solution.R deleted file mode 100644 index 8bd13d6..0000000 --- a/R/FOMC.solution.R +++ /dev/null @@ -1,4 +0,0 @@ -FOMC.solution <- function(t, parent.0, alpha, beta)
-{
- parent = parent.0 / (t/beta + 1)^alpha
-}
diff --git a/R/HS.solution.R b/R/HS.solution.R deleted file mode 100644 index 4651a6a..0000000 --- a/R/HS.solution.R +++ /dev/null @@ -1,6 +0,0 @@ -HS.solution <- function(t, parent.0, k1, k2, tb)
-{
- parent = ifelse(t <= tb,
- parent.0 * exp(-k1 * t),
- parent.0 * exp(-k1 * tb) * exp(-k2 * (t - tb)))
-}
diff --git a/R/SFO.solution.R b/R/SFO.solution.R deleted file mode 100644 index 3a376e4..0000000 --- a/R/SFO.solution.R +++ /dev/null @@ -1,4 +0,0 @@ -SFO.solution <- function(t, parent.0, k)
-{
- parent = parent.0 * exp(-k * t)
-}
diff --git a/R/SFORB.solution.R b/R/SFORB.solution.R deleted file mode 100644 index 45a4533..0000000 --- a/R/SFORB.solution.R +++ /dev/null @@ -1,9 +0,0 @@ -SFORB.solution = function(t, parent.0, k_12, k_21, k_1output) {
- sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 + k_12 * k_21 - (k_12 + k_1output) * k_21)
- b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp
- b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp
-
- parent = parent.0 *
- (((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) +
- ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t))
-}
diff --git a/R/endpoints.R b/R/endpoints.R deleted file mode 100644 index c9a6a51..0000000 --- a/R/endpoints.R +++ /dev/null @@ -1,108 +0,0 @@ -endpoints <- function(fit, pseudoDT50 = FALSE) {
- # Calculate dissipation times DT50 and DT90 and, if necessary, formation
- # fractions and SFORB eigenvalues from optimised parameters
- ep <- list()
- obs_vars <- fit$obs_vars
- parms.all <- fit$bparms.ode
- ep$distimes <- data.frame(DT50 = rep(NA, length(obs_vars)),
- DT90 = rep(NA, length(obs_vars)),
- row.names = obs_vars)
- ep$ff <- vector()
- ep$SFORB <- vector()
- for (obs_var in obs_vars) {
- type = names(fit$mkinmod$map[[obs_var]])[1]
-
- # Get formation fractions if directly fitted, and calculate remaining fraction to sink
- f_names = grep(paste("f", obs_var, sep = "_"), names(parms.all), value=TRUE)
- f_values = parms.all[f_names]
- f_to_sink = 1 - sum(f_values)
- names(f_to_sink) = ifelse(type == "SFORB",
- paste(obs_var, "free", "sink", sep = "_"),
- paste(obs_var, "sink", sep = "_"))
- for (f_name in f_names) {
- ep$ff[[sub("f_", "", sub("_to_", "_", f_name))]] = f_values[[f_name]]
- }
- ep$ff = append(ep$ff, f_to_sink)
-
- # Get the rest
- if (type == "SFO") {
- k_names = grep(paste("k", obs_var, sep="_"), names(parms.all), value=TRUE)
- k_tot = sum(parms.all[k_names])
- DT50 = log(2)/k_tot
- DT90 = log(10)/k_tot
- for (k_name in k_names)
- {
- ep$ff[[sub("k_", "", k_name)]] = parms.all[[k_name]] / k_tot
- }
- }
- if (type == "FOMC") {
- alpha = parms.all["alpha"]
- beta = parms.all["beta"]
- DT50 = beta * (2^(1/alpha) - 1)
- DT90 = beta * (10^(1/alpha) - 1)
- }
- if (type == "DFOP") {
- k1 = parms.all["k1"]
- k2 = parms.all["k2"]
- g = parms.all["g"]
- f <- function(t, x) {
- fraction <- g * exp( - k1 * t) + (1 - g) * exp( - k2 * t)
- (fraction - (1 - x/100))^2
- }
- DTmax <- 1000
- DT50.o <- optimize(f, c(0.001, DTmax), x=50)$minimum
- DT50 = ifelse(DTmax - DT50.o < 0.1, NA, DT50.o)
- DT90.o <- optimize(f, c(0.001, DTmax), x=90)$minimum
- DT90 = ifelse(DTmax - DT90.o < 0.1, NA, DT90.o)
- }
- if (type == "HS") {
- k1 = parms.all["k1"]
- k2 = parms.all["k2"]
- tb = parms.all["tb"]
- DTx <- function(x) {
- DTx.a <- (log(100/(100 - x)))/k1
- DTx.b <- tb + (log(100/(100 - x)) - k1 * tb)/k2
- if (DTx.a < tb) DTx <- DTx.a
- else DTx <- DTx.b
- return(DTx)
- }
- DT50 <- DTx(50)
- DT90 <- DTx(90)
- }
- if (type == "SFORB") {
- # FOCUS kinetics (2006), p. 60 f
- k_out_names = grep(paste("k", obs_var, "free", sep="_"), names(parms.all), value=TRUE)
- k_out_names = setdiff(k_out_names, paste("k", obs_var, "free", "bound", sep="_"))
- k_1output = sum(parms.all[k_out_names])
- k_12 = parms.all[paste("k", obs_var, "free", "bound", sep="_")]
- k_21 = parms.all[paste("k", obs_var, "bound", "free", sep="_")]
-
- sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 + k_12 * k_21 - (k_12 + k_1output) * k_21)
- b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp
- b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp
-
- SFORB_fraction = function(t) {
- ((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) +
- ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t)
- }
- f_50 <- function(t) (SFORB_fraction(t) - 0.5)^2
- max_DT <- 1000
- DT50.o <- optimize(f_50, c(0.01, max_DT))$minimum
- if (abs(DT50.o - max_DT) < 0.01) DT50 = NA else DT50 = DT50.o
- f_90 <- function(t) (SFORB_fraction(t) - 0.1)^2
- DT90.o <- optimize(f_90, c(0.01, max_DT))$minimum
- if (abs(DT90.o - max_DT) < 0.01) DT90 = NA else DT90 = DT90.o
-
- for (k_out_name in k_out_names)
- {
- ep$ff[[sub("k_", "", k_out_name)]] = parms.all[[k_out_name]] / k_1output
- }
-
- # Return the eigenvalues for comparison with DFOP rate constants
- ep$SFORB[[paste(obs_var, "b1", sep="_")]] = b1
- ep$SFORB[[paste(obs_var, "b2", sep="_")]] = b2
- }
- ep$distimes[obs_var, ] = c(DT50, DT90)
- }
- return(ep)
-}
@@ -1,11 +1,3 @@ gmkin <- function() { - if (require(gWidgetsWWW2)) load_app(system.file("GUI/gmkin.R", package = "mkin")) - else { - message( - "\nYou need to install gWidgetsWWW2 in order to run the mkin GUI gmkin.\n", - "This package is not currently on CRAN but can be installed from github\n", - "using the devtools package:\n", - "library(devtools)\n", - "install_github('gWidgetsWWW2', 'jverzani')") - } + load_app(system.file("GUI/gmkin.R", package = "gmkin")) } diff --git a/R/ilr.R b/R/ilr.R deleted file mode 100644 index 389653e..0000000 --- a/R/ilr.R +++ /dev/null @@ -1,51 +0,0 @@ -# $Id$ - -# Copyright (C) 2012 René Lehmann, Johannes Ranke -# Contact: jranke@uni-bremen.de - -# This file is part of the R package mkin - -# mkin is free software: you can redistribute it and/or modify it under the -# terms of the GNU General Public License as published by the Free Software -# Foundation, either version 3 of the License, or (at your option) any later -# version. - -# This program is distributed in the hope that it will be useful, but WITHOUT -# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more -# details. - -# You should have received a copy of the GNU General Public License along with -# this program. If not, see <http://www.gnu.org/licenses/> - -ilr <- function(x) { - z <- vector() - for (i in 1:(length(x) - 1)) { - z[i] <- sqrt(i/(i+1)) * log((prod(x[1:i]))^(1/i) / x[i+1]) - } - return(z) -} - -invilr<-function(x) { - D <- length(x) + 1 - z <- c(x, 0) - y <- rep(0, D) - s <- sqrt(1:D*2:(D+1)) - q <- z/s - y[1] <- sum(q[1:D]) - for (i in 2:D) { - y[i] <- sum(q[i:D]) - sqrt((i-1)/i) * z[i-1] - } - z <- vector() - for (i in 1:D) { - z[i] <- exp(y[i])/sum(exp(y)) - } - - # Work around a numerical problem with NaN values returned by the above - # Only works if there is only one NaN value: replace it with 1 - # if the sum of the other components is < 1e-10 - if (sum(is.na(z)) == 1 && sum(z[!is.na(z)]) < 1e-10) - z = ifelse(is.na(z), 1, z) - - return(z) -} diff --git a/R/mkin_long_to_wide.R b/R/mkin_long_to_wide.R deleted file mode 100644 index 87e1afd..0000000 --- a/R/mkin_long_to_wide.R +++ /dev/null @@ -1,30 +0,0 @@ -# $Id$
-
-# Copyright (C) 2010-2011 Johannes Ranke
-# Contact: mkin-devel@lists.berlios.de
-
-# This file is part of the R package mkin
-
-# mkin is free software: you can redistribute it and/or modify it under the
-# terms of the GNU General Public License as published by the Free Software
-# Foundation, either version 3 of the License, or (at your option) any later
-# version.
-
-# This program is distributed in the hope that it will be useful, but WITHOUT
-# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
-# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-# details.
-
-# You should have received a copy of the GNU General Public License along with
-# this program. If not, see <http://www.gnu.org/licenses/>
-
-mkin_long_to_wide <- function(long_data, time = "time", outtime = "time")
-{
- colnames <- unique(long_data$name)
- wide_data <- data.frame(time = subset(long_data, name == colnames[1], time))
- names(wide_data) <- outtime
- for (var in colnames) {
- wide_data[var] <- subset(long_data, name == var, value)
- }
- return(wide_data)
-}
diff --git a/R/mkin_wide_to_long.R b/R/mkin_wide_to_long.R deleted file mode 100644 index f1814fc..0000000 --- a/R/mkin_wide_to_long.R +++ /dev/null @@ -1,34 +0,0 @@ -# $Id$
-
-# Copyright (C) 2010-2013 Johannes Ranke
-# Contact: mkin-devel@lists.berlios.de
-
-# This file is part of the R package mkin
-
-# mkin is free software: you can redistribute it and/or modify it under the
-# terms of the GNU General Public License as published by the Free Software
-# Foundation, either version 3 of the License, or (at your option) any later
-# version.
-
-# This program is distributed in the hope that it will be useful, but WITHOUT
-# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
-# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-# details.
-
-# You should have received a copy of the GNU General Public License along with
-# this program. If not, see <http://www.gnu.org/licenses/>
-if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "value"))
-
-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
- long_data <- data.frame(
- name = rep(vars, each = length(wide_data[[time]])),
- time = as.numeric(rep(wide_data[[time]], n)),
- value = as.numeric(unlist(wide_data[vars])),
- row.names = NULL)
- return(long_data)
-}
diff --git a/R/mkinerrmin.R b/R/mkinerrmin.R deleted file mode 100644 index 074cc56..0000000 --- a/R/mkinerrmin.R +++ /dev/null @@ -1,85 +0,0 @@ -# $Id$
-
-# Copyright (C) 2010-2013 Johannes Ranke
-# Contact: jranke@uni-bremen.de
-
-# This file is part of the R package mkin
-
-# mkin is free software: you can redistribute it and/or modify it under the
-# terms of the GNU General Public License as published by the Free Software
-# Foundation, either version 3 of the License, or (at your option) any later
-# version.
-
-# This program is distributed in the hope that it will be useful, but WITHOUT
-# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
-# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-# details.
-
-# You should have received a copy of the GNU General Public License along with
-# this program. If not, see <http://www.gnu.org/licenses/>
-if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "value_mean"))
-
-mkinerrmin <- function(fit, alpha = 0.05)
-{
- parms.optim <- fit$par
-
- kinerrmin <- function(errdata, n.parms) {
- means.mean <- mean(errdata$value_mean, na.rm = TRUE)
- df = length(errdata$value_mean) - n.parms
-
- err.min <- sqrt((1 / qchisq(1 - alpha, df)) *
- sum((errdata$value_mean - errdata$value_pred)^2)/(means.mean^2))
-
- return(list(err.min = err.min, n.optim = n.parms, df = df))
- }
-
- means <- aggregate(value ~ time + name, data = fit$observed, mean, na.rm=TRUE)
- errdata <- merge(means, fit$predicted, by = c("time", "name"),
- 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))
-
- n.optim.overall <- length(parms.optim)
-
- errmin.overall <- kinerrmin(errdata, n.optim.overall)
- errmin <- data.frame(err.min = errmin.overall$err.min,
- n.optim = errmin.overall$n.optim, df = errmin.overall$df)
- rownames(errmin) <- "All data"
-
- for (obs_var in fit$obs_vars)
- {
- errdata.var <- subset(errdata, name == obs_var)
-
- # Check if initial value is optimised
- n.initials.optim <- length(grep(paste(obs_var, ".*", "_0", sep=""), names(parms.optim)))
-
- # Rate constants are attributed to the source variable
- n.k.optim <- length(grep(paste("^k", obs_var, sep="_"), names(parms.optim)))
-
- # Formation fractions are attributed to the target variable
- n.ff.optim <- length(grep(paste("^f", ".*", obs_var, "$", sep=""), names(parms.optim)))
-
- n.optim <- n.k.optim + n.initials.optim + n.ff.optim
-
- # FOMC, DFOP and HS parameters are only counted if we are looking at the
- # first variable in the model which is always the source variable
- if (obs_var == fit$obs_vars[[1]]) {
- if ("alpha" %in% names(parms.optim)) n.optim <- n.optim + 1
- if ("beta" %in% names(parms.optim)) n.optim <- n.optim + 1
- if ("k1" %in% names(parms.optim)) n.optim <- n.optim + 1
- if ("k2" %in% names(parms.optim)) n.optim <- n.optim + 1
- if ("g" %in% names(parms.optim)) n.optim <- n.optim + 1
- if ("tb" %in% names(parms.optim)) n.optim <- n.optim + 1
- }
-
- # Calculate and add a line to the results
- errmin.tmp <- kinerrmin(errdata.var, n.optim)
- errmin[obs_var, c("err.min", "n.optim", "df")] <- errmin.tmp
- }
-
- return(errmin)
-}
diff --git a/R/mkinfit.R b/R/mkinfit.R deleted file mode 100644 index f44c67d..0000000 --- a/R/mkinfit.R +++ /dev/null @@ -1,482 +0,0 @@ -# Copyright (C) 2010-2014 Johannes Ranke
-# Portions of this code are copyright (C) 2013 Eurofins Regulatory AG
-# Contact: jranke@uni-bremen.de
-# The summary function is an adapted and extended version of summary.modFit
-# from the FME package, v 1.1 by Soetart and Petzoldt, which was in turn
-# inspired by summary.nls.lm
-
-# This file is part of the R package mkin
-
-# mkin is free software: you can redistribute it and/or modify it under the
-# terms of the GNU General Public License as published by the Free Software
-# Foundation, either version 3 of the License, or (at your option) any later
-# version.
-
-# This program is distributed in the hope that it will be useful, but WITHOUT
-# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
-# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-# details.
-
-# You should have received a copy of the GNU General Public License along with
-# this program. If not, see <http://www.gnu.org/licenses/>
-if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "value"))
-
-mkinfit <- function(mkinmod, observed,
- parms.ini = "auto",
- state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)),
- fixed_parms = NULL,
- fixed_initials = names(mkinmod$diffs)[-1],
- solution_type = "auto",
- method.ode = "lsoda",
- method.modFit = "Marq",
- control.modFit = list(),
- transform_rates = TRUE,
- transform_fractions = TRUE,
- plot = FALSE, quiet = FALSE,
- err = NULL, weight = "none", scaleVar = FALSE,
- atol = 1e-8, rtol = 1e-10, n.outtimes = 100,
- reweight.method = NULL,
- reweight.tol = 1e-8, reweight.max.iter = 10,
- trace_parms = FALSE,
- ...)
-{
- # Get the names of the state variables in the model
- mod_vars <- names(mkinmod$diffs)
-
- # Get the names of observed variables
- obs_vars = names(mkinmod$spec)
-
- # Subset observed data with names of observed data in the model
- observed <- subset(observed, name %in% obs_vars)
-
- # Define starting values for parameters where not specified by the user
- if (parms.ini[[1]] == "auto") parms.ini = vector()
-
- # Prevent inital parameter specifications that are not in the model
- wrongpar.names <- setdiff(names(parms.ini), mkinmod$parms)
- if (length(wrongpar.names) > 0) {
- stop("Initial parameter(s) ", paste(wrongpar.names, collapse = ", "),
- " not used in the model")
- }
-
- k_salt = 0
- defaultpar.names <- setdiff(mkinmod$parms, names(parms.ini))
- for (parmname in defaultpar.names) {
- # Default values for rate constants, depending on the parameterisation
- if (substr(parmname, 1, 2) == "k_") {
- parms.ini[parmname] = 0.1 + k_salt
- k_salt = k_salt + 1e-4
- }
- # Default values for rate constants for reversible binding
- if (grepl("free_bound$", parmname)) parms.ini[parmname] = 0.1
- if (grepl("bound_free$", parmname)) parms.ini[parmname] = 0.02
- # Default values for formation fractions
- if (substr(parmname, 1, 2) == "f_") parms.ini[parmname] = 0.2
- # Default values for the FOMC, DFOP and HS models
- if (parmname == "alpha") parms.ini[parmname] = 1
- if (parmname == "beta") parms.ini[parmname] = 10
- if (parmname == "k1") parms.ini[parmname] = 0.1
- if (parmname == "k2") parms.ini[parmname] = 0.01
- if (parmname == "tb") parms.ini[parmname] = 5
- if (parmname == "g") parms.ini[parmname] = 0.5
- }
-
- # Name the inital state variable values if they are not named yet
- if(is.null(names(state.ini))) names(state.ini) <- mod_vars
-
- # Transform initial parameter values for fitting
- transparms.ini <- transform_odeparms(parms.ini, mod_vars,
- transform_rates = transform_rates,
- transform_fractions = transform_fractions)
-
- # Parameters to be optimised:
- # Kinetic parameters in parms.ini whose names are not in fixed_parms
- parms.fixed <- transparms.ini[fixed_parms]
- parms.optim <- transparms.ini[setdiff(names(transparms.ini), fixed_parms)]
-
- # Inital state variables in state.ini whose names are not in fixed_initials
- state.ini.fixed <- state.ini[fixed_initials]
- state.ini.optim <- state.ini[setdiff(names(state.ini), fixed_initials)]
-
- # Preserve names of state variables before renaming initial state variable
- # parameters
- state.ini.optim.boxnames <- names(state.ini.optim)
- state.ini.fixed.boxnames <- names(state.ini.fixed)
- if(length(state.ini.optim) > 0) {
- names(state.ini.optim) <- paste(names(state.ini.optim), "0", sep="_")
- }
- if(length(state.ini.fixed) > 0) {
- names(state.ini.fixed) <- paste(names(state.ini.fixed), "0", sep="_")
- }
-
- # Decide if the solution of the model can be based on a simple analytical
- # formula, the spectral decomposition of the matrix (fundamental system)
- # or a numeric ode solver from the deSolve package
- if (!solution_type %in% c("auto", "analytical", "eigen", "deSolve"))
- stop("solution_type must be auto, analytical, eigen or de Solve")
- if (solution_type == "analytical" && length(mkinmod$spec) > 1)
- stop("Analytical solution not implemented for models with metabolites.")
- if (solution_type == "eigen" && !is.matrix(mkinmod$coefmat))
- stop("Eigenvalue based solution not possible, coefficient matrix not present.")
- if (solution_type == "auto") {
- if (length(mkinmod$spec) == 1) {
- solution_type = "analytical"
- } else {
- if (is.matrix(mkinmod$coefmat)) {
- solution_type = "eigen"
- if (max(observed$value, na.rm = TRUE) < 0.1) {
- stop("The combination of small observed values (all < 0.1) and solution_type = eigen is error-prone")
- }
- } else {
- solution_type = "deSolve"
- }
- }
- }
-
- cost.old <- 1e100 # The first model cost should be smaller than this value
- calls <- 0 # Counter for number of model solutions
- out_predicted <- NA
- # Define the model cost function
- cost <- function(P)
- {
- assign("calls", calls+1, inherits=TRUE) # Increase the model solution counter
-
- # Trace parameter values if requested
- if(trace_parms) cat(P, "\n")
-
- # Time points at which observed data are available
- # Make sure we include time 0, so initial values for state variables are for time 0
- outtimes = sort(unique(c(observed$time, seq(min(observed$time),
- max(observed$time),
- length.out = n.outtimes))))
-
- if(length(state.ini.optim) > 0) {
- odeini <- c(P[1:length(state.ini.optim)], state.ini.fixed)
- names(odeini) <- c(state.ini.optim.boxnames, state.ini.fixed.boxnames)
- } else {
- odeini <- state.ini.fixed
- names(odeini) <- state.ini.fixed.boxnames
- }
-
- odeparms <- c(P[(length(state.ini.optim) + 1):length(P)], parms.fixed)
-
- parms <- backtransform_odeparms(odeparms, mod_vars,
- transform_rates = transform_rates,
- transform_fractions = transform_fractions)
-
- # Solve the system with current transformed parameter values
- out <- mkinpredict(mkinmod, parms, odeini, outtimes,
- solution_type = solution_type,
- method.ode = method.ode,
- atol = atol, rtol = rtol, ...)
-
- assign("out_predicted", out, inherits=TRUE)
-
- mC <- modCost(out, observed, y = "value",
- err = err, weight = weight, scaleVar = scaleVar)
-
- # Report and/or plot if the model is improved
- if (mC$model < cost.old) {
- if(!quiet) cat("Model cost at call ", calls, ": ", mC$model, "\n")
-
- # 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 <- mkinpredict(mkinmod, parms, odeini, outtimes_plot,
- solution_type = solution_type,
- method.ode = method.ode,
- atol = atol, rtol = rtol, ...)
-
- plot(0, type="n",
- xlim = range(observed$time), ylim = range(observed$value, na.rm=TRUE),
- xlab = "Time", ylab = "Observed")
- col_obs <- pch_obs <- 1:length(obs_vars)
- lty_obs <- rep(1, length(obs_vars))
- names(col_obs) <- names(pch_obs) <- names(lty_obs) <- obs_vars
- for (obs_var in obs_vars) {
- points(subset(observed, name == obs_var, c(time, value)),
- pch = pch_obs[obs_var], col = col_obs[obs_var])
- }
- matlines(out_plot$time, out_plot[-1], col = col_obs, lty = lty_obs)
- legend("topright", inset=c(0.05, 0.05), legend=obs_vars,
- col=col_obs, pch=pch_obs, lty=1:length(pch_obs))
- }
-
- assign("cost.old", mC$model, inherits=TRUE)
- }
- 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, lower = lower, ...)
-
- # Reiterate the fit until convergence of the variance components (IRLS)
- # if requested by the user
- weight.ini = weight
- if (!is.null(err)) weight.ini = "manual"
-
- if (!is.null(reweight.method)) {
- if (reweight.method != "obs") stop("Only reweighting method 'obs' is implemented")
- if(!quiet) {
- cat("IRLS based on variance estimates for each observed variable\n")
- }
- if (!quiet) {
- cat("Initial variance estimates are:\n")
- print(signif(fit$var_ms_unweighted, 8))
- }
- reweight.diff = 1
- n.iter <- 0
- if (!is.null(err)) observed$err.ini <- observed[[err]]
- err = "err.irls"
- while (reweight.diff > reweight.tol & n.iter < reweight.max.iter) {
- n.iter <- n.iter + 1
- 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, lower = lower, ...)
- reweight.diff = sum((sqrt(fit$var_ms_unweighted) - sigma.old)^2)
- if (!quiet) {
- cat("Iteration", n.iter, "yields variance estimates:\n")
- print(signif(fit$var_ms_unweighted, 8))
- cat("Sum of squared differences to last variance estimates:",
- signif(reweight.diff, 2), "\n")
- }
- }
- }
-
- # We need to return some more data for summary and plotting
- fit$solution_type <- solution_type
- fit$transform_rates <- transform_rates
- fit$transform_fractions <- transform_fractions
-
- # We also need the model for summary and plotting
- fit$mkinmod <- mkinmod
-
- # We need data and predictions for summary and plotting
- fit$observed <- observed
- fit$obs_vars <- obs_vars
- fit$predicted <- mkin_wide_to_long(out_predicted, time = "time")
-
- # Collect initial parameter values in two dataframes
- fit$start <- data.frame(value = c(state.ini.optim,
- backtransform_odeparms(parms.optim, mod_vars,
- transform_rates = transform_rates,
- transform_fractions = transform_fractions)))
- fit$start$type = c(rep("state", length(state.ini.optim)),
- rep("deparm", length(parms.optim)))
- fit$start$transformed = c(state.ini.optim, parms.optim)
-
- fit$fixed <- data.frame(value = c(state.ini.fixed,
- backtransform_odeparms(parms.fixed, mod_vars,
- transform_rates = transform_rates,
- transform_fractions = transform_fractions)))
- fit$fixed$type = c(rep("state", length(state.ini.fixed)),
- rep("deparm", length(parms.fixed)))
-
- bparms.optim = backtransform_odeparms(fit$par, mod_vars,
- transform_rates = transform_rates,
- transform_fractions = transform_fractions)
- bparms.fixed = backtransform_odeparms(c(state.ini.fixed, parms.fixed),
- mod_vars,
- transform_rates = transform_rates,
- transform_fractions = transform_fractions)
- bparms.all = c(bparms.optim, bparms.fixed)
-
- # Collect observed, predicted, residuals and weighting
- data <- merge(fit$observed, fit$predicted, by = c("time", "name"))
- data$name <- ordered(data$name, levels = obs_vars)
- data <- data[order(data$name, data$time), ]
-
- fit$data <- data.frame(time = data$time,
- variable = data$name,
- observed = data$value.x,
- predicted = data$value.y)
- fit$data$residual <- fit$data$observed - fit$data$predicted
- if (!is.null(data$err.ini)) fit$data$err.ini <- data$err.ini
- if (!is.null(err)) fit$data[[err]] <- data[[err]]
-
- fit$atol <- atol
- fit$rtol <- rtol
- fit$weight.ini <- weight.ini
- fit$reweight.method <- reweight.method
- fit$reweight.tol <- reweight.tol
- fit$reweight.max.iter <- reweight.max.iter
-
- # Return all backtransformed parameters for summary
- fit$bparms.optim <- bparms.optim
- fit$bparms.fixed <- bparms.fixed
- fit$bparms.ode <- bparms.all[mkinmod$parms] # Return ode parameters for further fitting
- fit$date <- date()
-
- class(fit) <- c("mkinfit", "modFit")
- return(fit)
-}
-
-summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, ...) {
- param <- object$par
- pnames <- names(param)
- p <- length(param)
- mod_vars <- names(object$mkinmod$diffs)
- covar <- try(solve(0.5*object$hessian), silent = TRUE) # unscaled covariance
- rdf <- object$df.residual
- resvar <- object$ssr / rdf
- if (!is.numeric(covar)) {
- covar <- NULL
- se <- lci <- uci <- tval <- pval1 <- pval2 <- rep(NA, p)
- } else {
- rownames(covar) <- colnames(covar) <- pnames
- se <- sqrt(diag(covar) * resvar)
- lci <- param + qt(alpha/2, rdf) * se
- uci <- param + qt(1-alpha/2, rdf) * se
- tval <- param/se
- pval1 <- 2 * pt(abs(tval), rdf, lower.tail = FALSE)
- pval2 <- pt(abs(tval), rdf, lower.tail = FALSE)
- }
-
- names(se) <- pnames
- modVariance <- object$ssr / length(object$residuals)
-
- param <- cbind(param, se, lci, uci, tval, pval1, pval2)
- dimnames(param) <- list(pnames, c("Estimate", "Std. Error", "Lower", "Upper",
- "t value", "Pr(>|t|)", "Pr(>t)"))
-
- blci <- buci <- numeric()
- # Only transform boundaries of CI for one parameter at a time
- for (pname in pnames) {
- par.lower <- par.upper <- object$par
- par.lower[pname] <- param[pname, "Lower"]
- par.upper[pname] <- param[pname, "Upper"]
- blci[pname] <- backtransform_odeparms(par.lower, mod_vars,
- object$transform_rates, object$transform_fractions)[pname]
- buci[pname] <- backtransform_odeparms(par.upper, mod_vars,
- object$transform_rates, object$transform_fractions)[pname]
- }
- bparam <- cbind(object$bparms.optim, blci, buci)
- dimnames(bparam) <- list(pnames, c("Estimate", "Lower", "Upper"))
-
- ans <- list(
- version = as.character(packageVersion("mkin")),
- Rversion = paste(R.version$major, R.version$minor, sep="."),
- date.fit = object$date,
- date.summary = date(),
- solution_type = object$solution_type,
- use_of_ff = object$mkinmod$use_of_ff,
- weight.ini = object$weight.ini,
- reweight.method = object$reweight.method,
- residuals = object$residuals,
- residualVariance = resvar,
- sigma = sqrt(resvar),
- modVariance = modVariance,
- df = c(p, rdf),
- cov.unscaled = covar,
- cov.scaled = covar * resvar,
- info = object$info,
- niter = object$iterations,
- stopmess = message,
- par = param,
- bpar = bparam)
-
- ans$diffs <- object$mkinmod$diffs
- if(data) ans$data <- object$data
- ans$start <- object$start
-
- ans$fixed <- object$fixed
-
- ans$errmin <- mkinerrmin(object, alpha = 0.05)
-
- ans$bparms.ode <- object$bparms.ode
- ep <- endpoints(object)
- if (length(ep$ff) != 0)
- ans$ff <- ep$ff
- if(distimes) ans$distimes <- ep$distimes
- if(length(ep$SFORB) != 0) ans$SFORB <- ep$SFORB
- class(ans) <- c("summary.mkinfit", "summary.modFit")
- return(ans)
-}
-
-# Expanded from print.summary.modFit
-print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), ...) {
- cat("mkin version: ", x$version, "\n")
- cat("R version: ", x$Rversion, "\n")
- cat("Date of fit: ", x$date.fit, "\n")
- cat("Date of summary:", x$date.summary, "\n")
-
- cat("\nEquations:\n")
- print(noquote(as.character(x[["diffs"]])))
- df <- x$df
- rdf <- df[2]
-
- cat("\nMethod used for solution of differential equation system:\n")
- cat(x$solution_type, "\n")
-
- cat("\nWeighting:", x$weight.ini)
- if(!is.null(x$reweight.method)) cat(" then iterative reweighting method",
- x$reweight.method)
- cat("\n")
-
- cat("\nStarting values for optimised parameters:\n")
- print(x$start)
-
- cat("\nFixed parameter values:\n")
- if(length(x$fixed$value) == 0) cat("None\n")
- else print(x$fixed)
-
- cat("\nOptimised, transformed parameters:\n")
- print(signif(x$par, digits = digits))
-
- cat("\nBacktransformed parameters:\n")
- print(signif(x$bpar, digits = digits))
-
- cat("\nResidual standard error:",
- format(signif(x$sigma, digits)), "on", rdf, "degrees of freedom\n")
-
- cat("\nChi2 error levels in percent:\n")
- x$errmin$err.min <- 100 * x$errmin$err.min
- print(x$errmin, digits=digits,...)
-
- printdistimes <- !is.null(x$distimes)
- if(printdistimes){
- cat("\nEstimated disappearance times:\n")
- print(x$distimes, digits=digits,...)
- }
-
- printff <- !is.null(x$ff)
- if(printff & x$use_of_ff == "min"){
- cat("\nEstimated formation fractions:\n")
- print(data.frame(ff = x$ff), digits=digits,...)
- }
-
- printSFORB <- !is.null(x$SFORB)
- if(printSFORB){
- cat("\nEstimated Eigenvalues of SFORB model(s):\n")
- print(x$SFORB, digits=digits,...)
- }
-
- cat("\nParameter correlation:\n")
- if (!is.null(x$cov.unscaled)){
- Corr <- cov2cor(x$cov.unscaled)
- rownames(Corr) <- colnames(Corr) <- rownames(x$par)
- print(Corr, digits = digits, ...)
- } else {
- cat("Could not estimate covariance matrix; singular system:\n")
- }
-
- printdata <- !is.null(x$data)
- if (printdata){
- cat("\nData:\n")
- print(format(x$data, digits = digits, ...), row.names = FALSE)
- }
-
- invisible(x)
-}
-# vim: set ts=2 sw=2 expandtab:
diff --git a/R/mkinmod.R b/R/mkinmod.R deleted file mode 100644 index 2bc6ae3..0000000 --- a/R/mkinmod.R +++ /dev/null @@ -1,242 +0,0 @@ -# Copyright (C) 2010-2013 Johannes Ranke {{{
-# Contact: jranke@uni-bremen.de
-
-# This file is part of the R package mkin
-
-# mkin is free software: you can redistribute it and/or modify it under the
-# terms of the GNU General Public License as published by the Free Software
-# Foundation, either version 3 of the License, or (at your option) any later
-# version.
-
-# This program is distributed in the hope that it will be useful, but WITHOUT
-# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
-# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-# details.
-
-# You should have received a copy of the GNU General Public License along with
-# this program. If not, see <http://www.gnu.org/licenses/> }}}
-
-mkinmod <- function(..., use_of_ff = "min", speclist = NULL)
-{
- if (is.null(speclist)) spec <- list(...)
- else spec <- speclist
- obs_vars <- names(spec)
-
- # Check if any of the names of the observed variables contains any other
- for (obs_var in obs_vars) {
- if (length(grep(obs_var, obs_vars)) > 1) stop("Sorry, variable names can not contain each other")
- if (grepl("_to_", obs_var)) stop("Sorry, names of observed variables can not contain _to_")
- }
-
- if (!use_of_ff %in% c("min", "max"))
- stop("The use of formation fractions 'use_of_ff' can only be 'min' or 'max'")
-
- # The returned model will be a list of character vectors, containing {{{
- # differential equations, parameter names and a mapping from model variables
- # to observed variables. If possible, a matrix representation of the
- # differential equations is included
- parms <- vector()
- diffs <- vector()
- map <- list()
- # }}}
-
- # Give a warning when a model with time dependent degradation uses formation {{{
- # fractions
- if(spec[[1]]$type %in% c("FOMC", "DFOP", "HS")) {
- mat = FALSE
- } else mat = TRUE
- #}}}
-
- # Establish a list of differential equations as well as a map from observed {{{
- # compartments to differential equations
- for (varname in obs_vars)
- {
- # Check the type component of the compartment specification {{{
- if(is.null(spec[[varname]]$type)) stop(
- "Every part of the model specification must be a list containing a type component")
- if(!spec[[varname]]$type %in% c("SFO", "FOMC", "DFOP", "HS", "SFORB")) stop(
- "Available types are SFO, FOMC, DFOP, HS and SFORB only")
- if(spec[[varname]]$type %in% c("FOMC", "DFOP", "HS") & match(varname, obs_vars) != 1) {
- stop(paste("Types FOMC, DFOP and HS are only implemented for the first compartment,",
- "which is assumed to be the source compartment"))
- } #}}}
- # New (sub)compartments (boxes) needed for the model type {{{
- new_boxes <- switch(spec[[varname]]$type,
- SFO = varname,
- FOMC = varname,
- DFOP = varname,
- HS = varname,
- SFORB = paste(varname, c("free", "bound"), sep="_")
- )
- map[[varname]] <- new_boxes
- names(map[[varname]]) <- rep(spec[[varname]]$type, length(new_boxes)) #}}}
- # Start a new differential equation for each new box {{{
- new_diffs <- paste("d_", new_boxes, " =", sep="")
- names(new_diffs) <- new_boxes
- diffs <- c(diffs, new_diffs) #}}}
- } #}}}
-
- # Create content of differential equations and build parameter list {{{
- for (varname in obs_vars)
- {
- # Get the name of the box(es) we are working on for the decline term(s)
- box_1 = map[[varname]][[1]] # This is the only box unless type is SFORB
- # Turn on sink if this is not explicitly excluded by the user by
- # specifying sink=FALSE
- if(is.null(spec[[varname]]$sink)) spec[[varname]]$sink <- TRUE
- if(spec[[varname]]$type %in% c("SFO", "SFORB")) { # {{{ Add SFO or SFORB decline term
- if (use_of_ff == "min") { # Minimum use of formation fractions
- if(spec[[varname]]$sink) {
- # If sink is required, add first-order sink term
- k_compound_sink <- paste("k", box_1, "sink", sep="_")
- parms <- c(parms, k_compound_sink)
- decline_term <- paste(k_compound_sink, "*", box_1)
- } else { # otherwise no decline term needed here
- decline_term = "0"
- }
- } else {
- k_compound <- paste("k", box_1, sep="_")
- parms <- c(parms, k_compound)
- decline_term <- paste(k_compound, "*", box_1)
- }
- } #}}}
- if(spec[[varname]]$type == "FOMC") { # {{{ Add FOMC decline term
- # From p. 53 of the FOCUS kinetics report
- decline_term <- paste("(alpha/beta) * ((time/beta) + 1)^-1 *", box_1)
- parms <- c(parms, "alpha", "beta")
- } #}}}
- if(spec[[varname]]$type == "DFOP") { # {{{ Add DFOP decline term
- # From p. 57 of the FOCUS kinetics report
- decline_term <- paste("((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) *", box_1)
- parms <- c(parms, "k1", "k2", "g")
- } #}}}
- if(spec[[varname]]$type == "HS") { # {{{ Add HS decline term
- # From p. 55 of the FOCUS kinetics report
- decline_term <- paste("ifelse(time <= tb, k1, k2)", "*", box_1)
- parms <- c(parms, "k1", "k2", "tb")
- } #}}}
- # Add origin decline term to box 1 (usually the only box, unless type is SFORB)#{{{
- diffs[[box_1]] <- paste(diffs[[box_1]], "-", decline_term)#}}}
- if(spec[[varname]]$type == "SFORB") { # {{{ Add SFORB reversible binding terms
- box_2 = map[[varname]][[2]]
- if (use_of_ff == "min") { # Minimum use of formation fractions
- k_free_bound <- paste("k", varname, "free", "bound", sep="_")
- k_bound_free <- paste("k", varname, "bound", "free", sep="_")
- parms <- c(parms, k_free_bound, k_bound_free)
- reversible_binding_term_1 <- paste("-", k_free_bound, "*", box_1, "+",
- k_bound_free, "*", box_2)
- reversible_binding_term_2 <- paste("+", k_free_bound, "*", box_1, "-",
- k_bound_free, "*", box_2)
- } else { # Use formation fractions also for the free compartment
- stop("The maximum use of formation fractions is not supported for SFORB models")
- # The problems were: Calculation of dissipation times did not work in this case
- # and the coefficient matrix is not generated correctly by the code present
- # in this file in this case
- f_free_bound <- paste("f", varname, "free", "bound", sep="_")
- k_bound_free <- paste("k", varname, "bound", "free", sep="_")
- parms <- c(parms, f_free_bound, k_bound_free)
- reversible_binding_term_1 <- paste("+", k_bound_free, "*", box_2)
- reversible_binding_term_2 <- paste("+", f_free_bound, "*", k_compound, "*", box_1, "-",
- k_bound_free, "*", box_2)
- }
- diffs[[box_1]] <- paste(diffs[[box_1]], reversible_binding_term_1)
- diffs[[box_2]] <- paste(diffs[[box_2]], reversible_binding_term_2)
- } #}}}
-
- # Transfer between compartments#{{{
- to <- spec[[varname]]$to
- if(!is.null(to)) {
- # Name of box from which transfer takes place
- origin_box <- box_1
-
- # Add transfer terms to listed compartments
- for (target in to) {
- target_box <- switch(spec[[target]]$type,
- SFO = target,
- SFORB = paste(target, "free", sep="_"))
- if (use_of_ff == "min" && spec[[varname]]$type %in% c("SFO",
- "SFORB")) {
- k_from_to <- paste("k", origin_box, target_box, sep="_")
- parms <- c(parms, k_from_to)
- diffs[[origin_box]] <- paste(diffs[[origin_box]], "-",
- k_from_to, "*", origin_box)
- diffs[[target_box]] <- paste(diffs[[target_box]], "+",
- k_from_to, "*", origin_box)
- } else {
- if (!spec[[varname]]$sink) {
- stop("Turning off the sink when using formation fractions is not supported")
- }
- fraction_to_target = paste("f", origin_box, "to", target, sep="_")
- parms <- c(parms, fraction_to_target)
- diffs[[target_box]] <- paste(diffs[[target_box]], "+",
- fraction_to_target, "*", decline_term)
- }
- }
- } #}}}
- } #}}}
-
- model <- list(diffs = diffs, parms = parms, map = map, spec = spec, use_of_ff = use_of_ff)
-
- # Create coefficient matrix if appropriate#{{{
- if (mat) {
- boxes <- names(diffs)
- n <- length(boxes)
- m <- matrix(nrow=n, ncol=n, dimnames=list(boxes, boxes))
-
- if (use_of_ff == "min") { # {{{ Minimum use of formation fractions
- for (from in boxes) {
- for (to in boxes) {
- if (from == to) { # diagonal elements
- k.candidate = paste("k", from, c(boxes, "sink"), sep="_")
- k.candidate = sub("free.*bound", "free_bound", k.candidate)
- k.candidate = sub("bound.*free", "bound_free", k.candidate)
- k.effective = intersect(model$parms, k.candidate)
- m[from,to] = ifelse(length(k.effective) > 0,
- paste("-", k.effective, collapse = " "), "0")
-
- } else { # off-diagonal elements
- k.candidate = paste("k", from, to, sep="_")
- if (sub("_free$", "", from) == sub("_bound$", "", to)) {
- k.candidate = paste("k", sub("_free$", "_free_bound", from), sep="_")
- }
- if (sub("_bound$", "", from) == sub("_free$", "", to)) {
- k.candidate = paste("k", sub("_bound$", "_bound_free", from), sep="_")
- }
- k.effective = intersect(model$parms, k.candidate)
- m[to, from] = ifelse(length(k.effective) > 0,
- k.effective, "0")
- }
- }
- } # }}}
- } else { # {{{ Use formation fractions where possible
- for (from in boxes) {
- for (to in boxes) {
- if (from == to) { # diagonal elements
- k.candidate = paste("k", from, sep="_")
- m[from,to] = ifelse(k.candidate %in% model$parms,
- paste("-", k.candidate), "0")
- if(grepl("_free", from)) { # add transfer to bound compartment for SFORB
- m[from,to] = paste(m[from,to], "-", paste("k", from, "bound", sep="_"))
- }
- if(grepl("_bound", from)) { # add backtransfer to free compartment for SFORB
- m[from,to] = paste("- k", from, "free", sep="_")
- }
- m[from,to] = m[from,to]
- } else { # off-diagonal elements
- f.candidate = paste("f", from, "to", to, sep="_")
- k.candidate = paste("k", from, to, sep="_")
- # SFORB with maximum use of formation fractions not implemented, see above
- m[to, from] = ifelse(f.candidate %in% model$parms,
- paste(f.candidate, " * k_", from, sep=""),
- ifelse(k.candidate %in% model$parms, k.candidate, "0"))
- }
- }
- }
- } # }}}
- model$coefmat <- m
- }#}}}
-
- class(model) <- "mkinmod"
- return(model)
-}
-# vim: set foldmethod=marker ts=2 sw=2 expandtab:
diff --git a/R/mkinparplot.R b/R/mkinparplot.R deleted file mode 100644 index 28c1d2a..0000000 --- a/R/mkinparplot.R +++ /dev/null @@ -1,62 +0,0 @@ -# Copyright (C) 2014 Johannes Ranke -# Contact: jranke@uni-bremen.de - -# This file is part of the R package mkin - -# mkin is free software: you can redistribute it and/or modify it under the -# terms of the GNU General Public License as published by the Free Software -# Foundation, either version 3 of the License, or (at your option) any later -# version. - -# This program is distributed in the hope that it will be useful, but WITHOUT -# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more -# details. - -# You should have received a copy of the GNU General Public License along with -# this program. If not, see <http://www.gnu.org/licenses/> -mkinparplot <- function(object) { - state.optim = rownames(subset(object$start, type == "state")) - deparms.optim = rownames(subset(object$start, type == "deparm")) - fractions.optim = grep("^f_", deparms.optim, value = TRUE) - if ("g" %in% deparms.optim) fractions.optim <- c("g", fractions.optim) - rates.optim.unsorted = setdiff(deparms.optim, fractions.optim) - rates.optim <- rownames(object$start[rates.optim.unsorted, ]) - n.plot <- c(state.optim = length(state.optim), - rates.optim = length(rates.optim), - fractions.optim = length(fractions.optim)) - n.plot <- n.plot[n.plot > 0] - - oldpars <- par(no.readonly = TRUE) - layout(matrix(1:length(n.plot), ncol = 1), heights = n.plot + 1) - - s <- summary(object) - bpar <- data.frame(t(s$bpar)) - par(mar = c(2.1, 2.1, 0.1, 2.1)) - par(cex = 1) - for (type in names(n.plot)) { - parnames <- get(type) - values <- bpar[parnames] - xlim = switch(type, - state.optim = range(c(0, unlist(values)), - na.rm = TRUE, finite = TRUE), - rates.optim = range(c(0, unlist(values)), - na.rm = TRUE, finite = TRUE), - fractions.optim = range(c(0, 1, unlist(values)), - na.rm = TRUE, finite = TRUE)) - stripchart(values["Estimate", ][length(parnames):1], - xlim = xlim, - ylim = c(0.5, length(get(type)) + 0.5), - yaxt = "n") - if (type %in% c("rates.optim", "fractions.optim")) abline(v = 0, lty = 2) - if (type %in% c("fractions.optim")) abline(v = 1, lty = 2) - position <- ifelse(values["Estimate", ] < mean(xlim)/2, "right", "left") - text(ifelse(position == "left", min(xlim), max(xlim)), - length(parnames):1, parnames, - pos = ifelse(position == "left", 4, 2)) - arrows(as.numeric(values["Lower", ]), length(parnames):1, - as.numeric(values["Upper", ]), length(parnames):1, - code = 3, angle = 90, length = 0.05) - } - par(oldpars) -} diff --git a/R/mkinplot.R b/R/mkinplot.R deleted file mode 100644 index b9becfd..0000000 --- a/R/mkinplot.R +++ /dev/null @@ -1,4 +0,0 @@ -mkinplot <- function(fit, ...) -{ - plot(fit, ...) -} diff --git a/R/mkinpredict.R b/R/mkinpredict.R deleted file mode 100644 index 14efc56..0000000 --- a/R/mkinpredict.R +++ /dev/null @@ -1,121 +0,0 @@ -# Copyright (C) 2010-2014 Johannes Ranke -# Some lines in this code are copyright (C) 2013 Eurofins Regulatory AG -# Contact: jranke@uni-bremen.de - -# This file is part of the R package mkin - -# mkin is free software: you can redistribute it and/or modify it under the -# terms of the GNU General Public License as published by the Free Software -# Foundation, either version 3 of the License, or (at your option) any later -# version. - -# This program is distributed in the hope that it will be useful, but WITHOUT -# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more -# details. - -# You should have received a copy of the GNU General Public License along with -# this program. If not, see <http://www.gnu.org/licenses/> - -mkinpredict <- function(mkinmod, odeparms, odeini, - outtimes, solution_type = "deSolve", - method.ode = "lsoda", atol = 1e-8, rtol = 1e-10, - map_output = TRUE, ...) { - - # Get the names of the state variables in the model - mod_vars <- names(mkinmod$diffs) - - # Order the inital values for state variables if they are named - if (!is.null(names(odeini))) { - odeini <- odeini[mod_vars] - } - - # Create function for evaluation of expressions with ode parameters and initial values - evalparse <- function(string) - { - eval(parse(text=string), as.list(c(odeparms, odeini))) - } - - # Create a function calculating the differentials specified by the model - # if necessary - if (solution_type == "analytical") { - parent.type = names(mkinmod$map[[1]])[1] - parent.name = names(mkinmod$diffs)[[1]] - o <- switch(parent.type, - SFO = SFO.solution(outtimes, - evalparse(parent.name), - ifelse(mkinmod$use_of_ff == "min", - evalparse(paste("k", parent.name, "sink", sep="_")), - evalparse(paste("k", parent.name, sep="_")))), - FOMC = FOMC.solution(outtimes, - evalparse(parent.name), - evalparse("alpha"), evalparse("beta")), - DFOP = DFOP.solution(outtimes, - evalparse(parent.name), - evalparse("k1"), evalparse("k2"), - evalparse("g")), - HS = HS.solution(outtimes, - evalparse(parent.name), - evalparse("k1"), evalparse("k2"), - evalparse("tb")), - SFORB = SFORB.solution(outtimes, - evalparse(parent.name), - evalparse(paste("k", parent.name, "bound", sep="_")), - evalparse(paste("k", sub("free", "bound", parent.name), "free", sep="_")), - evalparse(paste("k", parent.name, "sink", sep="_"))) - ) - out <- data.frame(outtimes, o) - names(out) <- c("time", sub("_free", "", parent.name)) - } - if (solution_type == "eigen") { - coefmat.num <- matrix(sapply(as.vector(mkinmod$coefmat), evalparse), - nrow = length(mod_vars)) - e <- eigen(coefmat.num) - c <- solve(e$vectors, odeini) - f.out <- function(t) { - e$vectors %*% diag(exp(e$values * t), nrow=length(mod_vars)) %*% c - } - o <- matrix(mapply(f.out, outtimes), - nrow = length(mod_vars), ncol = length(outtimes)) - out <- data.frame(outtimes, t(o)) - names(out) <- c("time", mod_vars) - } - if (solution_type == "deSolve") { - mkindiff <- function(t, state, parms) { - - time <- t - diffs <- vector() - for (box in names(mkinmod$diffs)) - { - diffname <- paste("d", box, sep="_") - diffs[diffname] <- with(as.list(c(time, state, parms)), - eval(parse(text=mkinmod$diffs[[box]]))) - } - return(list(c(diffs))) - } - out <- ode( - y = odeini, - times = outtimes, - func = mkindiff, - parms = odeparms, - method = method.ode, - atol = atol, - rtol = rtol, - ... - ) - } - if (map_output) { - # Output transformation for models with unobserved compartments like SFORB - out_mapped <- data.frame(time = out[,"time"]) - for (var in names(mkinmod$map)) { - if((length(mkinmod$map[[var]]) == 1) || solution_type == "analytical") { - out_mapped[var] <- out[, var] - } else { - out_mapped[var] <- rowSums(out[, mkinmod$map[[var]]]) - } - } - return(out_mapped) - } else { - return(out) - } -} diff --git a/R/mkinresplot.R b/R/mkinresplot.R deleted file mode 100644 index 07bd7df..0000000 --- a/R/mkinresplot.R +++ /dev/null @@ -1,53 +0,0 @@ -# Copyright (C) 2008-2014 Johannes Ranke
-# Contact: jranke@uni-bremen.de
-
-# This file is part of the R package mkin
-
-# mkin is free software: you can redistribute it and/or modify it under the
-# terms of the GNU General Public License as published by the Free Software
-# Foundation, either version 3 of the License, or (at your option) any later
-# version.
-
-# This program is distributed in the hope that it will be useful, but WITHOUT
-# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
-# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-# details.
-
-# You should have received a copy of the GNU General Public License along with
-# this program. If not, see <http://www.gnu.org/licenses/>
-if(getRversion() >= '2.15.1') utils::globalVariables(c("variable", "residual"))
-
-mkinresplot <- function (object,
- obs_vars = names(object$mkinmod$map),
- xlab = "Time", ylab = "Residual",
- maxabs = "auto", legend= TRUE, lpos = "topright", ...)
-{
- obs_vars_all <- as.character(unique(object$data$variable))
-
- if (length(obs_vars) > 0){
- obs_vars <- intersect(obs_vars_all, obs_vars)
- } else obs_vars <- obs_vars_all
-
- residuals <- subset(object$data, variable %in% obs_vars, residual)
-
- if (maxabs == "auto") maxabs = max(abs(residuals), na.rm = TRUE)
-
- col_obs <- pch_obs <- 1:length(obs_vars)
- names(col_obs) <- names(pch_obs) <- obs_vars
-
- plot(0, xlab = xlab, ylab = ylab,
- xlim = c(0, 1.1 * max(object$data$time)),
- ylim = c(-1.2 * maxabs, 1.2 * maxabs), ...)
-
- for(obs_var in obs_vars){
- residuals_plot <- subset(object$data, variable == obs_var, c("time", "residual"))
- points(residuals_plot, pch = pch_obs[obs_var], col = col_obs[obs_var])
- }
-
- abline(h = 0, lty = 2)
-
- if (legend == TRUE) {
- legend(lpos, inset = c(0.05, 0.05), legend = obs_vars,
- col = col_obs[obs_vars], pch = pch_obs[obs_vars])
- }
-}
diff --git a/R/plot.mkinfit.R b/R/plot.mkinfit.R deleted file mode 100644 index 4132d96..0000000 --- a/R/plot.mkinfit.R +++ /dev/null @@ -1,98 +0,0 @@ -# Copyright (C) 2010-2014 Johannes Ranke -# Contact: jranke@uni-bremen.de - -# This file is part of the R package mkin - -# mkin is free software: you can redistribute it and/or modify it under the -# terms of the GNU General Public License as published by the Free Software -# Foundation, either version 3 of the License, or (at your option) any later -# version. - -# This program is distributed in the hope that it will be useful, but WITHOUT -# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more -# details. - -# You should have received a copy of the GNU General Public License along with -# this program. If not, see <http://www.gnu.org/licenses/> -if(getRversion() >= '2.15.1') utils::globalVariables(c("type", "variable", "observed")) - -plot.mkinfit <- function(x, fit = x, - obs_vars = names(fit$mkinmod$map), - xlab = "Time", ylab = "Observed", - xlim = range(fit$data$time), - ylim = "default", - col_obs = 1:length(fit$mkinmod$map), - pch_obs = col_obs, - lty_obs = rep(1, length(fit$mkinmod$map)), - add = FALSE, legend = !add, - show_residuals = FALSE, maxabs = "auto", - lpos = "topright", inset = c(0.05, 0.05), ...) -{ - if (add && show_residuals) stop("If adding to an existing plot we can not show residuals") - - if (ylim == "default") { - ylim = c(0, max(subset(fit$data, variable %in% obs_vars)$observed, na.rm = TRUE)) - } - - solution_type = fit$solution_type - parms.all <- c(fit$bparms.optim, fit$bparms.fixed) - - ininames <- c( - rownames(subset(fit$start, type == "state")), - rownames(subset(fit$fixed, type == "state"))) - odeini <- parms.all[ininames] - - # Order initial state variables - names(odeini) <- sub("_0", "", names(odeini)) - odeini <- odeini[names(fit$mkinmod$diffs)] - - outtimes <- seq(xlim[1], xlim[2], length.out=100) - - odenames <- c( - rownames(subset(fit$start, type == "deparm")), - rownames(subset(fit$fixed, type == "deparm"))) - odeparms <- parms.all[odenames] - - out <- mkinpredict(fit$mkinmod, odeparms, odeini, outtimes, - solution_type = solution_type, atol = fit$atol, rtol = fit$rtol) - - # Set up the plot if not to be added to an existing plot - if (add == FALSE) { - if (show_residuals) { - oldpar <- par(no.readonly = TRUE) - layout(matrix(c(1, 2), 2, 1), heights = c(2, 1.3)) - par(mar = c(3, 4, 4, 2) + 0.1) - } - plot(0, type="n", - xlim = xlim, ylim = ylim, - xlab = xlab, ylab = ylab, ...) - } - # Plot the data and model output - names(col_obs) <- names(pch_obs) <- names(lty_obs) <- names(fit$mkinmod$map) - for (obs_var in obs_vars) { - points(subset(fit$data, variable == obs_var, c(time, observed)), - pch = pch_obs[obs_var], col = col_obs[obs_var]) - } - matlines(out$time, out[obs_vars], col = col_obs[obs_vars], lty = lty_obs[obs_vars]) - if (legend == TRUE) { - legend(lpos, inset= inset, legend = obs_vars, - col = col_obs[obs_vars], pch = pch_obs[obs_vars], lty = lty_obs[obs_vars]) - } - # Show residuals if requested - if (show_residuals) { - par(mar = c(5, 4, 0, 2) + 0.1) - residuals <- subset(fit$data, variable %in% obs_vars, residual) - if (maxabs == "auto") maxabs = max(abs(residuals), na.rm = TRUE) - plot(0, type="n", - xlim = xlim, - ylim = c(-1.2 * maxabs, 1.2 * maxabs), - xlab = xlab, ylab = "Residuals") - for(obs_var in obs_vars){ - residuals_plot <- subset(fit$data, variable == obs_var, c("time", "residual")) - points(residuals_plot, pch = pch_obs[obs_var], col = col_obs[obs_var]) - } - abline(h = 0, lty = 2) - par(oldpar, no.readonly = TRUE) - } -} diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R deleted file mode 100644 index a0c302f..0000000 --- a/R/transform_odeparms.R +++ /dev/null @@ -1,101 +0,0 @@ -# Copyright (C) 2010-2014 Johannes Ranke
-# Contact: jranke@uni-bremen.de
-
-# This file is part of the R package mkin
-
-# mkin is free software: you can redistribute it and/or modify it under the
-# terms of the GNU General Public License as published by the Free Software
-# Foundation, either version 3 of the License, or (at your option) any later
-# version.
-
-# This program is distributed in the hope that it will be useful, but WITHOUT
-# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
-# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-# details.
-
-# You should have received a copy of the GNU General Public License along with
-# this program. If not, see <http://www.gnu.org/licenses/>
-
-transform_odeparms <- function(parms, mod_vars,
- transform_rates = TRUE,
- transform_fractions = TRUE)
-{
- # Set up container for transformed parameters
- transparms <- parms
-
- # Log transformation for rate constants if requested
- index_k <- grep("^k_", names(transparms))
- if (length(index_k) > 0) {
- if(transform_rates) transparms[index_k] <- log(parms[index_k])
- else transparms[index_k] <- parms[index_k]
- }
-
- # Go through state variables and apply isotropic logratio transformation if requested
- for (box in mod_vars) {
- indices_f <- grep(paste("^f", box, sep = "_"), names(parms))
- f_names <- grep(paste("^f", box, sep = "_"), names(parms), value = TRUE)
- n_paths <- length(indices_f)
- if (n_paths > 0) {
- f <- parms[indices_f]
- trans_f <- ilr(c(f, 1 - sum(f)))
- names(trans_f) <- f_names
- if(transform_fractions) transparms[indices_f] <- trans_f
- else transparms[indices_f] <- f
- }
- }
-
- # Transform rates, fractions and tb also for FOMC, DFOP and HS models if requested
- for (pname in c("alpha", "beta", "k1", "k2", "tb")) {
- 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])
- }
- }
- if (!is.na(parms["g"])) {
- g <- parms["g"]
- transparms["g"] <- ifelse(transform_fractions, ilr(c(g, 1 - g)), g)
- }
-
- return(transparms)
-}
-
-backtransform_odeparms <- function(transparms, mod_vars,
- transform_rates = TRUE,
- transform_fractions = TRUE)
-{
- # Set up container for backtransformed parameters
- parms <- transparms
-
- # Exponential transformation for rate constants
- index_k <- grep("^k_", names(parms))
- if (length(index_k) > 0) {
- if(transform_rates) parms[index_k] <- exp(transparms[index_k])
- else parms[index_k] <- transparms[index_k]
- }
-
- # Go through state variables and apply inverse isotropic logratio transformation
- for (box in mod_vars) {
- indices_f <- grep(paste("^f", box, sep = "_"), names(transparms))
- f_names <- grep(paste("^f", box, sep = "_"), names(transparms), value = TRUE)
- n_paths <- length(indices_f)
- if (n_paths > 0) {
- f <- invilr(transparms[indices_f])[1:n_paths] # We do not need the last component
- names(f) <- f_names
- if(transform_fractions) parms[indices_f] <- f
- else parms[indices_f] <- transparms[indices_f]
- }
- }
-
- # Transform parameters also for FOMC, DFOP and HS models
- for (pname in c("alpha", "beta", "k1", "k2", "tb")) {
- if (!is.na(transparms[pname])) {
- parms[pname] <- ifelse(transform_rates, exp(transparms[pname]), transparms[pname])
- }
- }
- if (!is.na(transparms["g"])) {
- g <- transparms["g"]
- parms["g"] <- ifelse(transform_fractions, invilr(g)[1], g)
- }
-
- return(parms)
-}
@@ -1,121 +1,23 @@ -# mkin +# gmkin -The R package **mkin** provides calculation routines for the analysis of -chemical degradation data, including <b>m</b>ulticompartment <b>kin</b>etics as -needed for modelling the formation and decline of transformation products, or -if several compartments are involved. +The R package **gmkin** provides a browser based graphical user interface (GUI) for +fitting kinetic models to chemical degradation data based on R package +[mkin](http://github.com/jranke/mkin). The GUI is based on the +[gWidgetsWWW2](http://github.com/jverzani/gWidgetsWWW2) package developed by +John Verzani. The GUI elements are created by the JavaScript library +ExtJS which is bundled with gWidgetsWWW2. ## Installation -You can install the latest released version from -[CRAN](http://cran.r-project.org/package=mkin) from within R: - -```s -install.packages('mkin') -``` - -A development version is usually available from [R-Forge](http://r-forge.r-project.org/R/?group_id=615): - -```s -install.packages('mkin', repos = 'http://r-forge.r-project.org') -``` - -If R-Forge is lacking behind or if you prefer, you can install directly from -github using the `devtools` package. - ```s require(devtools) -install_github("mkin", "jranke", quick = TRUE) # quick = TRUE avoids vignette rebuilds +install_github("gWidgetsWWW2", "jverzani") ``` -## Background - -In the regulatory evaluation of chemical substances like plant protection -products (pesticides), biocides and other chemicals, degradation data play an -important role. For the evaluation of pesticide degradation experiments, -detailed guidance and helpful tools have been developed as detailed in -'Credits and historical remarks' below. - -## Usage - -A very simple usage example would be - - library("mkin") - example_data = data.frame( - name = rep("parent", 9), - time = c(0, 1, 3, 7, 14, 28, 63, 91, 119), - value = c(85.1, 57.9, 29.9, 14.6, 9.7, 6.6, 4, 3.9, 0.6) - ) - SFO <- mkinmod(parent = list(type = "SFO")) - SFO.fit <- mkinfit(SFO, example_data) - summary(SFO.fit) - plot(SFO.fit) - -For more examples have a look at the examples provided in the -[`mkinfit`](http://kinfit.r-forge.r-project.org/mkin_static/mkinfit.html) -documentation -or the package vignettes referenced from the -[mkin package documentation page](http://kinfit.r-forge.r-project.org/mkin_static/index.html) - -## Features - -* Highly flexible model specification using - [`mkinmod`](http://kinfit.r-forge.r-project.org/mkin_static/mkinmod.html), - including equilibrium reactions and using the single first-order - reversible binding (SFORB) model, which will automatically create - two latent state variables for the observed variable. -* Model solution (forward modelling) in the function - [`mkinpredict`](http://kinfit.r-forge.r-project.org/mkin_static/mkinpredict.html) - is performed either using the analytical solution for the case of - parent only degradation, an eigenvalue based solution if only simple - first-order (SFO) or SFORB kinetics are used in the model, or - using a numeric solver from the `deSolve` package (default is `lsoda`). - These have decreasing efficiency, and are automatically chosen - by default. -* Model optimisation with - [`mkinfit`](http://kinfit.r-forge.r-project.org/mkin_static/mkinfit.html) - internally using the `modFit` function from the `FME` package, - which uses the least-squares Levenberg-Marquardt algorithm from - `minpack.lm` per default. -* Kinetic rate constants and kinetic formation fractions are transformed - internally using - [`transform_odeparms`](http://kinfit.r-forge.r-project.org/mkin_static/transform_odeparms.html) - so their estimators can more reasonably be expected to follow - a normal distribution. This has the side effect that no constraints - are needed in the optimisation. Thanks to René Lehmann for the nice - cooperation on this, especially the isotropic logration transformation - that is now used for the formation fractions. -* A side effect of this is that when parameter estimates are backtransformed - to match the model definition, confidence intervals calculated from - standard errors are also backtransformed to the correct scale, and will - not include meaningless values like negative rate constants or - formation fractions adding up to more than 1, which can not occur in - a single experiment with a single defined radiolabel position. -* Summary and plotting functions. The `summary` of an `mkinfit` object is in - fact a full report that should give enough information to be able to - approximately reproduce the fit with other tools. -* The chi-squared error level as defined in the FOCUS kinetics guidance - (see below) is calculated for each observed variable. -* I recently added iteratively reweighted least squares in a similar way - it is done in KinGUII and CAKE (see below). Simply add the argument - `reweight = "obs"` to your call to `mkinfit` and a separate variance - componenent for each of the observed variables will be optimised - in a second stage after the primary optimisation algorithm has converged. - -## GUI - -There is a graphical user interface that I consider useful for real work. It -depends on the gWidgetsWWW2 package from John Verzani which also lives on -github. Installing gWidgetsWWW2 yields a lot of warnings concerning overly -long path names. This is because the JavaScript library ExtJS is installed +Installing gWidgetsWWW2 yields a lot of warnings concerning overly long path +names. This is because the JavaScript library ExtJS is installed along with it which has lots of files with long paths to be installed. - -```s -require(devtools) -install_github("gWidgetsWWW2", "jverzani") -``` - You start the GUI from your R terminal with latest mkin installed as shown below. You may also want to adapt the browser that R starts (using `options(browser="/usr/bin/firefox")` on linux, or setting the default browser @@ -138,57 +40,15 @@ save(FOCUS_2006_Z_gmkin, file = "FOCUS_2006_gmkin_Z.RData") ```  - -## Credits and historical remarks - -`mkin` would not be possible without the underlying software stack consisting -of R and the packages [deSolve](http://cran.r-project.org/package=deSolve), -[minpack.lm](http://cran.r-project.org/package=minpack.lm) and -[FME](http://cran.r-project.org/package=FME), to say the least. - -It could not have been written without me being introduced to regulatory fate -modelling of pesticides by Adrian Gurney during my time at Harlan Laboratories -Ltd (formerly RCC Ltd). `mkin` greatly profits from and largely follows -the work done by the -[FOCUS Degradation Kinetics Workgroup](http://focus.jrc.ec.europa.eu/dk), -as detailed in ther guidance document from 2006, slightly updated in 2011. - -Also, it was inspired by the first version of KinGUI developed by -BayerCropScience, which is based on the MatLab runtime environment. - -The companion package [kinfit](http://kinfit.r-forge.r-project.org/kinfit_static/index.html) was [started in 2008](https://r-forge.r-project.org/scm/viewvc.php?view=rev&root=kinfit&revision=2) and -[first published on -CRAN](http://cran.r-project.org/src/contrib/Archive/kinfit/) on 01 May -2010. - -The first `mkin` code was [published on 11 May 2010](https://r-forge.r-project.org/scm/viewvc.php?view=rev&root=kinfit&revision=8) and the -[first CRAN version](http://cran.r-project.org/src/contrib/Archive/mkin) -on 18 May 2010. - -After this, Bayer has developed an R based successor to KinGUI named KinGUII -whose R code is based on `mkin`, but which added, amongst other refinements, a -closed source graphical user interface (GUI), iteratively reweighted least -squares (IRLS) optimisation of the variance for each of the observed -variables, and Markov Chain Monte Carlo (MCMC) simulation functionality, -similar to what is available e.g. in the `FME` package. - -Somewhat in parallel, Syngenta has sponsored the development of an `mkin` (and -KinGUII?) based GUI application called CAKE, which also adds IRLS and MCMC, is -more limited in the model formulation, but puts more weight on usability. -CAKE is available for download from the [CAKE -website](http://projects.tessella.com/cake), where you can also -find a zip archive of the R scripts derived from `mkin`, published under the GPL -license. - -Finally, there is -[KineticEval](http://github.com/zhenglei-gao/KineticEval), which contains -a further development of the scripts used for KinGUII, so the different tools -will hopefully be able to learn from each other in the future as well. - - -## Contribute -Contributions are welcome! Your [mkin fork](https://help.github.com/articles/fork-a-repo) is just a mouse click away... This git repository is now the -master branch, but I figured out how to merge changes in both directions, -thanks to [this blog entry](http://cameron.bracken.bz/git-with-r-forge) -by Cameron Bracken, so contributors from r-forge are welcome as well. +## Status and known issues + +- Starting the GUI takes some time. Once it is started, it is reasonably responsive. +- Do not delete the last dataset or the last model from the respective lists, + this is not supported. +- The fit list is not always updated when using Firefox version 28 on Windows. This + works with Firefox 29 and with Chrome. +- gmkin was developed in the hope that it will be useful. However, no warranty can be + given that it will meet your expectations. There may be bugs, so please be + careful, check your results for plausibility and use your own expertise to judge + yourself. @@ -1,13 +1 @@ -TODO for version 1.0 -- Think about what a user would expect from version 1.0 -- Support model definitions without pathway to sink in combination with formation fractions -- Complete the main package vignette named mkin to include a method description -- Calculate pseudoDT50 values as recommended by FOCUS -- Improve formatting of differential equations in the summary - -Nice to have: -- Calculate confidence intervals for DT50 and DT90 values when only one parameter is involved -- Calculate transformation only DT50 values (exclude pathways to sink) as additional information -- Calculate DT50 values from smaller rate constant/eigenvalue of DFOP, HS and SFORB models -- Check possibility to include a stripped extJS in gWidgetsWWW2, or to publish a separate - R package for extJS +- Make plot of fit and residuals more configurable diff --git a/data/FOCUS_2006_DFOP_ref_A_to_B.rda b/data/FOCUS_2006_DFOP_ref_A_to_B.rda Binary files differdeleted file mode 100644 index 8b9dc9f..0000000 --- a/data/FOCUS_2006_DFOP_ref_A_to_B.rda +++ /dev/null diff --git a/data/FOCUS_2006_FOMC_ref_A_to_F.rda b/data/FOCUS_2006_FOMC_ref_A_to_F.rda Binary files differdeleted file mode 100644 index 246a53f..0000000 --- a/data/FOCUS_2006_FOMC_ref_A_to_F.rda +++ /dev/null diff --git a/data/FOCUS_2006_HS_ref_A_to_F.rda b/data/FOCUS_2006_HS_ref_A_to_F.rda Binary files differdeleted file mode 100644 index 1e3b386..0000000 --- a/data/FOCUS_2006_HS_ref_A_to_F.rda +++ /dev/null diff --git a/data/FOCUS_2006_SFO_ref_A_to_F.rda b/data/FOCUS_2006_SFO_ref_A_to_F.rda Binary files differdeleted file mode 100644 index 411c536..0000000 --- a/data/FOCUS_2006_SFO_ref_A_to_F.rda +++ /dev/null diff --git a/data/FOCUS_2006_datasets.RData b/data/FOCUS_2006_datasets.RData Binary files differdeleted file mode 100644 index eb9e526..0000000 --- a/data/FOCUS_2006_datasets.RData +++ /dev/null diff --git a/data/mccall81_245T.RData b/data/mccall81_245T.RData Binary files differdeleted file mode 100644 index 530d83b..0000000 --- a/data/mccall81_245T.RData +++ /dev/null diff --git a/data/schaefer07_complex_case.RData b/data/schaefer07_complex_case.RData Binary files differdeleted file mode 100644 index d58aedd..0000000 --- a/data/schaefer07_complex_case.RData +++ /dev/null diff --git a/inst/GUI/README b/inst/GUI/README deleted file mode 100644 index 89e3726..0000000 --- a/inst/GUI/README +++ /dev/null @@ -1 +0,0 @@ -These code fragments are still experimental, but some basic functionality works. diff --git a/inst/GUI/TODO b/inst/GUI/TODO deleted file mode 100644 index 1fbcca6..0000000 --- a/inst/GUI/TODO +++ /dev/null @@ -1 +0,0 @@ -- Make plot of fit and residuals configurable diff --git a/inst/GUI/gmkin.R b/inst/GUI/gmkin.R index ae9d88d..7662443 100644 --- a/inst/GUI/gmkin.R +++ b/inst/GUI/gmkin.R @@ -19,7 +19,6 @@ # You should have received a copy of the GNU General Public License along with
# this program. If not, see <http://www.gnu.org/licenses/>
-require(mkin) # {{{1
# Set the GUI title and create the basic widget layout {{{1
w <- gwindow("gmkin - Browser based GUI for kinetic evaluations using mkin")
sb <- gstatusbar(paste("Powered by gWidgetsWWW2, ExtJS, Rook, FME, deSolve",
diff --git a/inst/unitTests/Makefile b/inst/unitTests/Makefile deleted file mode 100644 index 8d13225..0000000 --- a/inst/unitTests/Makefile +++ /dev/null @@ -1,15 +0,0 @@ -TOP=../.. -PKG=${shell cd ${TOP};pwd} -SUITE=doRUnit.R -R=R - -all: inst test - -inst: # Install package - cd ${TOP}/..;\ - ${R} CMD INSTALL ${PKG} - -test: # Run unit tests - export RCMDCHECK=FALSE;\ - cd ${TOP}/tests;\ - ${R} --vanilla --slave < ${SUITE} diff --git a/inst/unitTests/runit.mkinfit.R b/inst/unitTests/runit.mkinfit.R deleted file mode 100644 index 9a6bd72..0000000 --- a/inst/unitTests/runit.mkinfit.R +++ /dev/null @@ -1,294 +0,0 @@ -# $Id: runit.mkinfit.R 68 2010-09-09 22:40:04Z jranke $
-
-# Copyright (C) 2010-2013 Johannes Ranke
-# Contact: jranke@uni-bremen.de
-
-# This file is part of the R package mkin
-
-# mkin is free software: you can redistribute it and/or modify it under the
-# terms of the GNU General Public License as published by the Free Software
-# Foundation, either version 3 of the License, or (at your option) any later
-# version.
-
-# This program is distributed in the hope that it will be useful, but WITHOUT
-# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
-# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-# details.
-
-# You should have received a copy of the GNU General Public License along with
-# this program. If not, see <http://www.gnu.org/licenses/>
-
-# Test SFO model to a relative tolerance of 1% # {{{
-test.FOCUS_2006_SFO <- function()
-{
- SFO.1 <- mkinmod(parent = list(type = "SFO"))
- SFO.2 <- mkinmod(parent = list(type = "SFO"), use_of_ff = "max")
-
- fit.A.SFO.1 <- mkinfit(SFO.1, FOCUS_2006_A, quiet=TRUE)
- fit.A.SFO.2 <- mkinfit(SFO.2, FOCUS_2006_A, quiet=TRUE)
-
- median.A.SFO <- as.numeric(lapply(subset(FOCUS_2006_SFO_ref_A_to_F, dataset == "A",
- c(M0, k, DT50, DT90)), "median"))
-
- fit.A.SFO.1.r <- as.numeric(c(fit.A.SFO.1$bparms.optim, endpoints(fit.A.SFO.1)$distimes))
- dev.A.SFO.1 <- abs(round(100 * ((median.A.SFO - fit.A.SFO.1.r)/median.A.SFO), digits=1))
- checkIdentical(dev.A.SFO.1 < 1, rep(TRUE, length(dev.A.SFO.1)))
-
- fit.A.SFO.2.r <- as.numeric(c(fit.A.SFO.2$bparms.optim, endpoints(fit.A.SFO.2)$distimes))
- dev.A.SFO.2 <- abs(round(100 * ((median.A.SFO - fit.A.SFO.2.r)/median.A.SFO), digits=1))
- checkIdentical(dev.A.SFO.2 < 1, rep(TRUE, length(dev.A.SFO.2)))
-
- fit.C.SFO.1 <- mkinfit(SFO.1, FOCUS_2006_C, quiet=TRUE)
- fit.C.SFO.2 <- mkinfit(SFO.2, FOCUS_2006_C, quiet=TRUE)
-
- median.C.SFO <- as.numeric(lapply(subset(FOCUS_2006_SFO_ref_A_to_F, dataset == "C",
- c(M0, k, DT50, DT90)), "median"))
-
- fit.C.SFO.1.r <- as.numeric(c(fit.C.SFO.1$bparms.optim, endpoints(fit.C.SFO.1)$distimes))
- dev.C.SFO.1 <- abs(round(100 * ((median.C.SFO - fit.C.SFO.1.r)/median.C.SFO), digits=1))
- checkIdentical(dev.C.SFO.1 < 1, rep(TRUE, length(dev.C.SFO.1)))
-
- fit.C.SFO.2.r <- as.numeric(c(fit.C.SFO.2$bparms.optim, endpoints(fit.C.SFO.2)$distimes))
- dev.C.SFO.2 <- abs(round(100 * ((median.C.SFO - fit.C.SFO.2.r)/median.C.SFO), digits=1))
- checkIdentical(dev.C.SFO.2 < 1, rep(TRUE, length(dev.C.SFO.2)))
-} # }}}
-
-# Test FOMC model to a relative tolerance of 1% {{{
-# See kinfit vignette for a discussion of FOMC fits to FOCUS_2006_A
-# In this case, only M0, DT50 and DT90 are checked
-test.FOCUS_2006_FOMC <- function()
-{
- FOMC <- mkinmod(parent = list(type = "FOMC"))
-
- # FOCUS_2006_A (compare kinfit vignette for discussion)
- fit.A.FOMC <- mkinfit(FOMC, FOCUS_2006_A, quiet=TRUE)
-
- median.A.FOMC <- as.numeric(lapply(subset(FOCUS_2006_FOMC_ref_A_to_F, dataset == "A",
- c(M0, alpha, beta, DT50, DT90)), "median"))
-
- fit.A.FOMC.r <- as.numeric(c(fit.A.FOMC$bparms.optim, endpoints(fit.A.FOMC)$distimes))
- dev.A.FOMC <- abs(round(100 * ((median.A.FOMC - fit.A.FOMC.r)/median.A.FOMC), digits=1))
- dev.A.FOMC <- dev.A.FOMC[c(1, 4, 5)]
- checkIdentical(dev.A.FOMC < 1, rep(TRUE, length(dev.A.FOMC)))
-
- # FOCUS_2006_B
- fit.B.FOMC <- mkinfit(FOMC, FOCUS_2006_B, quiet=TRUE)
-
- median.B.FOMC <- as.numeric(lapply(subset(FOCUS_2006_FOMC_ref_A_to_F, dataset == "B",
- c(M0, alpha, beta, DT50, DT90)), "median"))
-
- fit.B.FOMC.r <- as.numeric(c(fit.B.FOMC$bparms.optim, endpoints(fit.B.FOMC)$distimes))
- dev.B.FOMC <- abs(round(100 * ((median.B.FOMC - fit.B.FOMC.r)/median.B.FOMC), digits=1))
- dev.B.FOMC <- dev.B.FOMC[c(1, 4, 5)]
- checkIdentical(dev.B.FOMC < 1, rep(TRUE, length(dev.B.FOMC)))
-
- # FOCUS_2006_C
- fit.C.FOMC <- mkinfit(FOMC, FOCUS_2006_C, quiet=TRUE)
-
- median.C.FOMC <- as.numeric(lapply(subset(FOCUS_2006_FOMC_ref_A_to_F, dataset == "C",
- c(M0, alpha, beta, DT50, DT90)), "median"))
-
- fit.C.FOMC.r <- as.numeric(c(fit.C.FOMC$bparms.optim, endpoints(fit.C.FOMC)$distimes))
- dev.C.FOMC <- abs(round(100 * ((median.C.FOMC - fit.C.FOMC.r)/median.C.FOMC), digits=1))
- dev.C.FOMC <- dev.C.FOMC[c(1, 4, 5)]
- checkIdentical(dev.C.FOMC < 1, rep(TRUE, length(dev.C.FOMC)))
-} # }}}
-
-# Test DFOP model, tolerance of 1% with the exception of f parameter for A {{{
-test.FOCUS_2006_DFOP <- function()
-{
- DFOP <- mkinmod(parent = list(type = "DFOP"))
-
- # FOCUS_2006_A
- fit.A.DFOP <- mkinfit(DFOP, FOCUS_2006_A, quiet=TRUE)
-
- median.A.DFOP <- as.numeric(lapply(subset(FOCUS_2006_DFOP_ref_A_to_B, dataset == "A",
- c(M0, k1, k2, f, DT50, DT90)), "median"))
-
- fit.A.DFOP.r <- as.numeric(c(fit.A.DFOP$bparms.optim, endpoints(fit.A.DFOP)$distimes))
- dev.A.DFOP <- abs(round(100 * ((median.A.DFOP - fit.A.DFOP.r)/median.A.DFOP), digits=1))
- # about 6.7% deviation for parameter f, the others are < 0.1%
- checkIdentical(dev.A.DFOP < c(1, 1, 1, 10, 1, 1), rep(TRUE, length(dev.A.DFOP)))
-
- # FOCUS_2006_B
- fit.B.DFOP <- mkinfit(DFOP, FOCUS_2006_B, quiet=TRUE)
-
- median.B.DFOP <- as.numeric(lapply(subset(FOCUS_2006_DFOP_ref_A_to_B, dataset == "B",
- c(M0, k1, k2, f, DT50, DT90)), "median"))
-
- fit.B.DFOP.r <- as.numeric(c(fit.B.DFOP$bparms.optim, endpoints(fit.B.DFOP)$distimes))
- dev.B.DFOP <- abs(round(100 * ((median.B.DFOP - fit.B.DFOP.r)/median.B.DFOP), digits=1))
- # about 0.6% deviation for parameter f, the others are <= 0.1%
- checkIdentical(dev.B.DFOP < 1, rep(TRUE, length(dev.B.DFOP)))
-} # }}}
-
-# Test HS model to a relative tolerance of 1% excluding Mathematica values {{{
-# as they are unreliable
-test.FOCUS_2006_HS <- function()
-{
- HS <- mkinmod(parent = list(type = "HS"))
-
- # FOCUS_2006_A
- fit.A.HS <- mkinfit(HS, FOCUS_2006_A, quiet=TRUE)
-
- median.A.HS <- as.numeric(lapply(subset(FOCUS_2006_HS_ref_A_to_F, dataset == "A",
- c(M0, k1, k2, tb, DT50, DT90)), "median"))
-
- fit.A.HS.r <- as.numeric(c(fit.A.HS$bparms.optim, endpoints(fit.A.HS)$distimes))
- dev.A.HS <- abs(round(100 * ((median.A.HS - fit.A.HS.r)/median.A.HS), digits=1))
- # about 6.7% deviation for parameter f, the others are < 0.1%
- checkIdentical(dev.A.HS < 1, rep(TRUE, length(dev.A.HS)))
-
- # FOCUS_2006_B
- fit.B.HS <- mkinfit(HS, FOCUS_2006_B, quiet=TRUE)
-
- median.B.HS <- as.numeric(lapply(subset(FOCUS_2006_HS_ref_A_to_F, dataset == "B",
- c(M0, k1, k2, tb, DT50, DT90)), "median"))
-
- fit.B.HS.r <- as.numeric(c(fit.B.HS$bparms.optim, endpoints(fit.B.HS)$distimes))
- dev.B.HS <- abs(round(100 * ((median.B.HS - fit.B.HS.r)/median.B.HS), digits=1))
- # < 10% deviation for M0, k1, DT50 and DT90, others are problematic
- dev.B.HS <- dev.B.HS[c(1, 2, 5, 6)]
- checkIdentical(dev.B.HS < 10, rep(TRUE, length(dev.B.HS)))
-
- # FOCUS_2006_C
- fit.C.HS <- mkinfit(HS, FOCUS_2006_C, quiet=TRUE)
-
- median.C.HS <- as.numeric(lapply(subset(FOCUS_2006_HS_ref_A_to_F, dataset == "C",
- c(M0, k1, k2, tb, DT50, DT90)), "median"))
-
- fit.A.HS.r <- as.numeric(c(fit.A.HS$bparms.optim, endpoints(fit.A.HS)$distimes))
- dev.A.HS <- abs(round(100 * ((median.A.HS - fit.A.HS.r)/median.A.HS), digits=1))
- # deviation <= 0.1%
- checkIdentical(dev.A.HS < 1, rep(TRUE, length(dev.A.HS)))
-} # }}}
-
-# Test SFORB model against DFOP solutions to a relative tolerance of 1% # {{{
-test.FOCUS_2006_SFORB <- function()
-{
- SFORB <- mkinmod(parent = list(type = "SFORB"))
-
- # FOCUS_2006_A
- fit.A.SFORB.1 <- mkinfit(SFORB, FOCUS_2006_A, quiet=TRUE)
- fit.A.SFORB.2 <- mkinfit(SFORB, FOCUS_2006_A, solution_type = "deSolve", quiet=TRUE)
-
- median.A.SFORB <- as.numeric(lapply(subset(FOCUS_2006_DFOP_ref_A_to_B, dataset == "A",
- c(M0, k1, k2, DT50, DT90)), "median"))
-
- fit.A.SFORB.1.r <- as.numeric(c(
- parent_0 = fit.A.SFORB.1$bparms.optim[[1]],
- k1 = endpoints(fit.A.SFORB.1)$SFORB[[1]],
- k2 = endpoints(fit.A.SFORB.1)$SFORB[[2]],
- endpoints(fit.A.SFORB.1)$distimes))
- dev.A.SFORB.1 <- abs(round(100 * ((median.A.SFORB - fit.A.SFORB.1.r)/median.A.SFORB), digits=1))
- # The first Eigenvalue is a lot different from k1 in the DFOP fit
- # The explanation is that the dataset is simply SFO
- dev.A.SFORB.1 <- dev.A.SFORB.1[c(1, 3, 4, 5)]
- checkIdentical(dev.A.SFORB.1 < 1, rep(TRUE, length(dev.A.SFORB.1)))
-
- fit.A.SFORB.2.r <- as.numeric(c(
- parent_0 = fit.A.SFORB.2$bparms.optim[[1]],
- k1 = endpoints(fit.A.SFORB.2)$SFORB[[1]],
- k2 = endpoints(fit.A.SFORB.2)$SFORB[[2]],
- endpoints(fit.A.SFORB.2)$distimes))
- dev.A.SFORB.2 <- abs(round(100 * ((median.A.SFORB - fit.A.SFORB.2.r)/median.A.SFORB), digits=1))
- # The first Eigenvalue is a lot different from k1 in the DFOP fit
- # The explanation is that the dataset is simply SFO
- dev.A.SFORB.2 <- dev.A.SFORB.2[c(1, 3, 4, 5)]
- checkIdentical(dev.A.SFORB.2 < 1, rep(TRUE, length(dev.A.SFORB.2)))
-
- # FOCUS_2006_B
- fit.B.SFORB.1 <- mkinfit(SFORB, FOCUS_2006_B, quiet=TRUE)
- fit.B.SFORB.2 <- mkinfit(SFORB, FOCUS_2006_B, solution_type = "deSolve", quiet=TRUE)
-
- median.B.SFORB <- as.numeric(lapply(subset(FOCUS_2006_DFOP_ref_A_to_B, dataset == "B",
- c(M0, k1, k2, DT50, DT90)), "median"))
-
- fit.B.SFORB.1.r <- as.numeric(c(
- parent_0 = fit.B.SFORB.1$bparms.optim[[1]],
- k1 = endpoints(fit.B.SFORB.1)$SFORB[[1]],
- k2 = endpoints(fit.B.SFORB.1)$SFORB[[2]],
- endpoints(fit.B.SFORB.1)$distimes))
- dev.B.SFORB.1 <- abs(round(100 * ((median.B.SFORB - fit.B.SFORB.1.r)/median.B.SFORB), digits=1))
- checkIdentical(dev.B.SFORB.1 < 1, rep(TRUE, length(dev.B.SFORB.1)))
-
- fit.B.SFORB.2.r <- as.numeric(c(
- parent_0 = fit.B.SFORB.2$bparms.optim[[1]],
- k1 = endpoints(fit.B.SFORB.2)$SFORB[[1]],
- k2 = endpoints(fit.B.SFORB.2)$SFORB[[2]],
- endpoints(fit.B.SFORB.2)$distimes))
- dev.B.SFORB.2 <- abs(round(100 * ((median.B.SFORB - fit.B.SFORB.2.r)/median.B.SFORB), digits=1))
- 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()
-{
- schaefer07_complex_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"))
-
- fit <- mkinfit(schaefer07_complex_model,
- mkin_wide_to_long(schaefer07_complex_case, time = "time"))
- s <- summary(fit)
- r <- schaefer07_complex_results
- attach(as.list(fit$bparms.optim))
- k_parent <- sum(k_parent_A1, k_parent_B1, k_parent_C1)
- r$mkin <- c(
- k_parent,
- s$distimes["parent", "DT50"],
- s$ff["parent_A1"],
- sum(k_A1_sink, k_A1_A2),
- s$distimes["A1", "DT50"],
- s$ff["parent_B1"],
- k_B1_sink,
- s$distimes["B1", "DT50"],
- s$ff["parent_C1"],
- k_C1_sink,
- s$distimes["C1", "DT50"],
- s$ff["A1_A2"],
- k_A2_sink,
- s$distimes["A2", "DT50"])
- r$means <- (r$KinGUI + r$ModelMaker)/2
- r$mkin.deviation <- abs(round(100 * ((r$mkin - r$means)/r$means), digits=1))
- checkIdentical(r$mkin.deviation[1:11] < 10, rep(TRUE, 11))
-} # }}}
-
-# vim: set foldmethod=marker ts=2 sw=2 expandtab:
diff --git a/inst/unitTests/runit.mkinpredict.R b/inst/unitTests/runit.mkinpredict.R deleted file mode 100644 index 997857c..0000000 --- a/inst/unitTests/runit.mkinpredict.R +++ /dev/null @@ -1,86 +0,0 @@ -# $Id: jranke $
-
-# Copyright (C) 2012 Johannes Ranke
-# Contact: jranke@uni-bremen.de
-
-# This file is part of the R package mkin
-
-# mkin is free software: you can redistribute it and/or modify it under the
-# terms of the GNU General Public License as published by the Free Software
-# Foundation, either version 3 of the License, or (at your option) any later
-# version.
-
-# This program is distributed in the hope that it will be useful, but WITHOUT
-# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
-# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-# details.
-
-# You should have received a copy of the GNU General Public License along with
-# this program. If not, see <http://www.gnu.org/licenses/>
-
-# Check solution types for SFO {{{
-test.SFO_solution_types <- function()
-{
- ot = seq(0, 100, by = 1)
- SFO <- mkinmod(parent = list(type = "SFO"))
- SFO.analytical <- round(subset(mkinpredict(SFO, c(k_parent_sink = 0.1),
- c(parent = 100), ot, solution_type = "analytical"), time == 100), digits=5)
- SFO.deSolve <- round(subset(mkinpredict(SFO, c(k_parent_sink = 0.1),
- c(parent = 100), ot, solution_type = "deSolve"), time == 100), digits=5)
- SFO.eigen <- round(subset(mkinpredict(SFO, c(k_parent_sink = 0.1),
- c(parent = 100), ot, solution_type = "eigen"), time == 100), digits=5)
-
- checkEquals(SFO.analytical, SFO.deSolve)
- checkEquals(SFO.analytical, SFO.eigen)
-} # }}}
-
-# Check model specification and solution types for SFO_SFO {{{
-# Relative Tolerance is 0.01%
-# Do not use time 0, as eigenvalue based solution does not give 0 at time 0 for metabolites
-# and relative tolerance is thus not met
-test.SFO_solution_types <- function()
-{
- tol = 0.01
- 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")
-
- ot = seq(0, 100, by = 1)
- r.1.e <- subset(mkinpredict(SFO_SFO.1,
- c(k_parent_m1 = 0.1, k_parent_sink = 0.1, k_m1_sink = 0.1),
- c(parent = 100, m1 = 0), ot, solution_type = "eigen"),
- time %in% c(1, 10, 50, 100))
- r.1.d <- subset(mkinpredict(SFO_SFO.1,
- c(k_parent_m1 = 0.1, k_parent_sink = 0.1, k_m1_sink = 0.1),
- c(parent = 100, m1 = 0), ot, solution_type = "deSolve"),
- time %in% c(1, 10, 50, 100))
-
- r.2.e <- subset(mkinpredict(SFO_SFO.2, c(k_parent = 0.2, f_parent_to_m1 = 0.5, k_m1 = 0.1),
- c(parent = 100, m1 = 0), ot, solution_type = "eigen"),
- time %in% c(1, 10, 50, 100))
- r.2.d <- subset(mkinpredict(SFO_SFO.2, c(k_parent = 0.2, f_parent_to_m1 = 0.5, k_m1 = 0.1),
- c(parent = 100, m1 = 0), ot, solution_type = "deSolve"),
- time %in% c(1, 10, 50, 100))
-
- # Compare eigen and deSolve for minimum use of formation fractions
- dev.1.e_d.percent = 100 * (r.1.e[-1] - r.1.d[-1])/r.1.e[-1]
- dev.1.e_d.percent = as.numeric(unlist((dev.1.e_d.percent)))
- dev.1.e_d.percent = ifelse(is.na(dev.1.e_d.percent), 0, dev.1.e_d.percent)
- checkIdentical(dev.1.e_d.percent < tol, rep(TRUE, length(dev.1.e_d.percent)))
-
- # Compare eigen and deSolve for maximum use of formation fractions
- dev.2.e_d.percent = 100 * (r.1.e[-1] - r.1.d[-1])/r.1.e[-1]
- dev.2.e_d.percent = as.numeric(unlist((dev.2.e_d.percent)))
- dev.2.e_d.percent = ifelse(is.na(dev.2.e_d.percent), 0, dev.2.e_d.percent)
- checkIdentical(dev.2.e_d.percent < tol, rep(TRUE, length(dev.2.e_d.percent)))
-
- # Compare minimum and maximum use of formation fractions
- dev.1_2.e.percent = 100 * (r.1.e[-1] - r.2.e[-1])/r.1.e[-1]
- dev.1_2.e.percent = as.numeric(unlist((dev.1_2.e.percent)))
- dev.1_2.e.percent = ifelse(is.na(dev.1_2.e.percent), 0, dev.1_2.e.percent)
- checkIdentical(dev.1_2.e.percent < tol, rep(TRUE, length(dev.1_2.e.percent)))
-
-} # }}}
-
-# vim: set foldmethod=marker ts=2 sw=2 expandtab:
diff --git a/man/DFOP.solution.Rd b/man/DFOP.solution.Rd deleted file mode 100644 index d30cf7f..0000000 --- a/man/DFOP.solution.Rd +++ /dev/null @@ -1,36 +0,0 @@ -\name{DFOP.solution}
-\Rdversion{1.1}
-\alias{DFOP.solution}
-\title{
-Dual First-Order in Parallel kinetics
-}
-\description{
- Function describing decline from a defined starting value using the sum
- of two exponential decline functions.
-}
-\usage{
-DFOP.solution(t, parent.0, k1, k2, g)
-}
-\arguments{
- \item{t}{ Time. }
- \item{parent.0}{ Starting value for the response variable at time zero. }
- \item{k1}{ First kinetic constant. }
- \item{k2}{ Second kinetic constant. }
- \item{g}{ Fraction of the starting value declining according to the
- first kinetic constant.
- }
-}
-\value{
- The value of the response variable at time \code{t}.
-}
-\references{
- FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence and
- Degradation Kinetics from Environmental Fate Studies on Pesticides in EU
- Registration} Report of the FOCUS Work Group on Degradation Kinetics,
- EC Document Reference Sanco/10058/2005 version 2.0, 434 pp,
- \url{http://focus.jrc.ec.europa.eu/dk}
-}
-\examples{
- \dontrun{plot(function(x) DFOP.solution(x, 100, 5, 0.5, 0.3), 0, 4, ylim=c(0,100))}
-}
-\keyword{ manip }
diff --git a/man/FOCUS_2006_DFOP_ref_A_to_B.Rd b/man/FOCUS_2006_DFOP_ref_A_to_B.Rd deleted file mode 100644 index 88bd4ac..0000000 --- a/man/FOCUS_2006_DFOP_ref_A_to_B.Rd +++ /dev/null @@ -1,39 +0,0 @@ -\name{FOCUS_2006_DFOP_ref_A_to_B}
-\Rdversion{1.1}
-\alias{FOCUS_2006_DFOP_ref_A_to_B}
-\docType{data}
-\title{
-Results of fitting the DFOP model to Datasets A to B of FOCUS (2006)
-}
-\description{
-A table with the fitted parameters and the resulting DT50 and DT90 values
-generated with different software packages. Taken directly from FOCUS (2006).
-The results from fitting the data with the Topfit software was removed, as
-the initial concentration of the parent compound was fixed to a value of 100
-in this fit.
-}
-\usage{data(FOCUS_2006_DFOP_ref_A_to_B)}
-\format{
- A data frame containing the following variables.
- \describe{
- \item{\code{package}}{a factor giving the name of the software package}
- \item{\code{M0}}{The fitted initial concentration of the parent compound}
- \item{\code{f}}{The fitted f parameter}
- \item{\code{k1}}{The fitted k1 parameter}
- \item{\code{k2}}{The fitted k2 parameter}
- \item{\code{DT50}}{The resulting half-life of the parent compound}
- \item{\code{DT90}}{The resulting DT90 of the parent compound}
- \item{\code{dataset}}{The FOCUS dataset that was used}
- }
-}
-\source{
- FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence and
- Degradation Kinetics from Environmental Fate Studies on Pesticides in EU
- Registration} Report of the FOCUS Work Group on Degradation Kinetics,
- EC Document Reference Sanco/10058/2005 version 2.0, 434 pp,
- \url{http://focus.jrc.ec.europa.eu/dk}
-}
-\examples{
-data(FOCUS_2006_DFOP_ref_A_to_B)
-}
-\keyword{datasets}
diff --git a/man/FOCUS_2006_FOMC_ref_A_to_F.Rd b/man/FOCUS_2006_FOMC_ref_A_to_F.Rd deleted file mode 100644 index 2fcc2db..0000000 --- a/man/FOCUS_2006_FOMC_ref_A_to_F.Rd +++ /dev/null @@ -1,38 +0,0 @@ -\name{FOCUS_2006_FOMC_ref_A_to_F}
-\Rdversion{1.1}
-\alias{FOCUS_2006_FOMC_ref_A_to_F}
-\docType{data}
-\title{
-Results of fitting the FOMC model to Datasets A to F of FOCUS (2006)
-}
-\description{
-A table with the fitted parameters and the resulting DT50 and DT90 values
-generated with different software packages. Taken directly from FOCUS (2006).
-The results from fitting the data with the Topfit software was removed, as
-the initial concentration of the parent compound was fixed to a value of 100
-in this fit.
-}
-\usage{data(FOCUS_2006_FOMC_ref_A_to_F)}
-\format{
- A data frame containing the following variables.
- \describe{
- \item{\code{package}}{a factor giving the name of the software package}
- \item{\code{M0}}{The fitted initial concentration of the parent compound}
- \item{\code{alpha}}{The fitted alpha parameter}
- \item{\code{beta}}{The fitted beta parameter}
- \item{\code{DT50}}{The resulting half-life of the parent compound}
- \item{\code{DT90}}{The resulting DT90 of the parent compound}
- \item{\code{dataset}}{The FOCUS dataset that was used}
- }
-}
-\source{
- FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence and
- Degradation Kinetics from Environmental Fate Studies on Pesticides in EU
- Registration} Report of the FOCUS Work Group on Degradation Kinetics,
- EC Document Reference Sanco/10058/2005 version 2.0, 434 pp,
- \url{http://focus.jrc.ec.europa.eu/dk}
-}
-\examples{
-data(FOCUS_2006_FOMC_ref_A_to_F)
-}
-\keyword{datasets}
diff --git a/man/FOCUS_2006_HS_ref_A_to_F.Rd b/man/FOCUS_2006_HS_ref_A_to_F.Rd deleted file mode 100644 index 6fc9993..0000000 --- a/man/FOCUS_2006_HS_ref_A_to_F.Rd +++ /dev/null @@ -1,39 +0,0 @@ -\name{FOCUS_2006_HS_ref_A_to_F}
-\Rdversion{1.1}
-\alias{FOCUS_2006_HS_ref_A_to_F}
-\docType{data}
-\title{
-Results of fitting the HS model to Datasets A to F of FOCUS (2006)
-}
-\description{
-A table with the fitted parameters and the resulting DT50 and DT90 values
-generated with different software packages. Taken directly from FOCUS (2006).
-The results from fitting the data with the Topfit software was removed, as
-the initial concentration of the parent compound was fixed to a value of 100
-in this fit.
-}
-\usage{data(FOCUS_2006_HS_ref_A_to_F)}
-\format{
- A data frame containing the following variables.
- \describe{
- \item{\code{package}}{a factor giving the name of the software package}
- \item{\code{M0}}{The fitted initial concentration of the parent compound}
- \item{\code{tb}}{The fitted tb parameter}
- \item{\code{k1}}{The fitted k1 parameter}
- \item{\code{k2}}{The fitted k2 parameter}
- \item{\code{DT50}}{The resulting half-life of the parent compound}
- \item{\code{DT90}}{The resulting DT90 of the parent compound}
- \item{\code{dataset}}{The FOCUS dataset that was used}
- }
-}
-\source{
- FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence and
- Degradation Kinetics from Environmental Fate Studies on Pesticides in EU
- Registration} Report of the FOCUS Work Group on Degradation Kinetics,
- EC Document Reference Sanco/10058/2005 version 2.0, 434 pp,
- \url{http://focus.jrc.ec.europa.eu/dk}
-}
-\examples{
-data(FOCUS_2006_HS_ref_A_to_F)
-}
-\keyword{datasets}
diff --git a/man/FOCUS_2006_SFO_ref_A_to_F.Rd b/man/FOCUS_2006_SFO_ref_A_to_F.Rd deleted file mode 100644 index 19650ed..0000000 --- a/man/FOCUS_2006_SFO_ref_A_to_F.Rd +++ /dev/null @@ -1,37 +0,0 @@ -\name{FOCUS_2006_SFO_ref_A_to_F}
-\Rdversion{1.1}
-\alias{FOCUS_2006_SFO_ref_A_to_F}
-\docType{data}
-\title{
-Results of fitting the SFO model to Datasets A to F of FOCUS (2006)
-}
-\description{
-A table with the fitted parameters and the resulting DT50 and DT90 values
-generated with different software packages. Taken directly from FOCUS (2006).
-The results from fitting the data with the Topfit software was removed, as
-the initial concentration of the parent compound was fixed to a value of 100
-in this fit.
-}
-\usage{data(FOCUS_2006_SFO_ref_A_to_F)}
-\format{
- A data frame containing the following variables.
- \describe{
- \item{\code{package}}{a factor giving the name of the software package}
- \item{\code{M0}}{The fitted initial concentration of the parent compound}
- \item{\code{k}}{The fitted first-order degradation rate constant}
- \item{\code{DT50}}{The resulting half-life of the parent compound}
- \item{\code{DT90}}{The resulting DT90 of the parent compound}
- \item{\code{dataset}}{The FOCUS dataset that was used}
- }
-}
-\source{
- FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence and
- Degradation Kinetics from Environmental Fate Studies on Pesticides in EU
- Registration} Report of the FOCUS Work Group on Degradation Kinetics,
- EC Document Reference Sanco/10058/2005 version 2.0, 434 pp,
- \url{http://focus.jrc.ec.europa.eu/dk}
-}
-\examples{
-data(FOCUS_2006_SFO_ref_A_to_F)
-}
-\keyword{datasets}
diff --git a/man/FOCUS_2006_datasets.Rd b/man/FOCUS_2006_datasets.Rd deleted file mode 100644 index 5053b88..0000000 --- a/man/FOCUS_2006_datasets.Rd +++ /dev/null @@ -1,35 +0,0 @@ -\name{FOCUS_2006_datasets}
-\Rdversion{1.1}
-\alias{FOCUS_2006_A}
-\alias{FOCUS_2006_B}
-\alias{FOCUS_2006_C}
-\alias{FOCUS_2006_D}
-\alias{FOCUS_2006_E}
-\alias{FOCUS_2006_F}
-\docType{data}
-\title{
-Datasets A to F from the FOCUS Kinetics report from 2006
-}
-\description{
-Data taken from FOCUS (2006), p. 258.
-}
-\usage{FOCUS_2006_datasets}
-\format{
- 6 datasets with observations on the following variables.
- \describe{
- \item{\code{name}}{a factor containing the name of the observed variable}
- \item{\code{time}}{a numeric vector containing time points}
- \item{\code{value}}{a numeric vector containing concentrations in percent of applied radioactivity}
- }
-}
-\source{
- FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence and
- Degradation Kinetics from Environmental Fate Studies on Pesticides in EU
- Registration} Report of the FOCUS Work Group on Degradation Kinetics,
- EC Document Reference Sanco/10058/2005 version 2.0, 434 pp,
- \url{http://focus.jrc.ec.europa.eu/dk}
-}
-\examples{
-FOCUS_2006_C
-}
-\keyword{datasets}
diff --git a/man/FOMC.solution.Rd b/man/FOMC.solution.Rd deleted file mode 100644 index d04d34e..0000000 --- a/man/FOMC.solution.Rd +++ /dev/null @@ -1,49 +0,0 @@ -\name{FOMC.solution}
-\Rdversion{1.1}
-\alias{FOMC.solution}
-\title{ First-Order Multi-Compartment kinetics }
-\description{
- Function describing exponential decline from a defined starting value, with
- a decreasing rate constant.
-
- The form given here differs slightly from the original reference by Gustafson
- and Holden (1990). The parameter \code{beta} corresponds to 1/beta in the
- original equation.
-}
-\usage{
-FOMC.solution(t, parent.0, alpha, beta)
-}
-\arguments{
- \item{t}{ Time. }
- \item{parent.0}{ Starting value for the response variable at time zero. }
- \item{alpha}{
- Shape parameter determined by coefficient of variation of rate constant
- values. }
- \item{beta}{
- Location parameter.
-}
-}
-\note{
- The solution of the FOMC kinetic model reduces to the
- \code{\link{SFO.solution}} for large values of \code{alpha} and
- \code{beta} with
- \eqn{k = \frac{\beta}{\alpha}}{k = beta/alpha}.
-}
-\value{
- The value of the response variable at time \code{t}.
-}
-\references{
- FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence and
- Degradation Kinetics from Environmental Fate Studies on Pesticides in EU
- Registration} Report of the FOCUS Work Group on Degradation Kinetics,
- EC Document Reference Sanco/10058/2005 version 2.0, 434 pp,
- \url{http://focus.jrc.ec.europa.eu/dk}
-
- Gustafson DI and Holden LR (1990) Nonlinear pesticide dissipation in soil: A
- new model based on spatial variability. \emph{Environmental Science and
- Technology} \bold{24}, 1032-1038
-}
-\examples{
- \dontrun{plot(function(x) FOMC.solution(x, 100, 10, 2), 0, 2)}
-}
-\keyword{ manip }
diff --git a/man/HS.solution.Rd b/man/HS.solution.Rd deleted file mode 100644 index 71f68a1..0000000 --- a/man/HS.solution.Rd +++ /dev/null @@ -1,34 +0,0 @@ -\name{HS.solution}
-\Rdversion{1.1}
-\alias{HS.solution}
-\title{ Hockey-Stick kinetics }
-\description{
- Function describing two exponential decline functions with a break point
- between them.
-}
-\usage{
-HS.solution(t, parent.0, k1, k2, tb)
-}
-\arguments{
- \item{t}{ Time. }
- \item{parent.0}{ Starting value for the response variable at time zero. }
- \item{k1}{ First kinetic constant. }
- \item{k2}{ Second kinetic constant. }
- \item{tb}{ Break point. Before this time, exponential decline according
- to \code{k1} is calculated, after this time, exponential decline proceeds
- according to \code{k2}. }
-}
-\value{
- The value of the response variable at time \code{t}.
-}
-\references{
- FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence and
- Degradation Kinetics from Environmental Fate Studies on Pesticides in EU
- Registration} Report of the FOCUS Work Group on Degradation Kinetics,
- EC Document Reference Sanco/10058/2005 version 2.0, 434 pp,
- \url{http://focus.jrc.ec.europa.eu/dk}
-}
-\examples{
- \dontrun{plot(function(x) HS.solution(x, 100, 2, 0.3, 0.5), 0, 2, ylim=c(0,100))}
-}
-\keyword{ manip }
diff --git a/man/SFO.solution.Rd b/man/SFO.solution.Rd deleted file mode 100644 index 41a9ba9..0000000 --- a/man/SFO.solution.Rd +++ /dev/null @@ -1,29 +0,0 @@ -\name{SFO.solution}
-\Rdversion{1.1}
-\alias{SFO.solution}
-\title{ Single First-Order kinetics }
-\description{
- Function describing exponential decline from a defined starting value.
-}
-\usage{
- SFO.solution(t, parent.0, k)
-}
-\arguments{
- \item{t}{ Time. }
- \item{parent.0}{ Starting value for the response variable at time zero. }
- \item{k}{ Kinetic constant. }
-}
-\value{
- The value of the response variable at time \code{t}.
-}
-\references{
- FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence and
- Degradation Kinetics from Environmental Fate Studies on Pesticides in EU
- Registration} Report of the FOCUS Work Group on Degradation Kinetics,
- EC Document Reference Sanco/10058/2005 version 2.0, 434 pp,
- \url{http://focus.jrc.ec.europa.eu/dk}
-}
-\examples{
- \dontrun{plot(function(x) SFO.solution(x, 100, 3), 0, 2)}
-}
-\keyword{ manip }
diff --git a/man/SFORB.solution.Rd b/man/SFORB.solution.Rd deleted file mode 100644 index a935f69..0000000 --- a/man/SFORB.solution.Rd +++ /dev/null @@ -1,36 +0,0 @@ -\name{SFORB.solution}
-\Rdversion{1.1}
-\alias{SFORB.solution}
-\title{ Single First-Order Reversible Binding kinetics }
-\description{
- Function describing the solution of the differential equations describing
- the kinetic model with first-order terms for a two-way transfer from a free
- to a bound fraction, and a first-order degradation term for the free
- fraction. The initial condition is a defined amount in the free fraction and
- no substance in the bound fraction.
-}
-\usage{
- SFORB.solution(t, parent.0, k_12, k_21, k_1output)
-}
-\arguments{
- \item{t}{ Time. }
- \item{parent.0}{ Starting value for the response variable at time zero. }
- \item{k_12}{ Kinetic constant describing transfer from free to bound. }
- \item{k_21}{ Kinetic constant describing transfer from bound to free. }
- \item{k_1output}{ Kinetic constant describing degradation of the free fraction. }
-}
-\value{
- The value of the response variable, which is the sum of free and bound
- fractions at time \code{t}.
-}
-\references{
- FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence and
- Degradation Kinetics from Environmental Fate Studies on Pesticides in EU
- Registration} Report of the FOCUS Work Group on Degradation Kinetics,
- EC Document Reference Sanco/10058/2005 version 2.0, 434 pp,
- \url{http://focus.jrc.ec.europa.eu/dk}
-}
-\examples{
- \dontrun{plot(function(x) SFORB.solution(x, 100, 0.5, 2, 3), 0, 2)}
-}
-\keyword{ manip }
diff --git a/man/endpoints.Rd b/man/endpoints.Rd deleted file mode 100644 index 21316cf..0000000 --- a/man/endpoints.Rd +++ /dev/null @@ -1,33 +0,0 @@ -\name{endpoints} -\alias{endpoints} -\title{ -Function to calculate endpoints for further use from kinetic models fitted with mkinfit -} -\description{ -This function calculates DT50 and DT90 values as well as formation fractions from kinetic models -fitted with mkinfit. If the SFORB model was specified for one of the parents or metabolites, -the Eigenvalues are returned. These are equivalent to the rate constantes of the DFOP model, but -with the advantage that the SFORB model can also be used for metabolites. -} -\usage{ -endpoints(fit, pseudoDT50 = FALSE) -} -\arguments{ - \item{fit}{ - An object of class \code{\link{mkinfit}}. -} - \item{pseudoDT50}{ - Should pseudoDT50 values for FOMC, DFOP and SFORB models be reported, as - recommended by the FOCUS group? Currently not implemented. -} -} -\note{ - The function is used internally by \code{\link{summary.mkinfit}}. -} -\value{ - A list with the components mentioned above. -} -\author{ - Johannes Ranke -} -\keyword{ manip } diff --git a/man/ilr.Rd b/man/ilr.Rd deleted file mode 100644 index cedb49c..0000000 --- a/man/ilr.Rd +++ /dev/null @@ -1,56 +0,0 @@ -\name{ilr} -\alias{ilr} -\alias{invilr} -\title{ - Function to perform isotropic log-ratio transformation -} -\description{ - This implementation is a special case of the class of isotropic log-ratio transformations. -} -\usage{ - ilr(x) - invilr(x) -} -\arguments{ - \item{x}{ - A numeric vector. Naturally, the forward transformation is only sensible for - vectors with all elements being greater than zero. - } -} -\value{ - The result of the forward or backward transformation. The returned components always - sum to 1 for the case of the inverse log-ratio transformation. -} -\references{ - Peter Filzmoser, Karel Hron (2008) Outlier Detection for Compositional Data Using Robust Methods. Math Geosci 40 233-248 -} -\author{ - René Lehmann and Johannes Ranke -} -\seealso{ - Other implementations are in R packages \code{compositions} and \code{robCompositions}. -} -\examples{ -# Order matters -ilr(c(0.1, 1, 10)) -ilr(c(10, 1, 0.1)) -# Equal entries give ilr transformations with zeros as elements -ilr(c(3, 3, 3)) -# Almost equal entries give small numbers -ilr(c(0.3, 0.4, 0.3)) -# Only the ratio between the numbers counts, not their sum -invilr(ilr(c(0.7, 0.29, 0.01))) -invilr(ilr(2.1 * c(0.7, 0.29, 0.01))) -# Inverse transformation of larger numbers gives unequal elements -invilr(-10) -invilr(c(-10, 0)) -# The sum of the elements of the inverse ilr is 1 -sum(invilr(c(-10, 0))) -# This is why we do not need all elements of the inverse transformation to go back: -a <- c(0.1, 0.3, 0.5) -b <- invilr(a) -length(b) # Four elements -ilr(c(b[1:3], 1 - sum(b[1:3]))) # Gives c(0.1, 0.3, 0.5) -} - -\keyword{ manip } diff --git a/man/mccall81_245T.Rd b/man/mccall81_245T.Rd deleted file mode 100644 index 6c1a137..0000000 --- a/man/mccall81_245T.Rd +++ /dev/null @@ -1,42 +0,0 @@ -\name{mccall81_245T}
-\alias{mccall81_245T}
-\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}
-\format{
- A dataframe 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}
- \item{\code{soil}}{a factor containing the name of the soil}
- }
-}
-\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{
- SFO_SFO_SFO <- mkinmod(T245 = list(type = "SFO", to = "phenol"),
- phenol = list(type = "SFO", to = "anisole"),
- anisole = list(type = "SFO"))
- \dontrun{fit.1 <- mkinfit(SFO_SFO_SFO, subset(mccall81_245T, soil == "Commerce"))
- summary(fit.1, data = FALSE)}
- # No covariance matrix and k_phenol_sink is really small, therefore fix it to zero
- fit.2 <- mkinfit(SFO_SFO_SFO, subset(mccall81_245T, soil == "Commerce"),
- parms.ini = c(k_phenol_sink = 0),
- fixed_parms = "k_phenol_sink")
- summary(fit.2, data = FALSE)
-}
-\keyword{datasets}
diff --git a/man/mkin_long_to_wide.Rd b/man/mkin_long_to_wide.Rd deleted file mode 100644 index e583664..0000000 --- a/man/mkin_long_to_wide.Rd +++ /dev/null @@ -1,36 +0,0 @@ -\name{mkin_long_to_wide} -\alias{mkin_long_to_wide} -\title{ - Convert a dataframe from long to wide format. -} -\usage{ -mkin_long_to_wide(long_data, time = "time", outtime = "time") -} -\description{ - This function takes a dataframe in the long form as required by \code{\link{modCost}} - and converts it into a dataframe with one independent variable and several - dependent variables as columns. -} -\arguments{ - \item{long_data}{ - The dataframe must contain one variable called "time" with the time values specified by the - \code{time} argument, one column called "name" with the grouping of the observed values, and - finally one column of observed values called "value". -} - \item{time}{ - The name of the time variable in the long input data. -} - \item{outtime}{ - The name of the time variable in the wide output data. -} -} -\value{ - Dataframe in wide format. -} -\author{ - Johannes Ranke -} -\examples{ -mkin_long_to_wide(FOCUS_2006_D) -} -\keyword{ manip } diff --git a/man/mkin_wide_to_long.Rd b/man/mkin_wide_to_long.Rd deleted file mode 100644 index d3dd200..0000000 --- a/man/mkin_wide_to_long.Rd +++ /dev/null @@ -1,32 +0,0 @@ -\name{mkin_wide_to_long} -\alias{mkin_wide_to_long} -\title{ - Convert a dataframe with observations over time into long format. -} -\usage{ -mkin_wide_to_long(wide_data, time = "t") -} -\description{ - This function simply takes a dataframe with one independent variable and several - dependent variable and converts it into the long form as required by \code{\link{modCost}}. -} -\arguments{ - \item{wide_data}{ - The dataframe must contain one variable with the time values specified by the - \code{time} argument and usually more than one column of observed values. -} - \item{time}{ - The name of the time variable. -} -} -\value{ - Dataframe in long format as needed for \code{\link{modCost}}. -} -\author{ - Johannes Ranke -} -\examples{ -wide <- data.frame(t = c(1,2,3), x = c(1,4,7), y = c(3,4,5)) -mkin_wide_to_long(wide) -} -\keyword{ manip } diff --git a/man/mkinerrmin.Rd b/man/mkinerrmin.Rd deleted file mode 100644 index c43d87a..0000000 --- a/man/mkinerrmin.Rd +++ /dev/null @@ -1,44 +0,0 @@ -\name{mkinerrmin}
-\Rdversion{1.1}
-\alias{mkinerrmin}
-\title{
-Calculate the minimum error to assume in order to pass the variance test
-}
-\description{
-This function uses \code{\link{optimize}} in order to iteratively find the
-smallest relative error still resulting in passing the chi-squared test
-as defined in the FOCUS kinetics report from 2006.
-}
-\usage{
-mkinerrmin(fit, alpha = 0.05)
-}
-\arguments{
- \item{fit}{
- an object of class \code{\link{mkinfit}}.
- }
- \item{alpha}{
- The confidence level chosen for the chi-squared test.
-}
-}
-\value{
- A dataframe with the following components:
- \item{err.min}{The relative error, expressed as a fraction.}
- \item{n.optim}{The number of optimised parameters attributed to the data series.}
- \item{df}{The number of remaining degrees of freedom for the chi2 error level
- calculations. Note that mean values are used for the chi2 statistic and
- therefore every time point with observed values in the series only counts
- one time.}
- The dataframe has one row for the total dataset and one further row for
- each observed state variable in the model.
-}
-\details{
- This function is used internally by \code{\link{summary.mkinfit}}.
-}
-\references{
- FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence and
- Degradation Kinetics from Environmental Fate Studies on Pesticides in EU
- Registration} Report of the FOCUS Work Group on Degradation Kinetics,
- EC Document Reference Sanco/10058/2005 version 2.0, 434 pp,
- \url{http://focus.jrc.ec.europa.eu/dk}
-}
-\keyword{ manip }
diff --git a/man/mkinfit.Rd b/man/mkinfit.Rd deleted file mode 100644 index 823ceee..0000000 --- a/man/mkinfit.Rd +++ /dev/null @@ -1,254 +0,0 @@ -\name{mkinfit} -\alias{mkinfit} -\title{ - Fit a kinetic model to data with one or more state variables. -} -\description{ - This function uses the Flexible Modelling Environment package - \code{\link{FME}} to create a function calculating the model cost, i.e. the - deviation between the kinetic model and the observed data. This model cost is - then minimised using the Levenberg-Marquardt algorithm \code{\link{nls.lm}}, - using the specified initial or fixed parameters and starting values. - In each step of the optimsation, the kinetic model is solved using the - function \code{\link{mkinpredict}}. The variance of the residuals for each - observed variable can optionally be iteratively reweighted until convergence - using the argument \code{reweight.method = "obs"}. -} -\usage{ -mkinfit(mkinmod, observed, - parms.ini = "auto", - state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)), - fixed_parms = NULL, fixed_initials = names(mkinmod$diffs)[-1], - solution_type = "auto", - method.ode = "lsoda", - method.modFit = "Marq", - control.modFit = list(), - transform_rates = TRUE, - transform_fractions = TRUE, - plot = FALSE, quiet = FALSE, err = NULL, weight = "none", - scaleVar = FALSE, - atol = 1e-8, rtol = 1e-10, n.outtimes = 100, - reweight.method = NULL, - reweight.tol = 1e-8, reweight.max.iter = 10, - trace_parms = FALSE, ...) -} -\arguments{ - \item{mkinmod}{ - A list of class \code{\link{mkinmod}}, containing the kinetic model to be fitted to the data. - } - \item{observed}{ - The observed data. It has to be in the long format as described in - \code{\link{modFit}}, i.e. the first column called "name" must contain the - name of the observed variable for each data point. The second column must - contain the times of observation, named "time". The third column must be - named "value" and contain the observed values. Optionally, a further column - can contain weights for each data point. If it is not named "err", its name - must be passed as a further argument named \code{err} which is then passed - on to \code{\link{modFit}}. - } - \item{parms.ini}{ - A named vector of initial values for the parameters, including parameters - to be optimised and potentially also fixed parameters as indicated by - \code{fixed_parms}. If set to "auto", initial values for rate constants - are set to default values. Using parameter names that are not in the model - gives an error. - - It is possible to only specify a subset of the parameters that the model - needs. You can use the parameter lists "bparms.ode" from a previously - fitted model, which contains the differential equation parameters from this - model. This works nicely if the models are nested. An example is given - below. - } - \item{state.ini}{ - A named vector of initial values for the state variables of the model. In - case the observed variables are represented by more than one model - variable, the names will differ from the names of the observed variables - (see \code{map} component of \code{\link{mkinmod}}). The default is to set - the initial value of the first model variable to 100 and all others to 0. - } - \item{fixed_parms}{ - The names of parameters that should not be optimised but rather kept at the - values specified in \code{parms.ini}. - } - \item{fixed_initials}{ - The names of model variables for which the initial state at time 0 should - be excluded from the optimisation. Defaults to all state variables except - for the first one. - } - \item{solution_type}{ - If set to "eigen", the solution of the system of differential equations is - based on the spectral decomposition of the coefficient matrix in cases that - this is possible. If set to "deSolve", a numerical ode solver from package - \code{\link{deSolve}} is used. If set to "analytical", an analytical - solution of the model is used. This is only implemented for simple - degradation experiments with only one state variable, i.e. with no - metabolites. The default is "auto", which uses "analytical" if possible, - otherwise "eigen" if the model can be expressed using eigenvalues and - eigenvectors, and finally "deSolve" for the remaining models (time - dependence of degradation rates and metabolites). This argument is passed - on to the helper function \code{\link{mkinpredict}}. - } - \item{method.ode}{ - The solution method passed via \code{\link{mkinpredict}} to - \code{\link{ode}} in case the solution type is "deSolve". The default - "lsoda" is performant, but sometimes fails to converge. - } - \item{method.modFit}{ - The optimisation method passed to \code{\link{modFit}}. The default "Marq" - is the Levenberg Marquardt algorithm \code{\link{nls.lm}} from the package - \code{minpack.lm}. Often other methods need more iterations to find the - same result. When using "Pseudo", "upper" and "lower" need to be - specified for the transformed parameters. - } - \item{control.modFit}{ - Additional arguments passed to the optimisation method used by - \code{\link{modFit}}. - } - \item{transform_rates}{ - Boolean specifying if kinetic rate constants should be transformed in the - model specification used in the fitting for better compliance with the - 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 - model specification used in the fitting for better compliance with the - assumption of normal distribution of the estimator. The default (TRUE) is - to do transformations. The g parameter of the DFOP and HS models are also - transformed, as they can also be seen as compositional data. The - transformation used for these transformations is the \code{\link{ilr}} - transformation. - } - \item{plot}{ - Should the observed values and the numerical solutions be plotted at each - stage of the optimisation? - } - \item{quiet}{ - Suppress printing out the current model cost after each improvement? - } - \item{err }{either \code{NULL}, or the name of the column with the - \emph{error} estimates, used to weigh the residuals (see details of - \code{\link{modCost}}); if \code{NULL}, then the residuals are not weighed. - } - \item{weight}{ - only if \code{err}=\code{NULL}: how to weight the residuals, one of "none", - "std", "mean", see details of \code{\link{modCost}}. - } - \item{scaleVar}{ - Will be passed to \code{\link{modCost}}. Default is not to scale Variables - according to the number of observations. - } - \item{atol}{ - Absolute error tolerance, passed to \code{\link{ode}}. Default is 1e-8, - lower than in \code{\link{lsoda}}. - } - \item{rtol}{ - Absolute error tolerance, passed to \code{\link{ode}}. Default is 1e-10, - much lower than in \code{\link{lsoda}}. - } - \item{n.outtimes}{ - The length of the dataseries that is produced by the model prediction - function \code{\link{mkinpredict}}. This impacts the accuracy of - the numerical solver if that is used (see \code{solution_type} argument. - The default value is 100. - } - \item{reweight.method}{ - The method used for iteratively reweighting residuals, also known - as iteratively reweighted least squares (IRLS). Default is NULL, - the other method implemented is called "obs", meaning that each - observed variable is assumed to have its own variance, this is - estimated from the fit and used for weighting the residuals - in each iteration until convergence of this estimate up to - \code{reweight.tol} or up to the maximum number of iterations - specified by \code{reweight.max.iter}. - } - \item{reweight.tol}{ - Tolerance for convergence criterion for the variance components - in IRLS fits. - } - \item{reweight.max.iter}{ - Maximum iterations in IRLS fits. - } - \item{trace_parms}{ - Should a trace of the parameter values be listed? - } - \item{\dots}{ - Further arguments that will be passed to \code{\link{modFit}}. - } -} -\value{ - A list with "mkinfit" and "modFit" in the class attribute. - A summary can be obtained by \code{\link{summary.mkinfit}}. -} -\note{ - The implementation of iteratively reweighted least squares is inspired by the - work of the KinGUII team at Bayer Crop Science (Walter Schmitt and Zhenglei - Gao). A similar implemention can also be found in CAKE 2.0, which is the - other GUI derivative of mkin, sponsored by Syngenta. -} -\author{ - Johannes Ranke <jranke@uni-bremen.de> -} -\examples{ -# One parent compound, one metabolite, both single first order. -SFO_SFO <- mkinmod( - parent = list(type = "SFO", to = "m1", sink = TRUE), - m1 = list(type = "SFO")) -# Fit the model to the FOCUS example dataset D using defaults -fit <- mkinfit(SFO_SFO, FOCUS_2006_D) -summary(fit) - -# Use stepwise fitting, using optimised parameters from parent only fit, FOMC -\dontrun{ -FOMC <- mkinmod(parent = list(type = "FOMC")) -FOMC_SFO <- mkinmod( - parent = list(type = "FOMC", to = "m1", sink = TRUE), - m1 = list(type = "SFO")) -# Fit the model to the FOCUS example dataset D using defaults -fit.FOMC_SFO <- mkinfit(FOMC_SFO, FOCUS_2006_D) -# Use starting parameters from parent only FOMC fit -fit.FOMC = mkinfit(FOMC, FOCUS_2006_D, plot=TRUE) -fit.FOMC_SFO <- mkinfit(FOMC_SFO, FOCUS_2006_D, - parms.ini = fit.FOMC$bparms.ode, plot=TRUE) - -# Use stepwise fitting, using optimised parameters from parent only fit, SFORB -SFORB <- mkinmod(parent = list(type = "SFORB")) -SFORB_SFO <- mkinmod( - parent = list(type = "SFORB", to = "m1", sink = TRUE), - m1 = list(type = "SFO")) -# Fit the model to the FOCUS example dataset D using defaults -fit.SFORB_SFO <- mkinfit(SFORB_SFO, FOCUS_2006_D) -# Use starting parameters from parent only SFORB fit (not really needed in this case) -fit.SFORB = mkinfit(SFORB, FOCUS_2006_D) -fit.SFORB_SFO <- mkinfit(SFORB_SFO, FOCUS_2006_D, parms.ini = fit.SFORB$bparms.ode) -} - -# Weighted fits, including IRLS -SFO_SFO.ff <- mkinmod(parent = list(type = "SFO", to = "m1"), - m1 = list(type = "SFO"), use_of_ff = "max") -f.noweight <- mkinfit(SFO_SFO.ff, FOCUS_2006_D) -summary(f.noweight) -f.irls <- mkinfit(SFO_SFO.ff, FOCUS_2006_D, reweight.method = "obs") -summary(f.irls) -f.w.mean <- mkinfit(SFO_SFO.ff, FOCUS_2006_D, weight = "mean") -summary(f.w.mean) -f.w.mean.irls <- mkinfit(SFO_SFO.ff, FOCUS_2006_D, weight = "mean", - reweight.method = "obs") -summary(f.w.mean.irls) - -\dontrun{ -# Manual weighting -dw <- FOCUS_2006_D -errors <- c(parent = 2, m1 = 1) -dw$err.man <- errors[FOCUS_2006_D$name] -f.w.man <- mkinfit(SFO_SFO.ff, dw, err = "err.man") -summary(f.w.man) -f.w.man.irls <- mkinfit(SFO_SFO.ff, dw, err = "err.man", - reweight.method = "obs") -summary(f.w.man.irls) -} -} -\keyword{ models } -\keyword{ optimize } diff --git a/man/mkinmod.Rd b/man/mkinmod.Rd deleted file mode 100644 index 76127c5..0000000 --- a/man/mkinmod.Rd +++ /dev/null @@ -1,63 +0,0 @@ -\name{mkinmod}
-\alias{mkinmod}
-\title{
- Function to set up a kinetic model with one or more state variables.
-}
-\description{
- The function takes a specification, consisting of a list of the observed variables
- in the data. Each observed variable is again represented by a list, specifying the
- kinetic model type and reaction or transfer to other observed compartments.
-}
-\usage{
-mkinmod(..., use_of_ff = "min", speclist = NULL)
-}
-\arguments{
- \item{...}{
- For each observed variable, a list has to be specified as an argument, containing
- at least a component \code{type}, specifying the type of kinetics to use
- for the variable. Currently, single first order kinetics "SFO" or
- single first order with reversible binding "SFORB" are implemented for all
- variables, while
- "FOMC", "DFOP" and "HS" can additionally be chosen for the first
- variable which is assumed to be the source compartment.
- Additionally, each component of the list can include a character vector \code{to},
- specifying names of variables to which a transfer is to be assumed in the
- model.
- If the argument \code{use_of_ff} is set to "min" (default) and the model for
- the compartment is "SFO" or "SFORB", an additional component of the list
- can be "sink=FALSE" effectively fixing the flux to sink to zero.
- }
- \item{use_of_ff}{
- Specification of the use of formation fractions in the model equations and, if
- applicable, the coefficient matrix. If "min", a minimum use of formation
- fractions is made in order to avoid fitting the product of formation fractions
- and rate constants. If "max", formation fractions are always used.
- }
- \item{speclist}{
- The specification of the observed variables and their submodel types and
- pathways can be given as a single list using this argument. Default is NULL.
- }
-}
-\value{
- A list of class \code{mkinmod} for use with \code{\link{mkinfit}}, containing
- \item{diffs}{ A vector of string representations of differential equations,
- one for each modelling variable. }
- \item{parms}{ A vector of parameter names occurring in the differential equations. }
- \item{map}{ A list containing named character vectors for each observed variable, specifying
- the modelling variables by which it is represented. }
- \item{use_of_ff}{ The content of \code{use_of_ff} is passed on in this list component. }
- \item{coefmat}{ The coefficient matrix, if the system of differential equations can be represented by one. }
-}
-\author{
- Johannes Ranke
-}
-\examples{
-# Specify the SFO model
-SFO <- mkinmod(parent = list(type = "SFO"))
-
-# One parent compound, one metabolite, both single first order.
-SFO_SFO <- mkinmod(
- parent = list(type = "SFO", to = "m1"),
- m1 = list(type = "SFO"))
-}
-\keyword{ models }
diff --git a/man/mkinparplot.Rd b/man/mkinparplot.Rd deleted file mode 100644 index 60103b8..0000000 --- a/man/mkinparplot.Rd +++ /dev/null @@ -1,35 +0,0 @@ -\name{mkinparplot}
-\alias{mkinparplot}
-\title{
- Function to plot the confidence intervals obtained using
- \code{\link{mkinfit}}
-}
-\description{
- This function plots the confidence intervals for the parameters
- fitted using \code{\link{mkinfit}}.
-}
-\usage{
- mkinparplot(object)
-}
-\arguments{
- \item{object}{
- A fit represented in an \code{\link{mkinfit}} object.
- }
-}
-\value{
- Nothing is returned by this function, as it is called for its side effect, namely to produce a plot.
-}
-\author{
- Johannes Ranke
-}
-
-\examples{
-model <- mkinmod(
- T245 = list(type = "SFO", to = c("phenol")),
- phenol = list(type = "SFO", to = c("anisole")),
- anisole = list(type = "SFO"), use_of_ff = "max")
-fit <- mkinfit(model, subset(mccall81_245T, soil == "Commerce"),
- parms.ini = c(f_phenol_to_anisole = 1), fixed_parms = "f_phenol_to_anisole")
-\dontrun{mkinparplot(fit)}
-}
-\keyword{ hplot }
diff --git a/man/mkinplot.Rd b/man/mkinplot.Rd deleted file mode 100644 index 4b0fef4..0000000 --- a/man/mkinplot.Rd +++ /dev/null @@ -1,26 +0,0 @@ -\name{mkinplot} -\alias{mkinplot} -\title{ - Plot the observed data and the fitted model of an mkinfit. -} -\description{ - Deprecated function. It now only calls the plot method \code{\link{plot.mkinfit}}. -} -\usage{ - mkinplot(fit, ...) -} -\arguments{ - \item{fit}{ - an object of class \code{\link{mkinfit}}. - } - \item{\dots}{ - further arguments passed to \code{\link{plot.mkinfit}}. -} -} -\value{ - The function is called for its side effect. -} -\author{ - Johannes Ranke -} -\keyword{ hplot } diff --git a/man/mkinpredict.Rd b/man/mkinpredict.Rd deleted file mode 100644 index 97db90e..0000000 --- a/man/mkinpredict.Rd +++ /dev/null @@ -1,91 +0,0 @@ -\name{mkinpredict} -\alias{mkinpredict} -\title{ - Produce predictions from a kinetic model using specifc parameters -} -\description{ - This function produces a time series for all the observed variables in a - kinetic model as specified by \code{\link{mkinmod}}, using a specific set of - kinetic parameters and initial values for the state variables. -} -\usage{ - mkinpredict(mkinmod, odeparms, odeini, outtimes, solution_type = "deSolve", - method.ode = "lsoda", atol = 1e-08, rtol = 1e-10, map_output = TRUE, ...) -} -\arguments{ - \item{mkinmod}{ - A kinetic model as produced by \code{\link{mkinmod}}. - } - \item{odeparms}{ - A numeric vector specifying the parameters used in the kinetic model, which - is generally defined as a set of ordinary differential equations. - } - \item{odeini}{ - A numeric vectory containing the initial values of the state variables of - the model. Note that the state variables can differ from the observed - variables, for example in the case of the SFORB model. - } - \item{outtimes}{ - A numeric vector specifying the time points for which model predictions - should be generated. - } - \item{solution_type}{ - The method that should be used for producing the predictions. This should - generally be "analytical" if there is only one observed variable, and - usually "deSolve" in the case of several observed variables. The third - possibility "eigen" is faster but not applicable to some models e.g. - using FOMC for the parent compound. - } - \item{method.ode}{ - The solution method passed via \code{\link{mkinpredict}} to - \code{\link{ode}} in case the solution type is "deSolve". The default - "lsoda" is performant, but sometimes fails to converge. - } - \item{atol}{ - Absolute error tolerance, passed to \code{\link{ode}}. Default is 1e-8, - lower than in \code{\link{lsoda}}. - } - \item{rtol}{ - Absolute error tolerance, passed to \code{\link{ode}}. Default is 1e-10, - much lower than in \code{\link{lsoda}}. - } - \item{map_output}{ - Boolean to specify if the output should list values for the observed - variables (default) or for all state variables (if set to FALSE). - } - \item{\dots}{ - Further arguments passed to the ode solver in case such a solver is used. - } -} -\value{ - A matrix in the same format as the output of \code{\link{ode}}. -} -\author{ - Johannes Ranke -} -\examples{ - SFO <- mkinmod(degradinol = list(type = "SFO")) - mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, - solution_type = "analytical") - mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, - solution_type = "eigen") - - mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 1:20, - solution_type = "analytical")[20,] - mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, - atol = 1e-20)[20,] - - # The integration method does not make a lot of difference - mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, - atol = 1e-20, method = "ode45")[20,] - mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, - atol = 1e-20, method = "rk4")[20,] - - # The number of output times used to make a lot of difference until the - # default for atol was adjusted - mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), - seq(0, 20, by = 0.1))[201,] - mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), - seq(0, 20, by = 0.01))[2001,] -} -\keyword{ manip } diff --git a/man/mkinresplot.Rd b/man/mkinresplot.Rd deleted file mode 100644 index 3f53dd1..0000000 --- a/man/mkinresplot.Rd +++ /dev/null @@ -1,67 +0,0 @@ -\name{mkinresplot}
-\alias{mkinresplot}
-\title{
- Function to plot residuals stored in an mkin object
-}
-\description{
- This function plots the residuals for the specified subset of the
- observed variables from an mkinfit object. A combined plot of the fitted
- model and the residuals can be obtained using \code{\link{plot.mkinfit}}
- using the argument \code{show_residuals = TRUE}.
-}
-\usage{
- mkinresplot(object,
- obs_vars = names(object$mkinmod$map),
- xlab = "Time", ylab = "Residual",
- maxabs = "auto", legend = TRUE, lpos = "topright", ...)
-}
-\arguments{
- \item{object}{
- A fit represented in an \code{\link{mkinfit}} object.
-}
- \item{obs_vars}{
- A character vector of names of the observed variables for which residuals
- should be plotted. Defaults to all observed variables in the model
-}
- \item{xlab}{
- Label for the x axis. Defaults to "Time [days]".
-}
- \item{ylab}{
- Label for the y axis. Defaults to "Residual [\% of applied radioactivity]".
-}
- \item{maxabs}{
- Maximum absolute value of the residuals. This is used for the scaling of
- the y axis and defaults to "auto".
-}
- \item{legend}{
- Should a legend be plotted? Defaults to "TRUE".
-}
- \item{lpos}{
- Where should the legend be placed? Default is "topright". Will be passed on to
- \code{\link{legend}}. }
- \item{\dots}{
- further arguments passed to \code{\link{plot}}.
-}
-}
-\value{
- Nothing is returned by this function, as it is called for its side effect, namely to produce a plot.
-}
-\author{
- Johannes Ranke
-}
-
-\seealso{
- \code{\link{mkinplot}}, for a way to plot the data and the fitted lines of the
- mkinfit object. }
-\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{fit <- mkinfit(model, data, plot=TRUE)}
-\dontrun{mkinresplot(fit, "A1")}
-}
-\keyword{ hplot }
diff --git a/man/plot.mkinfit.Rd b/man/plot.mkinfit.Rd deleted file mode 100644 index 7009e7d..0000000 --- a/man/plot.mkinfit.Rd +++ /dev/null @@ -1,94 +0,0 @@ -\name{plot.mkinfit} -\alias{plot.mkinfit} -\title{ - Plot the observed data and the fitted model of an mkinfit object. -} -\description{ - Solves the differential equations with the optimised and fixed parameters - from a previous successful call to \code{\link{mkinfit}} and plots - the observed data together with the solution of the fitted model. -} -\usage{ -\method{plot}{mkinfit}(x, fit = x, - obs_vars = names(fit$mkinmod$map), - xlab = "Time", ylab = "Observed", - xlim = range(fit$data$time), - ylim = "default", - col_obs = 1:length(fit$mkinmod$map), pch_obs = col_obs, - lty_obs = rep(1, length(fit$mkinmod$map)), - add = FALSE, legend = !add, - show_residuals = FALSE, maxabs = "auto", - lpos = "topright", inset = c(0.05, 0.05), \dots) -} -\arguments{ - \item{x}{ - Alias for fit introduced for compatibility with the generic S3 method. - } - \item{fit}{ - an object of class \code{\link{mkinfit}}. - } - \item{obs_vars}{ - A character vector of names of the observed variables for which the - data and the model should be plotted. Defauls to all observed variables - in the model. - } - \item{xlab}{ - label for the x axis. - } - \item{ylab}{ - label for the y axis. - } - \item{xlim}{ - plot range in x direction. - } - \item{ylim}{ - plot range in y direction. - } - \item{col_obs}{ - colors used for plotting the observed data and the corresponding model prediction lines. - } - \item{pch_obs}{ - symbols to be used for plotting the data. - } - \item{lty_obs}{ - line types to be used for the model predictions. - } - \item{add}{ - should the plot be added to an existing plot? - } - \item{legend}{ - should a legend be added to the plot? - } - \item{show_residuals}{ - should residuals be shown in the lower third of the plot? - } - \item{maxabs}{ - Maximum absolute value of the residuals. This is used for the scaling of - the y axis and defaults to "auto". - } - \item{lpos}{ - position of the legend. Passed to \code{\link{legend}} as the first argument. - } - \item{inset}{ - Passed to \code{\link{legend}} if applicable. - } - \item{\dots}{ - further arguments passed to \code{\link{plot}}. - } -} -\value{ - The function is called for its side effect. -} -\examples{ -# One parent compound, one metabolite, both single first order. -SFO_SFO <- mkinmod( - parent = list(type = "SFO", to = "m1", sink = TRUE), - m1 = list(type = "SFO")) -# Fit the model to the FOCUS example dataset D using defaults -fit <- mkinfit(SFO_SFO, FOCUS_2006_D) -\dontrun{plot(fit)} -} -\author{ - Johannes Ranke -} -\keyword{ hplot } diff --git a/man/schaefer07_complex_case.Rd b/man/schaefer07_complex_case.Rd deleted file mode 100644 index ed4f694..0000000 --- a/man/schaefer07_complex_case.Rd +++ /dev/null @@ -1,42 +0,0 @@ -\name{schaefer07_complex_case}
-\alias{schaefer07_complex_case}
-\alias{schaefer07_complex_results}
-\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).
- The results from the fitting are also included.
-}
-\usage{data(schaefer07_complex_case)}
-\format{
- The data set is 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}
- }
- The results are a data frame with 14 results for different parameter values
-}
-\references{
- 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}
diff --git a/man/summary.mkinfit.Rd b/man/summary.mkinfit.Rd deleted file mode 100644 index 472b5de..0000000 --- a/man/summary.mkinfit.Rd +++ /dev/null @@ -1,75 +0,0 @@ -\name{summary.mkinfit}
-\alias{summary.mkinfit}
-\alias{print.summary.mkinfit}
-\title{
- Summary method for class "mkinfit".
-}
-\description{
- Lists model equations, the summary as returned by \code{\link{summary.modFit}},
- the chi2 error levels calculated according to FOCUS guidance (2006) as far
- as defined therein, and optionally the data, consisting of observed, predicted
- and residual values.
-}
-\usage{
-\method{summary}{mkinfit}(object, data = TRUE, distimes = TRUE, alpha = 0.05, ...)
-\method{print}{summary.mkinfit}(x, digits = max(3, getOption("digits") - 3), ...)
-}
-
-\arguments{
- \item{object}{
- an object of class \code{\link{mkinfit}}.
-}
- \item{x}{
- an object of class \code{summary.mkinfit}.
-}
- \item{data}{
- logical, indicating whether the data should be included in the summary.
-}
- \item{distimes}{
- logical, indicating whether DT50 and DT90 values should be included.
-}
- \item{alpha}{
- error level for confidence interval estimation from t distribution
-}
- \item{digits}{
- Number of digits to use for printing
-}
- \item{\dots}{
- optional arguments passed to methods like \code{print}.
-}
-}
-\value{
- The summary function returns a list derived from
- \code{\link{summary.modFit}}, with components, among others
- \item{version, Rversion}{The mkin and R versions used}
- \item{date.fit, date.summary}{The dates where the fit and the summary were produced}
- \item{use_of_ff}{Was maximum or minimum use made of formation fractions}
- \item{residuals, residualVariance, sigma, modVariance, df}{As in summary.modFit}
- \item{cov.unscaled, cov.scaled, info, niter, stopmess, par}{As in summary.modFit}
- \item{bpar}{Optimised and backtransformed parameters}
- \item{diffs }{The differential equations used in the model}
- \item{data }{The data (see Description above).}
- \item{start }{The starting values and bounds, if applicable, for optimised parameters.}
- \item{fixed }{The values of fixed parameters.}
- \item{errmin }{The chi2 error levels for each observed variable.}
- \item{bparms.ode }{All backtransformed ODE parameters, for use as starting parameters for
- related models.}
- \item{ff }{The estimated formation fractions derived from the fitted model.}
- \item{distimes }{The DT50 and DT90 values for each observed variable.}
- \item{SFORB}{If applicable, eigenvalues of SFORB components of the model.}
- The print method is called for its side effect, i.e. printing the summary.
-}
-\references{
- FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence and
- Degradation Kinetics from Environmental Fate Studies on Pesticides in EU
- Registration} Report of the FOCUS Work Group on Degradation Kinetics,
- EC Document Reference Sanco/10058/2005 version 2.0, 434 pp,
- \url{http://focus.jrc.ec.europa.eu/dk}
-}
-\author{
- Johannes Ranke <jranke@{harlan.com,uni-bremen.de}>
-}
-\examples{
- summary(mkinfit(mkinmod(parent = list(type = "SFO")), FOCUS_2006_A))
-}
-\keyword{ utilities }
diff --git a/man/transform_odeparms.Rd b/man/transform_odeparms.Rd deleted file mode 100644 index c52fb4f..0000000 --- a/man/transform_odeparms.Rd +++ /dev/null @@ -1,69 +0,0 @@ -\name{transform_odeparms} -\alias{transform_odeparms} -\alias{backtransform_odeparms} -\title{ - Functions to transform and backtransform kinetic parameters for fitting -} -\description{ - The transformations are intended to map parameters that should only take - on restricted values to the full scale of real numbers. For kinetic rate - constants and other paramters that can only take on positive values, a - simple log transformation is used. For compositional parameters, such as - the formations fractions that should always sum up to 1 and can not be - negative, the \code{\link{ilr}} transformation is used. -} -\usage{ -transform_odeparms(parms, mod_vars, transform_rates = TRUE, transform_fractions = TRUE) -backtransform_odeparms(transparms, mod_vars, - transform_rates = TRUE, transform_fractions = TRUE) -} -\arguments{ - \item{parms}{ - Parameters of kinetic models as used in the differential equations. - } - \item{transparms}{ - Transformed parameters of kinetic models as used in the fitting procedure. - } - \item{mod_vars}{ - Names of the state variables in the kinetic model. These are used for - grouping the formation fractions before \code{\link{ilr}} transformation. - } - \item{transform_rates}{ - Boolean specifying if kinetic rate constants should be transformed in the - model specification used in the fitting for better compliance with the - 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 and the break point tb - of the HS model - } - \item{transform_fractions}{ - Boolean specifying if formation fractions constants should be transformed in the - model specification used in the fitting for better compliance with the - assumption of normal distribution of the estimator. The default (TRUE) is - to do transformations. The g parameter of the DFOP and HS models are also - transformed, as they can also be seen as compositional data. The - transformation used for these transformations is the \code{\link{ilr}} - transformation. - } -} -\value{ - A vector of transformed or backtransformed parameters with the same names - as the original parameters. -} -\author{ - Johannes Ranke -} -\examples{ -SFO_SFO <- mkinmod( - parent = list(type = "SFO", to = "m1", sink = TRUE), - m1 = list(type = "SFO")) -# Fit the model to the FOCUS example dataset D using defaults -fit <- mkinfit(SFO_SFO, FOCUS_2006_D) -summary(fit, data=FALSE) # See transformed and backtransformed parameters -initials <- fit$start$value -transformed <- fit$start$transformed -names(initials) <- names(transformed) <- rownames(fit$start) -transform_odeparms(initials, c("parent", "m1")) -backtransform_odeparms(transformed, c("parent", "m1")) -} -\keyword{ manip } diff --git a/tests/doRUnit.R b/tests/doRUnit.R deleted file mode 100644 index f0f8281..0000000 --- a/tests/doRUnit.R +++ /dev/null @@ -1,43 +0,0 @@ -# $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)) { - - pkg <- "mkin" - path <- system.file(package=pkg, "unitTests") - - cat("\nRunning unit tests\n") - print(list(pkg=pkg, getwd=getwd(), pathToUnitTests=path)) - - library(package=pkg, character.only=TRUE) - - ## Define tests - testSuite <- defineTestSuite(name=paste(pkg, " Unit Tests"), - dirs=path) - ## Run - tests <- runTestSuite(testSuite) - - ## Default report name - pathReport <- file.path(path, "report") - - ## Report to stdout and text files - cat("------------------- UNIT TEST SUMMARY ---------------------\n\n") - printTextProtocol(tests, showDetails=FALSE) - printTextProtocol(tests, showDetails=FALSE, - fileName=paste(pathReport, "Summary.txt", sep="")) - printTextProtocol(tests, showDetails=TRUE, - fileName=paste(pathReport, ".txt", sep="")) - - ## Report to HTML file - printHTMLProtocol(tests, fileName=paste(pathReport, ".html", sep="")) - - ## Return stop() to cause R CMD check stop in case of - ## - failures i.e. FALSE to unit tests or - ## - errors i.e. R errors - tmp <- getErrors(tests) - if(tmp$nFail > 0 | tmp$nErr > 0) { - stop(paste("\n\nunit testing failed (#test failures: ", tmp$nFail, - ", #R errors: ", tmp$nErr, ")\n\n", sep="")) - } -} else { - warning("cannot run unit tests -- package RUnit is not available") -} diff --git a/vignettes/FOCUS_L.Rmd b/vignettes/FOCUS_L.Rmd deleted file mode 100644 index 957b34a..0000000 --- a/vignettes/FOCUS_L.Rmd +++ /dev/null @@ -1,243 +0,0 @@ -<!--
-%\VignetteEngine{knitr::knitr}
-%\VignetteIndexEntry{Example evaluation of FOCUS Laboratory Data L1 to L3}
--->
-
-# Example evaluation of FOCUS Laboratory Data L1 to L3
-
-## Laboratory Data L1
-
-The following code defines example dataset L1 from the FOCUS kinetics
-report, p. 284
-
-```{r}
-library("mkin")
-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,
- 27.7, 27.3, 10.0, 10.4, 2.9, 4.0))
-FOCUS_2006_L1_mkin <- mkin_wide_to_long(FOCUS_2006_L1)
-```
-
-The next step is to set up the models used for the kinetic analysis. Note that
-the model definitions contain the names of the observed variables in the data.
-In this case, there is only one variable called `parent`.
-
-```{r}
-SFO <- mkinmod(parent = list(type = "SFO"))
-FOMC <- mkinmod(parent = list(type = "FOMC"))
-DFOP <- mkinmod(parent = list(type = "DFOP"))
-```
-
-The three models cover the first assumption of simple first order (SFO),
-the case of declining rate constant over time (FOMC) and the case of two
-different phases of the kinetics (DFOP). For a more detailed discussion
-of the models, please see the FOCUS kinetics report.
-
-The following two lines fit the model and produce the summary report
-of the model fit. This covers the numerical analysis given in the
-FOCUS report.
-
-```{r}
-m.L1.SFO <- mkinfit(SFO, FOCUS_2006_L1_mkin, quiet=TRUE)
-summary(m.L1.SFO)
-```
-
-A plot of the fit is obtained with the plot function for mkinfit objects.
-
-```{r fig.width=7, fig.height = 5}
-plot(m.L1.SFO)
-```
-The residual plot can be easily obtained by
-
-```{r fig.width=7, fig.height = 5}
-mkinresplot(m.L1.SFO, ylab = "Observed", xlab = "Time")
-```
-
-For comparison, the FOMC model is fitted as well, and the chi^2 error level
-is checked.
-
-```{r}
-m.L1.FOMC <- mkinfit(FOMC, FOCUS_2006_L1_mkin, quiet=TRUE)
-summary(m.L1.FOMC, data = FALSE)
-```
-
-Due to the higher number of parameters, and the lower number of degrees of
-freedom of the fit, the chi^2 error level is actually higher for the FOMC
-model (3.6%) than for the SFO model (3.4%). Additionally, the covariance
-matrix can not be obtained, indicating overparameterisation of the model.
-As a consequence, no standard errors for transformed parameters nor
-confidence intervals for backtransformed parameters are available.
-
-The chi^2 error levels reported in Appendix 3 and Appendix 7 to the FOCUS
-kinetics report are rounded to integer percentages and partly deviate by one
-percentage point from the results calculated by mkin. The reason for
-this is not known. However, mkin gives the same chi^2 error levels
-as the kinfit package.
-
-Furthermore, the calculation routines of the kinfit package have been extensively
-compared to the results obtained by the KinGUI software, as documented in the
-kinfit package vignette. KinGUI is a widely used standard package in this field.
-Therefore, the reason for the difference was not investigated further.
-
-## Laboratory Data L2
-
-The following code defines example dataset L2 from the FOCUS kinetics
-report, p. 287
-
-```{r}
-FOCUS_2006_L2 = data.frame(
- t = rep(c(0, 1, 3, 7, 14, 28), each = 2),
- parent = c(96.1, 91.8, 41.4, 38.7,
- 19.3, 22.3, 4.6, 4.6,
- 2.6, 1.2, 0.3, 0.6))
-FOCUS_2006_L2_mkin <- mkin_wide_to_long(FOCUS_2006_L2)
-```
-
-Again, the SFO model is fitted and a summary is obtained.
-
-```{r}
-m.L2.SFO <- mkinfit(SFO, FOCUS_2006_L2_mkin, quiet=TRUE)
-summary(m.L2.SFO)
-```
-
-The chi^2 error level of 14% suggests that the model does not fit very well.
-This is also obvious from the plots of the fit and the residuals.
-
-```{r fig.height = 8}
-par(mfrow = c(2, 1))
-plot(m.L2.SFO)
-mkinresplot(m.L2.SFO)
-```
-
-In the FOCUS kinetics report, it is stated that there is no apparent systematic
-error observed from the residual plot up to the measured DT90 (approximately at
-day 5), and there is an underestimation beyond that point.
-
-We may add that it is difficult to judge the random nature of the residuals just
-from the three samplings at days 0, 1 and 3. Also, it is not clear _a
-priori_ why a consistent underestimation after the approximate DT90 should be
-irrelevant. However, this can be rationalised by the fact that the FOCUS fate
-models generally only implement SFO kinetics.
-
-For comparison, the FOMC model is fitted as well, and the chi^2 error level
-is checked.
-
-```{r fig.height = 8}
-m.L2.FOMC <- mkinfit(FOMC, FOCUS_2006_L2_mkin, quiet = TRUE)
-par(mfrow = c(2, 1))
-plot(m.L2.FOMC)
-mkinresplot(m.L2.FOMC)
-summary(m.L2.FOMC, data = FALSE)
-```
-
-The error level at which the chi^2 test passes is much lower in this case.
-Therefore, the FOMC model provides a better description of the data, as less
-experimental error has to be assumed in order to explain the data.
-
-Fitting the four parameter DFOP model further reduces the chi^2 error level.
-
-```{r fig.height = 5}
-m.L2.DFOP <- mkinfit(DFOP, FOCUS_2006_L2_mkin, quiet = TRUE)
-plot(m.L2.DFOP)
-```
-
-Here, the default starting parameters for the DFOP model obviously do not lead
-to a reasonable solution. Therefore the fit is repeated with different starting
-parameters.
-
-```{r fig.height = 5}
-m.L2.DFOP <- mkinfit(DFOP, FOCUS_2006_L2_mkin,
- parms.ini = c(k1 = 1, k2 = 0.01, g = 0.8),
- quiet=TRUE)
-plot(m.L2.DFOP)
-summary(m.L2.DFOP, data = FALSE)
-```
-
-Here, the DFOP model is clearly the best-fit model for dataset L2 based on the
-chi^2 error level criterion. However, the failure to calculate the covariance
-matrix indicates that the parameter estimates correlate excessively. Therefore,
-the FOMC model may be preferred for this dataset.
-
-## Laboratory Data L3
-
-The following code defines example dataset L3 from the FOCUS kinetics report,
-p. 290.
-
-```{r}
-FOCUS_2006_L3 = data.frame(
- t = c(0, 3, 7, 14, 30, 60, 91, 120),
- parent = c(97.8, 60, 51, 43, 35, 22, 15, 12))
-FOCUS_2006_L3_mkin <- mkin_wide_to_long(FOCUS_2006_L3)
-```
-
-SFO model, summary and plot:
-
-```{r fig.height = 5}
-m.L3.SFO <- mkinfit(SFO, FOCUS_2006_L3_mkin, quiet = TRUE)
-plot(m.L3.SFO)
-summary(m.L3.SFO)
-```
-
-The chi^2 error level of 21% as well as the plot suggest that the model
-does not fit very well.
-
-The FOMC model performs better:
-
-```{r fig.height = 5}
-m.L3.FOMC <- mkinfit(FOMC, FOCUS_2006_L3_mkin, quiet = TRUE)
-plot(m.L3.FOMC)
-summary(m.L3.FOMC, data = FALSE)
-```
-
-The error level at which the chi^2 test passes is 7% in this case.
-
-Fitting the four parameter DFOP model further reduces the chi^2 error level
-considerably:
-
-```{r fig.height = 5}
-m.L3.DFOP <- mkinfit(DFOP, FOCUS_2006_L3_mkin, quiet = TRUE)
-plot(m.L3.DFOP)
-summary(m.L3.DFOP, data = FALSE)
-```
-
-Here, a look to the model plot, the confidence intervals of the parameters
-and the correlation matrix suggest that the parameter estimates are reliable, and
-the DFOP model can be used as the best-fit model based on the chi^2 error
-level criterion for laboratory data L3.
-
-## Laboratory Data L4
-
-The following code defines example dataset L4 from the FOCUS kinetics
-report, p. 293
-
-```{r}
-FOCUS_2006_L4 = data.frame(
- t = c(0, 3, 7, 14, 30, 60, 91, 120),
- parent = c(96.6, 96.3, 94.3, 88.8, 74.9, 59.9, 53.5, 49.0))
-FOCUS_2006_L4_mkin <- mkin_wide_to_long(FOCUS_2006_L4)
-```
-
-SFO model, summary and plot:
-
-```{r fig.height = 5}
-m.L4.SFO <- mkinfit(SFO, FOCUS_2006_L4_mkin, quiet = TRUE)
-plot(m.L4.SFO)
-summary(m.L4.SFO, data = FALSE)
-```
-
-The chi^2 error level of 3.3% as well as the plot suggest that the model
-fits very well.
-
-The FOMC model for comparison
-
-```{r fig.height = 5}
-m.L4.FOMC <- mkinfit(FOMC, FOCUS_2006_L4_mkin, quiet = TRUE)
-plot(m.L4.FOMC)
-summary(m.L4.FOMC, data = FALSE)
-```
-
-The error level at which the chi^2 test passes is slightly lower for the FOMC
-model. However, the difference appears negligible.
-
diff --git a/vignettes/FOCUS_L.md b/vignettes/FOCUS_L.md deleted file mode 100644 index 6c43889..0000000 --- a/vignettes/FOCUS_L.md +++ /dev/null @@ -1,931 +0,0 @@ -<!-- -%\VignetteEngine{knitr::knitr} -%\VignetteIndexEntry{Example evaluation of FOCUS Laboratory Data L1 to L3} ---> - -# Example evaluation of FOCUS Laboratory Data L1 to L3 - -## Laboratory Data L1 - -The following code defines example dataset L1 from the FOCUS kinetics -report, p. 284 - - -```r -library("mkin") -``` - -``` -## Loading required package: FME -## Loading required package: deSolve -## Loading required package: rootSolve -## Loading required package: minpack.lm -## Loading required package: MASS -## Loading required package: coda -## Loading required package: lattice -``` - -```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, 71.9, 50.3, 59.4, 47, - 45.1, 27.7, 27.3, 10, 10.4, 2.9, 4)) -FOCUS_2006_L1_mkin <- mkin_wide_to_long(FOCUS_2006_L1) -``` - - -The next step is to set up the models used for the kinetic analysis. Note that -the model definitions contain the names of the observed variables in the data. -In this case, there is only one variable called `parent`. - - -```r -SFO <- mkinmod(parent = list(type = "SFO")) -FOMC <- mkinmod(parent = list(type = "FOMC")) -DFOP <- mkinmod(parent = list(type = "DFOP")) -``` - - -The three models cover the first assumption of simple first order (SFO), -the case of declining rate constant over time (FOMC) and the case of two -different phases of the kinetics (DFOP). For a more detailed discussion -of the models, please see the FOCUS kinetics report. - -The following two lines fit the model and produce the summary report -of the model fit. This covers the numerical analysis given in the -FOCUS report. - - -```r -m.L1.SFO <- mkinfit(SFO, FOCUS_2006_L1_mkin, quiet = TRUE) -summary(m.L1.SFO) -``` - -``` -## mkin version: 0.9.25 -## R version: 3.0.2 -## Date of fit: Sun Nov 17 15:02:54 2013 -## Date of summary: Sun Nov 17 15:02:54 2013 -## -## Equations: -## [1] d_parent = - k_parent_sink * parent -## -## Method used for solution of differential equation system: -## analytical -## -## Weighting: none -## -## Starting values for optimised parameters: -## value type transformed -## parent_0 100.0 state 100.000 -## k_parent_sink 0.1 deparm -2.303 -## -## Fixed parameter values: -## None -## -## Optimised, transformed parameters: -## Estimate Std. Error Lower Upper -## parent_0 92.50 1.3700 89.60 95.40 -## k_parent_sink -2.35 0.0406 -2.43 -2.26 -## -## Backtransformed parameters: -## Estimate Lower Upper -## parent_0 92.5000 89.6000 95.400 -## k_parent_sink 0.0956 0.0877 0.104 -## -## Residual standard error: 2.95 on 16 degrees of freedom -## -## Chi2 error levels in percent: -## err.min n.optim df -## All data 3.42 2 7 -## parent 3.42 2 7 -## -## Estimated disappearance times: -## DT50 DT90 -## parent 7.25 24.1 -## -## Estimated formation fractions: -## ff -## parent_sink 1 -## -## Parameter correlation: -## parent_0 k_parent_sink -## parent_0 1.000 0.625 -## k_parent_sink 0.625 1.000 -## -## Data: -## time variable observed predicted residual -## 0 parent 88.3 92.47 -4.171 -## 0 parent 91.4 92.47 -1.071 -## 1 parent 85.6 84.04 1.561 -## 1 parent 84.5 84.04 0.461 -## 2 parent 78.9 76.38 2.524 -## 2 parent 77.6 76.38 1.224 -## 3 parent 72.0 69.41 2.588 -## 3 parent 71.9 69.41 2.488 -## 5 parent 50.3 57.33 -7.030 -## 5 parent 59.4 57.33 2.070 -## 7 parent 47.0 47.35 -0.352 -## 7 parent 45.1 47.35 -2.252 -## 14 parent 27.7 24.25 3.453 -## 14 parent 27.3 24.25 3.053 -## 21 parent 10.0 12.42 -2.416 -## 21 parent 10.4 12.42 -2.016 -## 30 parent 2.9 5.25 -2.351 -## 30 parent 4.0 5.25 -1.251 -``` - - -A plot of the fit is obtained with the plot function for mkinfit objects. - - -```r -plot(m.L1.SFO) -``` - - - -The residual plot can be easily obtained by - - -```r -mkinresplot(m.L1.SFO, ylab = "Observed", xlab = "Time") -``` - - - - -For comparison, the FOMC model is fitted as well, and the chi^2 error level -is checked. - - -```r -m.L1.FOMC <- mkinfit(FOMC, FOCUS_2006_L1_mkin, quiet = TRUE) -summary(m.L1.FOMC, data = FALSE) -``` - -``` -## mkin version: 0.9.25 -## R version: 3.0.2 -## Date of fit: Sun Nov 17 15:02:55 2013 -## Date of summary: Sun Nov 17 15:02:55 2013 -## -## Equations: -## [1] d_parent = - (alpha/beta) * ((time/beta) + 1)^-1 * parent -## -## Method used for solution of differential equation system: -## analytical -## -## Weighting: none -## -## Starting values for optimised parameters: -## value type transformed -## parent_0 100 state 100.000 -## alpha 1 deparm 0.000 -## beta 10 deparm 2.303 -## -## Fixed parameter values: -## None -## -## Optimised, transformed parameters: -## Estimate Std. Error Lower Upper -## parent_0 92.5 NA NA NA -## alpha 25.6 NA NA NA -## beta 28.0 NA NA NA -## -## Backtransformed parameters: -## Estimate Lower Upper -## parent_0 9.25e+01 NA NA -## alpha 1.35e+11 NA NA -## beta 1.41e+12 NA NA -## -## Residual standard error: 3.05 on 15 degrees of freedom -## -## Chi2 error levels in percent: -## err.min n.optim df -## All data 3.62 3 6 -## parent 3.62 3 6 -## -## Estimated disappearance times: -## DT50 DT90 -## parent 7.25 24.1 -## -## Estimated formation fractions: -## ff -## parent_sink 1 -## -## Parameter correlation: -## Could not estimate covariance matrix; singular system: -``` - - -Due to the higher number of parameters, and the lower number of degrees of -freedom of the fit, the chi^2 error level is actually higher for the FOMC -model (3.6%) than for the SFO model (3.4%). Additionally, the covariance -matrix can not be obtained, indicating overparameterisation of the model. -As a consequence, no standard errors for transformed parameters nor -confidence intervals for backtransformed parameters are available. - -The chi^2 error levels reported in Appendix 3 and Appendix 7 to the FOCUS -kinetics report are rounded to integer percentages and partly deviate by one -percentage point from the results calculated by mkin. The reason for -this is not known. However, mkin gives the same chi^2 error levels -as the kinfit package. - -Furthermore, the calculation routines of the kinfit package have been extensively -compared to the results obtained by the KinGUI software, as documented in the -kinfit package vignette. KinGUI is a widely used standard package in this field. -Therefore, the reason for the difference was not investigated further. - -## Laboratory Data L2 - -The following code defines example dataset L2 from the FOCUS kinetics -report, p. 287 - - -```r -FOCUS_2006_L2 = data.frame(t = rep(c(0, 1, 3, 7, 14, 28), each = 2), parent = c(96.1, - 91.8, 41.4, 38.7, 19.3, 22.3, 4.6, 4.6, 2.6, 1.2, 0.3, 0.6)) -FOCUS_2006_L2_mkin <- mkin_wide_to_long(FOCUS_2006_L2) -``` - - -Again, the SFO model is fitted and a summary is obtained. - - -```r -m.L2.SFO <- mkinfit(SFO, FOCUS_2006_L2_mkin, quiet = TRUE) -summary(m.L2.SFO) -``` - -``` -## mkin version: 0.9.25 -## R version: 3.0.2 -## Date of fit: Sun Nov 17 15:02:55 2013 -## Date of summary: Sun Nov 17 15:02:55 2013 -## -## Equations: -## [1] d_parent = - k_parent_sink * parent -## -## Method used for solution of differential equation system: -## analytical -## -## Weighting: none -## -## Starting values for optimised parameters: -## value type transformed -## parent_0 100.0 state 100.000 -## k_parent_sink 0.1 deparm -2.303 -## -## Fixed parameter values: -## None -## -## Optimised, transformed parameters: -## Estimate Std. Error Lower Upper -## parent_0 91.500 3.810 83.000 99.900 -## k_parent_sink -0.411 0.107 -0.651 -0.172 -## -## Backtransformed parameters: -## Estimate Lower Upper -## parent_0 91.500 83.000 99.900 -## k_parent_sink 0.663 0.522 0.842 -## -## Residual standard error: 5.51 on 10 degrees of freedom -## -## Chi2 error levels in percent: -## err.min n.optim df -## All data 14.4 2 4 -## parent 14.4 2 4 -## -## Estimated disappearance times: -## DT50 DT90 -## parent 1.05 3.47 -## -## Estimated formation fractions: -## ff -## parent_sink 1 -## -## Parameter correlation: -## parent_0 k_parent_sink -## parent_0 1.00 0.43 -## k_parent_sink 0.43 1.00 -## -## Data: -## time variable observed predicted residual -## 0 parent 96.1 9.15e+01 4.634 -## 0 parent 91.8 9.15e+01 0.334 -## 1 parent 41.4 4.71e+01 -5.740 -## 1 parent 38.7 4.71e+01 -8.440 -## 3 parent 19.3 1.25e+01 6.779 -## 3 parent 22.3 1.25e+01 9.779 -## 7 parent 4.6 8.83e-01 3.717 -## 7 parent 4.6 8.83e-01 3.717 -## 14 parent 2.6 8.53e-03 2.591 -## 14 parent 1.2 8.53e-03 1.191 -## 28 parent 0.3 7.96e-07 0.300 -## 28 parent 0.6 7.96e-07 0.600 -``` - - -The chi^2 error level of 14% suggests that the model does not fit very well. -This is also obvious from the plots of the fit and the residuals. - - -```r -par(mfrow = c(2, 1)) -plot(m.L2.SFO) -mkinresplot(m.L2.SFO) -``` - - - - -In the FOCUS kinetics report, it is stated that there is no apparent systematic -error observed from the residual plot up to the measured DT90 (approximately at -day 5), and there is an underestimation beyond that point. - -We may add that it is difficult to judge the random nature of the residuals just -from the three samplings at days 0, 1 and 3. Also, it is not clear _a -priori_ why a consistent underestimation after the approximate DT90 should be -irrelevant. However, this can be rationalised by the fact that the FOCUS fate -models generally only implement SFO kinetics. - -For comparison, the FOMC model is fitted as well, and the chi^2 error level -is checked. - - -```r -m.L2.FOMC <- mkinfit(FOMC, FOCUS_2006_L2_mkin, quiet = TRUE) -par(mfrow = c(2, 1)) -plot(m.L2.FOMC) -mkinresplot(m.L2.FOMC) -``` - - - -```r -summary(m.L2.FOMC, data = FALSE) -``` - -``` -## mkin version: 0.9.25 -## R version: 3.0.2 -## Date of fit: Sun Nov 17 15:02:56 2013 -## Date of summary: Sun Nov 17 15:02:56 2013 -## -## Equations: -## [1] d_parent = - (alpha/beta) * ((time/beta) + 1)^-1 * parent -## -## Method used for solution of differential equation system: -## analytical -## -## Weighting: none -## -## Starting values for optimised parameters: -## value type transformed -## parent_0 100 state 100.000 -## alpha 1 deparm 0.000 -## beta 10 deparm 2.303 -## -## Fixed parameter values: -## None -## -## Optimised, transformed parameters: -## Estimate Std. Error Lower Upper -## parent_0 93.800 1.860 89.600 98.000 -## alpha 0.318 0.187 -0.104 0.740 -## beta 0.210 0.294 -0.456 0.876 -## -## Backtransformed parameters: -## Estimate Lower Upper -## parent_0 93.80 89.600 98.0 -## alpha 1.37 0.901 2.1 -## beta 1.23 0.634 2.4 -## -## Residual standard error: 2.63 on 9 degrees of freedom -## -## Chi2 error levels in percent: -## err.min n.optim df -## All data 6.2 3 3 -## parent 6.2 3 3 -## -## Estimated disappearance times: -## DT50 DT90 -## parent 0.809 5.36 -## -## Estimated formation fractions: -## ff -## parent_sink 1 -## -## Parameter correlation: -## parent_0 alpha beta -## parent_0 1.0000 -0.0955 -0.186 -## alpha -0.0955 1.0000 0.976 -## beta -0.1863 0.9757 1.000 -``` - - -The error level at which the chi^2 test passes is much lower in this case. -Therefore, the FOMC model provides a better description of the data, as less -experimental error has to be assumed in order to explain the data. - -Fitting the four parameter DFOP model further reduces the chi^2 error level. - - -```r -m.L2.DFOP <- mkinfit(DFOP, FOCUS_2006_L2_mkin, quiet = TRUE) -plot(m.L2.DFOP) -``` - - - - -Here, the default starting parameters for the DFOP model obviously do not lead -to a reasonable solution. Therefore the fit is repeated with different starting -parameters. - - -```r -m.L2.DFOP <- mkinfit(DFOP, FOCUS_2006_L2_mkin, parms.ini = c(k1 = 1, k2 = 0.01, - g = 0.8), quiet = TRUE) -plot(m.L2.DFOP) -``` - - - -```r -summary(m.L2.DFOP, data = FALSE) -``` - -``` -## mkin version: 0.9.25 -## R version: 3.0.2 -## Date of fit: Sun Nov 17 15:02:57 2013 -## Date of summary: Sun Nov 17 15:02:57 2013 -## -## 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 -## -## Method used for solution of differential equation system: -## analytical -## -## Weighting: none -## -## Starting values for optimised parameters: -## value type transformed -## parent_0 1e+02 state 100.0000 -## k1 1e+00 deparm 0.0000 -## k2 1e-02 deparm -4.6052 -## g 8e-01 deparm 0.9803 -## -## Fixed parameter values: -## None -## -## Optimised, transformed parameters: -## Estimate Std. Error Lower Upper -## parent_0 93.900 NA NA NA -## k1 4.960 NA NA NA -## k2 -1.090 NA NA NA -## g -0.282 NA NA NA -## -## Backtransformed parameters: -## Estimate Lower Upper -## parent_0 93.900 NA NA -## k1 142.000 NA NA -## k2 0.337 NA NA -## g 0.402 NA NA -## -## Residual standard error: 1.73 on 8 degrees of freedom -## -## Chi2 error levels in percent: -## err.min n.optim df -## All data 2.53 4 2 -## parent 2.53 4 2 -## -## Estimated disappearance times: -## DT50 DT90 -## parent NA NA -## -## Estimated formation fractions: -## ff -## parent_sink 1 -## -## Parameter correlation: -## Could not estimate covariance matrix; singular system: -``` - - -Here, the DFOP model is clearly the best-fit model for dataset L2 based on the -chi^2 error level criterion. However, the failure to calculate the covariance -matrix indicates that the parameter estimates correlate excessively. Therefore, -the FOMC model may be preferred for this dataset. - -## Laboratory Data L3 - -The following code defines example dataset L3 from the FOCUS kinetics report, -p. 290. - - -```r -FOCUS_2006_L3 = data.frame(t = c(0, 3, 7, 14, 30, 60, 91, 120), parent = c(97.8, - 60, 51, 43, 35, 22, 15, 12)) -FOCUS_2006_L3_mkin <- mkin_wide_to_long(FOCUS_2006_L3) -``` - - -SFO model, summary and plot: - - -```r -m.L3.SFO <- mkinfit(SFO, FOCUS_2006_L3_mkin, quiet = TRUE) -plot(m.L3.SFO) -``` - - - -```r -summary(m.L3.SFO) -``` - -``` -## mkin version: 0.9.25 -## R version: 3.0.2 -## Date of fit: Sun Nov 17 15:02:57 2013 -## Date of summary: Sun Nov 17 15:02:57 2013 -## -## Equations: -## [1] d_parent = - k_parent_sink * parent -## -## Method used for solution of differential equation system: -## analytical -## -## Weighting: none -## -## Starting values for optimised parameters: -## value type transformed -## parent_0 100.0 state 100.000 -## k_parent_sink 0.1 deparm -2.303 -## -## Fixed parameter values: -## None -## -## Optimised, transformed parameters: -## Estimate Std. Error Lower Upper -## parent_0 74.90 8.460 54.20 95.60 -## k_parent_sink -3.68 0.326 -4.48 -2.88 -## -## Backtransformed parameters: -## Estimate Lower Upper -## parent_0 74.9000 54.2000 95.6000 -## k_parent_sink 0.0253 0.0114 0.0561 -## -## Residual standard error: 12.9 on 6 degrees of freedom -## -## Chi2 error levels in percent: -## err.min n.optim df -## All data 21.2 2 6 -## parent 21.2 2 6 -## -## Estimated disappearance times: -## DT50 DT90 -## parent 27.4 91.1 -## -## Estimated formation fractions: -## ff -## parent_sink 1 -## -## Parameter correlation: -## parent_0 k_parent_sink -## parent_0 1.000 0.548 -## k_parent_sink 0.548 1.000 -## -## Data: -## time variable observed predicted residual -## 0 parent 97.8 74.87 22.9273 -## 3 parent 60.0 69.41 -9.4065 -## 7 parent 51.0 62.73 -11.7340 -## 14 parent 43.0 52.56 -9.5634 -## 30 parent 35.0 35.08 -0.0828 -## 60 parent 22.0 16.44 5.5614 -## 91 parent 15.0 7.51 7.4896 -## 120 parent 12.0 3.61 8.3908 -``` - - -The chi^2 error level of 21% as well as the plot suggest that the model -does not fit very well. - -The FOMC model performs better: - - -```r -m.L3.FOMC <- mkinfit(FOMC, FOCUS_2006_L3_mkin, quiet = TRUE) -plot(m.L3.FOMC) -``` - - - -```r -summary(m.L3.FOMC, data = FALSE) -``` - -``` -## mkin version: 0.9.25 -## R version: 3.0.2 -## Date of fit: Sun Nov 17 15:02:57 2013 -## Date of summary: Sun Nov 17 15:02:57 2013 -## -## Equations: -## [1] d_parent = - (alpha/beta) * ((time/beta) + 1)^-1 * parent -## -## Method used for solution of differential equation system: -## analytical -## -## Weighting: none -## -## Starting values for optimised parameters: -## value type transformed -## parent_0 100 state 100.000 -## alpha 1 deparm 0.000 -## beta 10 deparm 2.303 -## -## Fixed parameter values: -## None -## -## Optimised, transformed parameters: -## Estimate Std. Error Lower Upper -## parent_0 97.000 4.550 85.3 109.000 -## alpha -0.862 0.170 -1.3 -0.424 -## beta 0.619 0.474 -0.6 1.840 -## -## Backtransformed parameters: -## Estimate Lower Upper -## parent_0 97.000 85.300 109.000 -## alpha 0.422 0.273 0.655 -## beta 1.860 0.549 6.290 -## -## Residual standard error: 4.57 on 5 degrees of freedom -## -## Chi2 error levels in percent: -## err.min n.optim df -## All data 7.32 3 5 -## parent 7.32 3 5 -## -## Estimated disappearance times: -## DT50 DT90 -## parent 7.73 431 -## -## Estimated formation fractions: -## ff -## parent_sink 1 -## -## Parameter correlation: -## parent_0 alpha beta -## parent_0 1.000 -0.151 -0.427 -## alpha -0.151 1.000 0.911 -## beta -0.427 0.911 1.000 -``` - - -The error level at which the chi^2 test passes is 7% in this case. - -Fitting the four parameter DFOP model further reduces the chi^2 error level -considerably: - - -```r -m.L3.DFOP <- mkinfit(DFOP, FOCUS_2006_L3_mkin, quiet = TRUE) -plot(m.L3.DFOP) -``` - - - -```r -summary(m.L3.DFOP, data = FALSE) -``` - -``` -## mkin version: 0.9.25 -## R version: 3.0.2 -## Date of fit: Sun Nov 17 15:02:58 2013 -## Date of summary: Sun Nov 17 15:02:58 2013 -## -## 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 -## -## Method used for solution of differential equation system: -## analytical -## -## Weighting: none -## -## Starting values for optimised parameters: -## value type transformed -## parent_0 1e+02 state 100.000 -## k1 1e-01 deparm -2.303 -## k2 1e-02 deparm -4.605 -## g 5e-01 deparm 0.000 -## -## Fixed parameter values: -## None -## -## Optimised, transformed parameters: -## Estimate Std. Error Lower Upper -## parent_0 97.700 1.4400 93.800 102.0000 -## k1 -0.661 0.1330 -1.030 -0.2910 -## k2 -4.290 0.0590 -4.450 -4.1200 -## g -0.123 0.0512 -0.265 0.0193 -## -## Backtransformed parameters: -## Estimate Lower Upper -## parent_0 97.7000 93.8000 102.0000 -## k1 0.5160 0.3560 0.7480 -## k2 0.0138 0.0117 0.0162 -## g 0.4570 0.4070 0.5070 -## -## Residual standard error: 1.44 on 4 degrees of freedom -## -## Chi2 error levels in percent: -## err.min n.optim df -## All data 2.23 4 4 -## parent 2.23 4 4 -## -## Estimated disappearance times: -## DT50 DT90 -## parent 7.46 123 -## -## Estimated formation fractions: -## ff -## parent_sink 1 -## -## Parameter correlation: -## parent_0 k1 k2 g -## parent_0 1.0000 0.164 0.0131 0.425 -## k1 0.1640 1.000 0.4648 -0.553 -## k2 0.0131 0.465 1.0000 -0.663 -## g 0.4253 -0.553 -0.6631 1.000 -``` - - -Here, a look to the model plot, the confidence intervals of the parameters -and the correlation matrix suggest that the parameter estimates are reliable, and -the DFOP model can be used as the best-fit model based on the chi^2 error -level criterion for laboratory data L3. - -## Laboratory Data L4 - -The following code defines example dataset L4 from the FOCUS kinetics -report, p. 293 - - -```r -FOCUS_2006_L4 = data.frame(t = c(0, 3, 7, 14, 30, 60, 91, 120), parent = c(96.6, - 96.3, 94.3, 88.8, 74.9, 59.9, 53.5, 49)) -FOCUS_2006_L4_mkin <- mkin_wide_to_long(FOCUS_2006_L4) -``` - - -SFO model, summary and plot: - - -```r -m.L4.SFO <- mkinfit(SFO, FOCUS_2006_L4_mkin, quiet = TRUE) -plot(m.L4.SFO) -``` - - - -```r -summary(m.L4.SFO, data = FALSE) -``` - -``` -## mkin version: 0.9.25 -## R version: 3.0.2 -## Date of fit: Sun Nov 17 15:02:58 2013 -## Date of summary: Sun Nov 17 15:02:58 2013 -## -## Equations: -## [1] d_parent = - k_parent_sink * parent -## -## Method used for solution of differential equation system: -## analytical -## -## Weighting: none -## -## Starting values for optimised parameters: -## value type transformed -## parent_0 100.0 state 100.000 -## k_parent_sink 0.1 deparm -2.303 -## -## Fixed parameter values: -## None -## -## Optimised, transformed parameters: -## Estimate Std. Error Lower Upper -## parent_0 96.40 1.95 91.70 101.00 -## k_parent_sink -5.03 0.08 -5.23 -4.83 -## -## Backtransformed parameters: -## Estimate Lower Upper -## parent_0 96.40000 91.70000 1.01e+02 -## k_parent_sink 0.00654 0.00538 7.95e-03 -## -## Residual standard error: 3.65 on 6 degrees of freedom -## -## Chi2 error levels in percent: -## err.min n.optim df -## All data 3.29 2 6 -## parent 3.29 2 6 -## -## Estimated disappearance times: -## DT50 DT90 -## parent 106 352 -## -## Estimated formation fractions: -## ff -## parent_sink 1 -## -## Parameter correlation: -## parent_0 k_parent_sink -## parent_0 1.000 0.587 -## k_parent_sink 0.587 1.000 -``` - - -The chi^2 error level of 3.3% as well as the plot suggest that the model -fits very well. - -The FOMC model for comparison - - -```r -m.L4.FOMC <- mkinfit(FOMC, FOCUS_2006_L4_mkin, quiet = TRUE) -plot(m.L4.FOMC) -``` - - - -```r -summary(m.L4.FOMC, data = FALSE) -``` - -``` -## mkin version: 0.9.25 -## R version: 3.0.2 -## Date of fit: Sun Nov 17 15:02:59 2013 -## Date of summary: Sun Nov 17 15:02:59 2013 -## -## Equations: -## [1] d_parent = - (alpha/beta) * ((time/beta) + 1)^-1 * parent -## -## Method used for solution of differential equation system: -## analytical -## -## Weighting: none -## -## Starting values for optimised parameters: -## value type transformed -## parent_0 100 state 100.000 -## alpha 1 deparm 0.000 -## beta 10 deparm 2.303 -## -## Fixed parameter values: -## None -## -## Optimised, transformed parameters: -## Estimate Std. Error Lower Upper -## parent_0 99.100 1.680 94.80 103.000 -## alpha -0.351 0.372 -1.31 0.607 -## beta 4.170 0.564 2.73 5.620 -## -## Backtransformed parameters: -## Estimate Lower Upper -## parent_0 99.100 94.80 103.00 -## alpha 0.704 0.27 1.83 -## beta 65.000 15.30 277.00 -## -## Residual standard error: 2.31 on 5 degrees of freedom -## -## Chi2 error levels in percent: -## err.min n.optim df -## All data 2.03 3 5 -## parent 2.03 3 5 -## -## Estimated disappearance times: -## DT50 DT90 -## parent 109 1644 -## -## Estimated formation fractions: -## ff -## parent_sink 1 -## -## Parameter correlation: -## parent_0 alpha beta -## parent_0 1.000 -0.536 -0.608 -## alpha -0.536 1.000 0.991 -## beta -0.608 0.991 1.000 -``` - - -The error level at which the chi^2 test passes is slightly lower for the FOMC -model. However, the difference appears negligible. - diff --git a/vignettes/FOCUS_Z.Rnw b/vignettes/FOCUS_Z.Rnw deleted file mode 100644 index b551007..0000000 --- a/vignettes/FOCUS_Z.Rnw +++ /dev/null @@ -1,310 +0,0 @@ -%\VignetteIndexEntry{Examples evaluation of FOCUS dataset Z} -%\VignetteEngine{knitr::knitr} -\documentclass[12pt,a4paper]{article} -\usepackage{a4wide} -\input{header} -\hypersetup{ - pdftitle = {Example evaluation of FOCUS dataset Z}, - pdfsubject = {Manuscript}, - pdfauthor = {Johannes Ranke}, - colorlinks = {true}, - linkcolor = {blue}, - citecolor = {blue}, - urlcolor = {red}, - hyperindex = {true}, - linktocpage = {true}, -} - -\begin{document} - -<<include=FALSE>>= -require(knitr) -opts_chunk$set(engine='R', tidy=FALSE) -@ - -\title{Example evaluation of FOCUS dataset Z} -\author{\textbf{Johannes Ranke} \\[0.5cm] -%EndAName -Wissenschaftlicher Berater\\ -Kronacher Str. 8, 79639 Grenzach-Wyhlen, Germany\\[0.5cm] -and\\[0.5cm] -University of Bremen\\ -} -\maketitle - -\thispagestyle{empty} \setcounter{page}{0} - -\clearpage - -\tableofcontents - -\textbf{Key words}: Kinetics, FOCUS, nonlinear optimisation - -\section{The data} - -The following code defines the example dataset from Appendix 7 to the FOCUS kinetics -report \citep{FOCUSkinetics2011}, p.350. - -<<FOCUS_2006_Z_data, echo=TRUE, eval=TRUE>>= -require(mkin) -LOD = 0.5 -FOCUS_2006_Z = data.frame( - t = c(0, 0.04, 0.125, 0.29, 0.54, 1, 2, 3, 4, 7, 10, 14, 21, - 42, 61, 96, 124), - Z0 = c(100, 81.7, 70.4, 51.1, 41.2, 6.6, 4.6, 3.9, 4.6, 4.3, 6.8, - 2.9, 3.5, 5.3, 4.4, 1.2, 0.7), - Z1 = c(0, 18.3, 29.6, 46.3, 55.1, 65.7, 39.1, 36, 15.3, 5.6, 1.1, - 1.6, 0.6, 0.5 * LOD, NA, NA, NA), - Z2 = c(0, NA, 0.5 * LOD, 2.6, 3.8, 15.3, 37.2, 31.7, 35.6, 14.5, - 0.8, 2.1, 1.9, 0.5 * LOD, NA, NA, NA), - Z3 = c(0, NA, NA, NA, NA, 0.5 * LOD, 9.2, 13.1, 22.3, 28.4, 32.5, - 25.2, 17.2, 4.8, 4.5, 2.8, 4.4)) - -FOCUS_2006_Z_mkin <- mkin_wide_to_long(FOCUS_2006_Z) -@ - -\section{Parent compound and one metabolite} - -The next step is to set up the models used for the kinetic analysis. As the -simultaneous fit of parent and the first metabolite is usually straightforward, -Step 1 (SFO for parent only) is skipped here. We start with the model 2a, -with formation and decline of metabolite Z1 and the pathway from parent -directly to sink included (default in mkin). - -<<FOCUS_2006_Z_fits_1, echo=TRUE, fig.height=4>>= -Z.2a <- mkinmod(Z0 = list(type = "SFO", to = "Z1"), - Z1 = list(type = "SFO")) -m.Z.2a <- mkinfit(Z.2a, FOCUS_2006_Z_mkin, quiet = TRUE) -plot(m.Z.2a) -summary(m.Z.2a, data = FALSE) -@ - -As obvious from the summary, the kinetic rate constant from parent compound Z to sink -is negligible. Accordingly, the exact magnitude of the fitted parameter -\texttt{log k\_Z\_sink} is ill-defined and the covariance matrix is not returned. -This suggests, in agreement with the analysis in the FOCUS kinetics report, to simplify -the model by removing the pathway to sink. - -A similar result can be obtained when formation fractions are used in the model -formulation: - -<<FOCUS_2006_Z_fits_2, echo=TRUE, fig.height=4>>= -Z.2a.ff <- mkinmod(Z0 = list(type = "SFO", to = "Z1"), - Z1 = list(type = "SFO"), - use_of_ff = "max") - -m.Z.2a.ff <- mkinfit(Z.2a.ff, FOCUS_2006_Z_mkin, quiet = TRUE) -plot(m.Z.2a.ff) -summary(m.Z.2a.ff, data = FALSE) -@ - -Here, the ilr transformed formation fraction fitted in the model takes a very -large value, and the backtransformed formation fraction from parent Z to Z1 is -practically unity. Again, the covariance matrix is not returned as the model is -overparameterised. - -The simplified model is obtained by setting the list component \texttt{sink} to -\texttt{FALSE}. - -<<FOCUS_2006_Z_fits_3, echo=TRUE, fig.height=4>>= -Z.3 <- mkinmod(Z0 = list(type = "SFO", to = "Z1", sink = FALSE), - Z1 = list(type = "SFO")) -m.Z.3 <- mkinfit(Z.3, FOCUS_2006_Z_mkin, quiet = TRUE) -plot(m.Z.3) -summary(m.Z.3, data = FALSE) -@ - -This model definition is not supported when formation fractions -are used, but the formation fraction can be fixed to unity. - -<<FOCUS_2006_Z_fits_3.ff, echo=TRUE, fig.height=4>>= -Z.3.ff <- mkinmod(Z0 = list(type = "SFO", to = "Z1"), - Z1 = list(type = "SFO"), use_of_ff = "max") -m.Z.3.ff <- mkinfit(Z.3.ff, FOCUS_2006_Z_mkin, - parms.ini = c(f_Z0_to_Z1 = 1), - fixed_parms = "f_Z0_to_Z1", - quiet = TRUE) -summary(m.Z.3.ff, data = FALSE) -@ - -\section{Including metabolites Z2 and Z3} - -As suggested in the FOCUS report, the pathway to sink was removed for metabolite Z1 as -well in the next step. While this step appears questionable on the basis of the above results, it -is followed here for the purpose of comparison. Also, in the FOCUS report, it is -assumed that there is additional empirical evidence that Z1 quickly and exclusively -hydrolyses to Z2. - -<<FOCUS_2006_Z_fits_5, echo=TRUE, fig.height=4>>= -Z.5 <- mkinmod(Z0 = list(type = "SFO", to = "Z1", sink = FALSE), - Z1 = list(type = "SFO", to = "Z2", sink = FALSE), - Z2 = list(type = "SFO")) -m.Z.5 <- mkinfit(Z.5, FOCUS_2006_Z_mkin, quiet = TRUE) -plot(m.Z.5) -summary(m.Z.5, data = FALSE) -@ - -Finally, metabolite Z3 is added to the model. The fit is accellerated -by using the starting parameters from the previous fit. - -<<FOCUS_2006_Z_fits_6, echo=TRUE, fig.height=4>>= -Z.FOCUS <- mkinmod(Z0 = list(type = "SFO", to = "Z1", sink = FALSE), - Z1 = list(type = "SFO", to = "Z2", sink = FALSE), - Z2 = list(type = "SFO", to = "Z3"), - Z3 = list(type = "SFO")) -m.Z.FOCUS <- mkinfit(Z.FOCUS, FOCUS_2006_Z_mkin, - parms.ini = m.Z.5$bparms.ode, - quiet = TRUE) -plot(m.Z.FOCUS) -summary(m.Z.FOCUS, data = FALSE) -@ - -This is the fit corresponding to the final result chosen in Appendix 7 of the -FOCUS report. The residual plots can be obtained by - -<<FOCUS_2006_Z_residuals_6, echo=TRUE>>= -par(mfrow = c(2, 2)) -mkinresplot(m.Z.FOCUS, "Z0", lpos = "bottomright") -mkinresplot(m.Z.FOCUS, "Z1", lpos = "bottomright") -mkinresplot(m.Z.FOCUS, "Z2", lpos = "bottomright") -mkinresplot(m.Z.FOCUS, "Z3", lpos = "bottomright") -@ - -We can also investigate the confidence interval for the formation -fraction from Z1 to Z2 by specifying the model using formation -fractions, and fixing only the formation fraction from Z0 to Z1 -to unity. - -<<FOCUS_2006_Z_fits_6_ff, echo=TRUE, fig.height=4>>= -Z.FOCUS.ff <- mkinmod(Z0 = list(type = "SFO", to = "Z1"), - Z1 = list(type = "SFO", to = "Z2"), - Z2 = list(type = "SFO", to = "Z3"), - Z3 = list(type = "SFO"), use_of_ff = "max") -m.Z.FOCUS.ff <- mkinfit(Z.FOCUS.ff, FOCUS_2006_Z_mkin, - parms.ini = c(f_Z0_to_Z1 = 1), - fixed_parms = c("f_Z0_to_Z1"), quiet = TRUE) -plot(m.Z.FOCUS.ff) -summary(m.Z.FOCUS.ff, data = FALSE) -@ - -\section{Using the SFORB model for parent and metabolites} - -As the FOCUS report states, there is a certain tailing of the time course of metabolite -Z3. Also, the time course of the parent compound is not fitted very well using the -SFO model, as residues at a certain low level remain. - -Therefore, an additional model is offered here, using the single first-order -reversible binding (SFORB) model for metabolite Z3. As expected, the $\chi^2$ -error level is lower for metabolite Z3 using this model and the graphical -fit for Z3 is improved. However, the covariance matrix is not returned. - -<<FOCUS_2006_Z_fits_7, echo=TRUE, fig.height=4>>= -Z.mkin.1 <- mkinmod(Z0 = list(type = "SFO", to = "Z1", sink = FALSE), - Z1 = list(type = "SFO", to = "Z2", sink = FALSE), - Z2 = list(type = "SFO", to = "Z3"), - Z3 = list(type = "SFORB")) -m.Z.mkin.1 <- mkinfit(Z.mkin.1, FOCUS_2006_Z_mkin, - parms.ini = c(k_Z0_Z1 = 0.5, k_Z1_Z2 = 0.3), - quiet = TRUE) -plot(m.Z.mkin.1) -summary(m.Z.mkin.1, data = FALSE) -@ - -Therefore, a further stepwise model building is performed starting from the -stage of parent and one metabolite, starting from the assumption that the model -fit for the parent compound can be improved by using the SFORB model. - -<<FOCUS_2006_Z_fits_8, echo=TRUE, fig.height=4>>= -Z.mkin.2 <- mkinmod(Z0 = list(type = "SFORB", to = "Z1", sink = FALSE), - Z1 = list(type = "SFO")) -m.Z.mkin.2 <- mkinfit(Z.mkin.2, FOCUS_2006_Z_mkin, quiet = TRUE) -plot(m.Z.mkin.2) -summary(m.Z.mkin.2, data = FALSE) -@ - -When metabolite Z2 is added, the additional sink for Z1 is turned off again, -for the same reasons as in the original analysis. - -<<FOCUS_2006_Z_fits_9, echo=TRUE, fig.height=4>>= -Z.mkin.3 <- mkinmod(Z0 = list(type = "SFORB", to = "Z1", sink = FALSE), - Z1 = list(type = "SFO", to = "Z2", sink = FALSE), - Z2 = list(type = "SFO")) -m.Z.mkin.3 <- mkinfit(Z.mkin.3, FOCUS_2006_Z_mkin, quiet = TRUE) -plot(m.Z.mkin.3) -summary(m.Z.mkin.3, data = FALSE) -@ - -This results in a much better representation of the behaviour of the parent -compound Z0. - -Finally, Z3 is added as well. These models appear overparameterised (no -covariance matrix returned) if the sink for Z1 is left in the models. - -<<FOCUS_2006_Z_fits_10, echo=TRUE, fig.height=4>>= -Z.mkin.4 <- mkinmod(Z0 = list(type = "SFORB", to = "Z1", sink = FALSE), - Z1 = list(type = "SFO", to = "Z2", sink = FALSE), - Z2 = list(type = "SFO", to = "Z3"), - Z3 = list(type = "SFO")) -m.Z.mkin.4 <- mkinfit(Z.mkin.4, FOCUS_2006_Z_mkin, - parms.ini = c(k_Z1_Z2 = 0.05), - quiet = TRUE) -plot(m.Z.mkin.4) -summary(m.Z.mkin.4, data = FALSE) -@ - -The error level of the fit, but especially of metabolite Z3, can be improved if -the SFORB model is chosen for this metabolite, as this model is capable of -representing the tailing of the metabolite decline phase. - -<<FOCUS_2006_Z_fits_11, echo=TRUE, fig.height=4>>= -Z.mkin.5 <- mkinmod(Z0 = list(type = "SFORB", to = "Z1", sink = FALSE), - Z1 = list(type = "SFO", to = "Z2", sink = FALSE), - Z2 = list(type = "SFO", to = "Z3"), - Z3 = list(type = "SFORB")) -m.Z.mkin.5 <- mkinfit(Z.mkin.5, FOCUS_2006_Z_mkin, - parms.ini = m.Z.mkin.4$bparms.ode[1:5], - quiet = TRUE) -plot(m.Z.mkin.5) -summary(m.Z.mkin.5, data = FALSE)$bpar -@ - -The summary view of the backtransformed parameters shows that we get no -confidence intervals due to overparameterisation. As the optimized -\texttt{k\_Z3\_bound\_free} is excessively small, it is reasonable to fix it to -zero. - -<<FOCUS_2006_Z_fits_11a, echo=TRUE>>= -m.Z.mkin.5a <- mkinfit(Z.mkin.5, FOCUS_2006_Z_mkin, - parms.ini = c(m.Z.mkin.4$bparms.ode[1:5], - k_Z3_bound_free = 0), - fixed_parms = "k_Z3_bound_free", - quiet = TRUE) -summary(m.Z.mkin.5a, data = FALSE)$bpar -@ - -A graphical representation of the confidence intervals can finally be obtained. -<<FOCUS_2006_Z_fits_11b, echo=TRUE>>= -mkinparplot(m.Z.mkin.5a) -@ - -It is clear that nothing can be said about the degradation rate of Z3 towards -the end of the experiment. However, this appears to be a feature of the data. - -<<FOCUS_2006_Z_residuals_11>>= -par(mfrow = c(2, 2)) -mkinresplot(m.Z.mkin.5, "Z0", lpos = "bottomright") -mkinresplot(m.Z.mkin.5, "Z1", lpos = "bottomright") -mkinresplot(m.Z.mkin.5, "Z2", lpos = "bottomright") -mkinresplot(m.Z.mkin.5, "Z3", lpos = "bottomright") -@ - -As expected, the residual plots are much more random than in the case of the -all SFO model for which they were shown above. In conclusion, the model -\texttt{Z.mkin.5} is proposed as the best-fit model for the dataset from -Appendix 7 of the FOCUS report. - -\bibliographystyle{plainnat} -\bibliography{references} - -\end{document} -% vim: set foldmethod=syntax: diff --git a/vignettes/FOCUS_Z.pdf b/vignettes/FOCUS_Z.pdf Binary files differdeleted file mode 100644 index e6138d0..0000000 --- a/vignettes/FOCUS_Z.pdf +++ /dev/null diff --git a/vignettes/figure/FOCUS_2006_Z_fits_1.pdf b/vignettes/figure/FOCUS_2006_Z_fits_1.pdf Binary files differdeleted file mode 100644 index 7533073..0000000 --- a/vignettes/figure/FOCUS_2006_Z_fits_1.pdf +++ /dev/null diff --git a/vignettes/figure/FOCUS_2006_Z_fits_10.pdf b/vignettes/figure/FOCUS_2006_Z_fits_10.pdf Binary files differdeleted file mode 100644 index 5ebd381..0000000 --- a/vignettes/figure/FOCUS_2006_Z_fits_10.pdf +++ /dev/null diff --git a/vignettes/figure/FOCUS_2006_Z_fits_11.pdf b/vignettes/figure/FOCUS_2006_Z_fits_11.pdf Binary files differdeleted file mode 100644 index 1ccff78..0000000 --- a/vignettes/figure/FOCUS_2006_Z_fits_11.pdf +++ /dev/null diff --git a/vignettes/figure/FOCUS_2006_Z_fits_11b.pdf b/vignettes/figure/FOCUS_2006_Z_fits_11b.pdf Binary files differdeleted file mode 100644 index 4224b93..0000000 --- a/vignettes/figure/FOCUS_2006_Z_fits_11b.pdf +++ /dev/null diff --git a/vignettes/figure/FOCUS_2006_Z_fits_2.pdf b/vignettes/figure/FOCUS_2006_Z_fits_2.pdf Binary files differdeleted file mode 100644 index fef7b27..0000000 --- a/vignettes/figure/FOCUS_2006_Z_fits_2.pdf +++ /dev/null diff --git a/vignettes/figure/FOCUS_2006_Z_fits_3.pdf b/vignettes/figure/FOCUS_2006_Z_fits_3.pdf Binary files differdeleted file mode 100644 index 4896136..0000000 --- a/vignettes/figure/FOCUS_2006_Z_fits_3.pdf +++ /dev/null diff --git a/vignettes/figure/FOCUS_2006_Z_fits_5.pdf b/vignettes/figure/FOCUS_2006_Z_fits_5.pdf Binary files differdeleted file mode 100644 index 91a06ba..0000000 --- a/vignettes/figure/FOCUS_2006_Z_fits_5.pdf +++ /dev/null diff --git a/vignettes/figure/FOCUS_2006_Z_fits_6.pdf b/vignettes/figure/FOCUS_2006_Z_fits_6.pdf Binary files differdeleted file mode 100644 index 967da63..0000000 --- a/vignettes/figure/FOCUS_2006_Z_fits_6.pdf +++ /dev/null diff --git a/vignettes/figure/FOCUS_2006_Z_fits_6_ff.pdf b/vignettes/figure/FOCUS_2006_Z_fits_6_ff.pdf Binary files differdeleted file mode 100644 index c4fe094..0000000 --- a/vignettes/figure/FOCUS_2006_Z_fits_6_ff.pdf +++ /dev/null diff --git a/vignettes/figure/FOCUS_2006_Z_fits_7.pdf b/vignettes/figure/FOCUS_2006_Z_fits_7.pdf Binary files differdeleted file mode 100644 index f40ea63..0000000 --- a/vignettes/figure/FOCUS_2006_Z_fits_7.pdf +++ /dev/null diff --git a/vignettes/figure/FOCUS_2006_Z_fits_8.pdf b/vignettes/figure/FOCUS_2006_Z_fits_8.pdf Binary files differdeleted file mode 100644 index e0f9b94..0000000 --- a/vignettes/figure/FOCUS_2006_Z_fits_8.pdf +++ /dev/null diff --git a/vignettes/figure/FOCUS_2006_Z_fits_9.pdf b/vignettes/figure/FOCUS_2006_Z_fits_9.pdf Binary files differdeleted file mode 100644 index 9b6f570..0000000 --- a/vignettes/figure/FOCUS_2006_Z_fits_9.pdf +++ /dev/null diff --git a/vignettes/figure/FOCUS_2006_Z_residuals_11.pdf b/vignettes/figure/FOCUS_2006_Z_residuals_11.pdf Binary files differdeleted file mode 100644 index 5a41904..0000000 --- a/vignettes/figure/FOCUS_2006_Z_residuals_11.pdf +++ /dev/null diff --git a/vignettes/figure/FOCUS_2006_Z_residuals_6.pdf b/vignettes/figure/FOCUS_2006_Z_residuals_6.pdf Binary files differdeleted file mode 100644 index 3e6649d..0000000 --- a/vignettes/figure/FOCUS_2006_Z_residuals_6.pdf +++ /dev/null diff --git a/vignettes/figure/unnamed-chunk-10.png b/vignettes/figure/unnamed-chunk-10.png Binary files differdeleted file mode 100644 index cd3a700..0000000 --- a/vignettes/figure/unnamed-chunk-10.png +++ /dev/null diff --git a/vignettes/figure/unnamed-chunk-11.png b/vignettes/figure/unnamed-chunk-11.png Binary files differdeleted file mode 100644 index ca488e6..0000000 --- a/vignettes/figure/unnamed-chunk-11.png +++ /dev/null diff --git a/vignettes/figure/unnamed-chunk-12.png b/vignettes/figure/unnamed-chunk-12.png Binary files differdeleted file mode 100644 index 3a64413..0000000 --- a/vignettes/figure/unnamed-chunk-12.png +++ /dev/null diff --git a/vignettes/figure/unnamed-chunk-14.png b/vignettes/figure/unnamed-chunk-14.png Binary files differdeleted file mode 100644 index 46d9c50..0000000 --- a/vignettes/figure/unnamed-chunk-14.png +++ /dev/null diff --git a/vignettes/figure/unnamed-chunk-15.png b/vignettes/figure/unnamed-chunk-15.png Binary files differdeleted file mode 100644 index 0eddbc6..0000000 --- a/vignettes/figure/unnamed-chunk-15.png +++ /dev/null diff --git a/vignettes/figure/unnamed-chunk-16.png b/vignettes/figure/unnamed-chunk-16.png Binary files differdeleted file mode 100644 index 4d1b738..0000000 --- a/vignettes/figure/unnamed-chunk-16.png +++ /dev/null diff --git a/vignettes/figure/unnamed-chunk-18.png b/vignettes/figure/unnamed-chunk-18.png Binary files differdeleted file mode 100644 index b109a11..0000000 --- a/vignettes/figure/unnamed-chunk-18.png +++ /dev/null diff --git a/vignettes/figure/unnamed-chunk-19.png b/vignettes/figure/unnamed-chunk-19.png Binary files differdeleted file mode 100644 index af84c2a..0000000 --- a/vignettes/figure/unnamed-chunk-19.png +++ /dev/null diff --git a/vignettes/figure/unnamed-chunk-4.png b/vignettes/figure/unnamed-chunk-4.png Binary files differdeleted file mode 100644 index 04187f8..0000000 --- a/vignettes/figure/unnamed-chunk-4.png +++ /dev/null diff --git a/vignettes/figure/unnamed-chunk-5.png b/vignettes/figure/unnamed-chunk-5.png Binary files differdeleted file mode 100644 index f40ba5c..0000000 --- a/vignettes/figure/unnamed-chunk-5.png +++ /dev/null diff --git a/vignettes/figure/unnamed-chunk-9.png b/vignettes/figure/unnamed-chunk-9.png Binary files differdeleted file mode 100644 index 76fd0c3..0000000 --- a/vignettes/figure/unnamed-chunk-9.png +++ /dev/null diff --git a/vignettes/header.tex b/vignettes/header.tex deleted file mode 100644 index 476415e..0000000 --- a/vignettes/header.tex +++ /dev/null @@ -1,23 +0,0 @@ -\usepackage{booktabs}
-\usepackage{amsfonts}
-\usepackage{latexsym}
-\usepackage{amsmath}
-\usepackage{amssymb}
-\usepackage{graphicx}
-\usepackage{parskip}
-\usepackage[round]{natbib}
-\usepackage{amstext}
-\usepackage{hyperref}
-\usepackage[utf8]{inputenc}
-
-\newcommand{\Rpackage}[1]{{\normalfont\fontseries{b}\selectfont #1}}
-\newcommand{\Robject}[1]{\texttt{#1}}
-\newcommand{\Rclass}[1]{\textit{#1}}
-\newcommand{\Rcmd}[1]{\texttt{#1}}
-
-\newcommand{\RR}{\textsf{R}}
-
-\RequirePackage[T1]{fontenc}
-\RequirePackage{graphicx,ae,fancyvrb}
-\IfFileExists{upquote.sty}{\RequirePackage{upquote}}{}
-\usepackage{relsize}
diff --git a/vignettes/mkin.Rnw b/vignettes/mkin.Rnw deleted file mode 100644 index 0ac114e..0000000 --- a/vignettes/mkin.Rnw +++ /dev/null @@ -1,74 +0,0 @@ -%\VignetteIndexEntry{Routines for fitting kinetic models with one or more state variables to chemical degradation data} -%\VignetteEngine{knitr::knitr} -\documentclass[12pt,a4paper]{article} -\usepackage{a4wide} -\input{header} -\hypersetup{ - pdftitle = {mkin - Routines for fitting kinetic models with one or more state variables to chemical degradation data}, - pdfsubject = {Manuscript}, - pdfauthor = {Johannes Ranke}, - colorlinks = {true}, - linkcolor = {blue}, - citecolor = {blue}, - urlcolor = {red}, - linktocpage = {true}, -} - -\begin{document} - -<<include=FALSE>>= -require(knitr) -opts_chunk$set(engine='R', tidy=FALSE) -@ - -\title{mkin -\\ -Routines for fitting kinetic models with one or more state variables to -chemical degradation data} -\author{\textbf{Johannes Ranke} \\[0.5cm] -%EndAName -Wissenschaftlicher Berater\\ -Kronacher Str. 8, 79639 Grenzach-Wyhlen, Germany\\[0.5cm] -and\\[0.5cm] -University of Bremen\\ -} -\maketitle - -\begin{abstract} -In the regulatory evaluation of chemical substances like plant protection -products (pesticides), biocides and other chemicals, degradation data play an -important role. For the evaluation of pesticide degradation experiments, -detailed guidance has been developed, based on nonlinear optimisation. -The \RR{} add-on package \Rpackage{mkin} implements fitting some of the models -recommended in this guidance from within R and calculates some statistical -measures for data series within one or more compartments, for parent and -metabolites. -\end{abstract} - - -\thispagestyle{empty} \setcounter{page}{0} - -\clearpage - -\tableofcontents - -\textbf{Key words}: Kinetics, FOCUS, nonlinear optimisation - -\section{Introduction} -\label{intro} - -Many approaches are possible regarding the evaluation of chemical degradation -data. The \Rpackage{kinfit} package \citep{pkg:kinfit} in \RR{} -\citep{rcore2013} implements the approach recommended in the kinetics report -provided by the FOrum for Co-ordination of pesticide fate models and their -USe \citep{FOCUS2006, FOCUSkinetics2011} for simple data series for one parent -compound in one compartment. - -The \Rpackage{mkin} package \citep{pkg:mkin} extends this approach to data series -with metabolites and more than one compartment and includes the possibility -for back reactions. - -\bibliographystyle{plainnat} -\bibliography{references} - -\end{document} -% vim: set foldmethod=syntax: diff --git a/vignettes/mkin.pdf b/vignettes/mkin.pdf Binary files differdeleted file mode 100644 index 42a44d6..0000000 --- a/vignettes/mkin.pdf +++ /dev/null diff --git a/vignettes/references.bib b/vignettes/references.bib deleted file mode 100644 index e069daf..0000000 --- a/vignettes/references.bib +++ /dev/null @@ -1,79 +0,0 @@ -% This file was originally created with JabRef 2.7b. -% Encoding: ISO8859_1 - -@BOOK{bates88, - title = {Nonlinear regression and its applications}, - publisher = {Wiley-Interscience}, - year = {1988}, - author = {D. Bates and D. Watts} -} - -@MANUAL{FOCUSkinetics2011, - title = {Generic guidance for estimating persistence and degradation kinetics - from environmental fate studies on pesticides in EU registration}, - author = {{FOCUS Work Group on Degradation Kinetics}}, - edition = {1.0}, - month = {November}, - year = {2011}, - file = {FOCUS kinetics 2011 Generic guidance:/home/ranke/dok/orgs/focus/FOCUSkineticsvc_1_0_Nov23.pdf:PDF}, - url = {http://focus.jrc.ec.europa.eu/dk} -} - -@MANUAL{FOCUS2006, - title = {Guidance Document on Estimating Persistence and Degradation Kinetics - from Environmental Fate Studies on Pesticides in EU Registration. - Report of the FOCUS Work Group on Degradation Kinetics}, - author = {{FOCUS Work Group on Degradation Kinetics}}, - year = {2006}, - note = {EC Document Reference Sanco/10058/2005 version 2.0}, - url = {http://focus.jrc.ec.europa.eu/dk} -} - -@MANUAL{rcore2013, - title = {\textsf{R}: A Language and Environment for Statistical Computing}, - author = {{R Development Core Team}}, - organization = {R Foundation for Statistical Computing}, - address = {Vienna, Austria}, - year = {2013}, - note = {{ISBN} 3-900051-07-0}, - url = {http://www.R-project.org} -} - -@MANUAL{pkg:kinfit, - title = {kinfit: {R}outines for fitting simple kinetic models to chemical - degradation data}, - author = {Johannes Ranke}, - year = {2013}, - url = {http://CRAN.R-project.org/package=kinfit} -} - -@MANUAL{pkg:mkin, - title = {mkin: {R}outines for fitting kinetic models with one or more state - variables to chemical degradation data}, - author = {Johannes Ranke}, - year = {2013}, - url = {http://CRAN.R-project.org/package/kinfit} -} - -@INPROCEEDINGS{schaefer2007, - author = {D. Sch\"{a}fer and M. Mikolasch and P. Rainbird and B. Harvey}, - title = {{KinGUI}: a new kinetic software tool for evaluations according to - FOCUS degradation kinetics}, - booktitle = {Proceedings of the XIII Symposium Pesticide Chemistry}, - year = {2007}, - editor = {Del Re A. A. M. and Capri E. and Fragoulis G. and Trevisan M.}, - pages = {916--923}, - address = {Piacenza} -} - -@ARTICLE{soetaert10, - author = {Karline Soetaert and Thomas Petzoldt}, - title = {Inverse Modelling, Sensitivity and Monte Carlo Analysis in {R} Using - Package {FME}}, - journal = {Journal of Statistical Software}, - year = {2010}, - volume = {33}, - pages = {1--28}, - number = {3}, - url = {http://www.jstatsoft.org/v33/i03/} -} |
