From e959fde98f95f3595e01490b67892678bbcd1b27 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 7 May 2014 14:47:28 +0200 Subject: Fork the gmkin GUI from mkin. See ChangeLog for details --- ChangeLog | 8 + DESCRIPTION | 31 +- GNUmakefile | 52 +- NAMESPACE | 7 +- R/DFOP.solution.R | 5 - R/FOMC.solution.R | 4 - R/HS.solution.R | 6 - R/SFO.solution.R | 4 - R/SFORB.solution.R | 9 - R/endpoints.R | 108 --- R/gmkin.R | 10 +- R/ilr.R | 51 -- R/mkin_long_to_wide.R | 30 - R/mkin_wide_to_long.R | 34 - R/mkinerrmin.R | 85 --- R/mkinfit.R | 482 ------------- R/mkinmod.R | 242 ------- R/mkinparplot.R | 62 -- R/mkinplot.R | 4 - R/mkinpredict.R | 121 ---- R/mkinresplot.R | 53 -- R/plot.mkinfit.R | 98 --- R/transform_odeparms.R | 101 --- README.md | 182 +---- TODO | 14 +- data/FOCUS_2006_DFOP_ref_A_to_B.rda | Bin 723 -> 0 bytes data/FOCUS_2006_FOMC_ref_A_to_F.rda | Bin 1227 -> 0 bytes data/FOCUS_2006_HS_ref_A_to_F.rda | Bin 1360 -> 0 bytes data/FOCUS_2006_SFO_ref_A_to_F.rda | Bin 1097 -> 0 bytes data/FOCUS_2006_datasets.RData | Bin 905 -> 0 bytes data/mccall81_245T.RData | Bin 759 -> 0 bytes data/schaefer07_complex_case.RData | Bin 828 -> 0 bytes inst/GUI/README | 1 - inst/GUI/TODO | 1 - inst/GUI/gmkin.R | 1 - inst/unitTests/Makefile | 15 - inst/unitTests/runit.mkinfit.R | 294 -------- inst/unitTests/runit.mkinpredict.R | 86 --- man/DFOP.solution.Rd | 36 - man/FOCUS_2006_DFOP_ref_A_to_B.Rd | 39 -- man/FOCUS_2006_FOMC_ref_A_to_F.Rd | 38 - man/FOCUS_2006_HS_ref_A_to_F.Rd | 39 -- man/FOCUS_2006_SFO_ref_A_to_F.Rd | 37 - man/FOCUS_2006_datasets.Rd | 35 - man/FOMC.solution.Rd | 49 -- man/HS.solution.Rd | 34 - man/SFO.solution.Rd | 29 - man/SFORB.solution.Rd | 36 - man/endpoints.Rd | 33 - man/ilr.Rd | 56 -- man/mccall81_245T.Rd | 42 -- man/mkin_long_to_wide.Rd | 36 - man/mkin_wide_to_long.Rd | 32 - man/mkinerrmin.Rd | 44 -- man/mkinfit.Rd | 254 ------- man/mkinmod.Rd | 63 -- man/mkinparplot.Rd | 35 - man/mkinplot.Rd | 26 - man/mkinpredict.Rd | 91 --- man/mkinresplot.Rd | 67 -- man/plot.mkinfit.Rd | 94 --- man/schaefer07_complex_case.Rd | 42 -- man/summary.mkinfit.Rd | 75 -- man/transform_odeparms.Rd | 69 -- tests/doRUnit.R | 43 -- vignettes/FOCUS_L.Rmd | 243 ------- vignettes/FOCUS_L.md | 931 ------------------------- vignettes/FOCUS_Z.Rnw | 310 -------- vignettes/FOCUS_Z.pdf | Bin 218511 -> 0 bytes vignettes/figure/FOCUS_2006_Z_fits_1.pdf | Bin 5670 -> 0 bytes vignettes/figure/FOCUS_2006_Z_fits_10.pdf | Bin 7587 -> 0 bytes vignettes/figure/FOCUS_2006_Z_fits_11.pdf | Bin 7598 -> 0 bytes vignettes/figure/FOCUS_2006_Z_fits_11b.pdf | Bin 5021 -> 0 bytes vignettes/figure/FOCUS_2006_Z_fits_2.pdf | Bin 5670 -> 0 bytes vignettes/figure/FOCUS_2006_Z_fits_3.pdf | Bin 5670 -> 0 bytes vignettes/figure/FOCUS_2006_Z_fits_5.pdf | Bin 6155 -> 0 bytes vignettes/figure/FOCUS_2006_Z_fits_6.pdf | Bin 6966 -> 0 bytes vignettes/figure/FOCUS_2006_Z_fits_6_ff.pdf | Bin 6937 -> 0 bytes vignettes/figure/FOCUS_2006_Z_fits_7.pdf | Bin 6967 -> 0 bytes vignettes/figure/FOCUS_2006_Z_fits_8.pdf | Bin 6158 -> 0 bytes vignettes/figure/FOCUS_2006_Z_fits_9.pdf | Bin 6785 -> 0 bytes vignettes/figure/FOCUS_2006_Z_residuals_11.pdf | Bin 5940 -> 0 bytes vignettes/figure/FOCUS_2006_Z_residuals_6.pdf | Bin 5959 -> 0 bytes vignettes/figure/unnamed-chunk-10.png | Bin 8140 -> 0 bytes vignettes/figure/unnamed-chunk-11.png | Bin 5089 -> 0 bytes vignettes/figure/unnamed-chunk-12.png | Bin 5518 -> 0 bytes vignettes/figure/unnamed-chunk-14.png | Bin 5690 -> 0 bytes vignettes/figure/unnamed-chunk-15.png | Bin 5742 -> 0 bytes vignettes/figure/unnamed-chunk-16.png | Bin 5691 -> 0 bytes vignettes/figure/unnamed-chunk-18.png | Bin 5373 -> 0 bytes vignettes/figure/unnamed-chunk-19.png | Bin 5333 -> 0 bytes vignettes/figure/unnamed-chunk-4.png | Bin 5950 -> 0 bytes vignettes/figure/unnamed-chunk-5.png | Bin 4665 -> 0 bytes vignettes/figure/unnamed-chunk-9.png | Bin 7851 -> 0 bytes vignettes/header.tex | 23 - vignettes/mkin.Rnw | 74 -- vignettes/mkin.pdf | Bin 124442 -> 0 bytes vignettes/references.bib | 79 --- 98 files changed, 43 insertions(+), 5292 deletions(-) delete mode 100644 R/DFOP.solution.R delete mode 100644 R/FOMC.solution.R delete mode 100644 R/HS.solution.R delete mode 100644 R/SFO.solution.R delete mode 100644 R/SFORB.solution.R delete mode 100644 R/endpoints.R delete mode 100644 R/ilr.R delete mode 100644 R/mkin_long_to_wide.R delete mode 100644 R/mkin_wide_to_long.R delete mode 100644 R/mkinerrmin.R delete mode 100644 R/mkinfit.R delete mode 100644 R/mkinmod.R delete mode 100644 R/mkinparplot.R delete mode 100644 R/mkinplot.R delete mode 100644 R/mkinpredict.R delete mode 100644 R/mkinresplot.R delete mode 100644 R/plot.mkinfit.R delete mode 100644 R/transform_odeparms.R delete mode 100644 data/FOCUS_2006_DFOP_ref_A_to_B.rda delete mode 100644 data/FOCUS_2006_FOMC_ref_A_to_F.rda delete mode 100644 data/FOCUS_2006_HS_ref_A_to_F.rda delete mode 100644 data/FOCUS_2006_SFO_ref_A_to_F.rda delete mode 100644 data/FOCUS_2006_datasets.RData delete mode 100644 data/mccall81_245T.RData delete mode 100644 data/schaefer07_complex_case.RData delete mode 100644 inst/GUI/README delete mode 100644 inst/GUI/TODO delete mode 100644 inst/unitTests/Makefile delete mode 100644 inst/unitTests/runit.mkinfit.R delete mode 100644 inst/unitTests/runit.mkinpredict.R delete mode 100644 man/DFOP.solution.Rd delete mode 100644 man/FOCUS_2006_DFOP_ref_A_to_B.Rd delete mode 100644 man/FOCUS_2006_FOMC_ref_A_to_F.Rd delete mode 100644 man/FOCUS_2006_HS_ref_A_to_F.Rd delete mode 100644 man/FOCUS_2006_SFO_ref_A_to_F.Rd delete mode 100644 man/FOCUS_2006_datasets.Rd delete mode 100644 man/FOMC.solution.Rd delete mode 100644 man/HS.solution.Rd delete mode 100644 man/SFO.solution.Rd delete mode 100644 man/SFORB.solution.Rd delete mode 100644 man/endpoints.Rd delete mode 100644 man/ilr.Rd delete mode 100644 man/mccall81_245T.Rd delete mode 100644 man/mkin_long_to_wide.Rd delete mode 100644 man/mkin_wide_to_long.Rd delete mode 100644 man/mkinerrmin.Rd delete mode 100644 man/mkinfit.Rd delete mode 100644 man/mkinmod.Rd delete mode 100644 man/mkinparplot.Rd delete mode 100644 man/mkinplot.Rd delete mode 100644 man/mkinpredict.Rd delete mode 100644 man/mkinresplot.Rd delete mode 100644 man/plot.mkinfit.Rd delete mode 100644 man/schaefer07_complex_case.Rd delete mode 100644 man/summary.mkinfit.Rd delete mode 100644 man/transform_odeparms.Rd delete mode 100644 tests/doRUnit.R delete mode 100644 vignettes/FOCUS_L.Rmd delete mode 100644 vignettes/FOCUS_L.md delete mode 100644 vignettes/FOCUS_Z.Rnw delete mode 100644 vignettes/FOCUS_Z.pdf delete mode 100644 vignettes/figure/FOCUS_2006_Z_fits_1.pdf delete mode 100644 vignettes/figure/FOCUS_2006_Z_fits_10.pdf delete mode 100644 vignettes/figure/FOCUS_2006_Z_fits_11.pdf delete mode 100644 vignettes/figure/FOCUS_2006_Z_fits_11b.pdf delete mode 100644 vignettes/figure/FOCUS_2006_Z_fits_2.pdf delete mode 100644 vignettes/figure/FOCUS_2006_Z_fits_3.pdf delete mode 100644 vignettes/figure/FOCUS_2006_Z_fits_5.pdf delete mode 100644 vignettes/figure/FOCUS_2006_Z_fits_6.pdf delete mode 100644 vignettes/figure/FOCUS_2006_Z_fits_6_ff.pdf delete mode 100644 vignettes/figure/FOCUS_2006_Z_fits_7.pdf delete mode 100644 vignettes/figure/FOCUS_2006_Z_fits_8.pdf delete mode 100644 vignettes/figure/FOCUS_2006_Z_fits_9.pdf delete mode 100644 vignettes/figure/FOCUS_2006_Z_residuals_11.pdf delete mode 100644 vignettes/figure/FOCUS_2006_Z_residuals_6.pdf delete mode 100644 vignettes/figure/unnamed-chunk-10.png delete mode 100644 vignettes/figure/unnamed-chunk-11.png delete mode 100644 vignettes/figure/unnamed-chunk-12.png delete mode 100644 vignettes/figure/unnamed-chunk-14.png delete mode 100644 vignettes/figure/unnamed-chunk-15.png delete mode 100644 vignettes/figure/unnamed-chunk-16.png delete mode 100644 vignettes/figure/unnamed-chunk-18.png delete mode 100644 vignettes/figure/unnamed-chunk-19.png delete mode 100644 vignettes/figure/unnamed-chunk-4.png delete mode 100644 vignettes/figure/unnamed-chunk-5.png delete mode 100644 vignettes/figure/unnamed-chunk-9.png delete mode 100644 vignettes/header.tex delete mode 100644 vignettes/mkin.Rnw delete mode 100644 vignettes/mkin.pdf delete mode 100644 vignettes/references.bib diff --git a/ChangeLog b/ChangeLog index 2bbc857..70ba083 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,11 @@ +2014-05-07 Johannes Ranke 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 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 diff --git a/NAMESPACE b/NAMESPACE index dfa9e9b..4ddb2b5 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -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) -} diff --git a/R/gmkin.R b/R/gmkin.R index 912d31c..31e8cf4 100644 --- a/R/gmkin.R +++ b/R/gmkin.R @@ -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 - -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 - -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 -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 -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 -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 }}} - -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 -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 - -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 -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 -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 - -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) -} diff --git a/README.md b/README.md index 98fc276..9d1f793 100644 --- a/README.md +++ b/README.md @@ -1,121 +1,23 @@ -# mkin +# gmkin -The R package **mkin** provides calculation routines for the analysis of -chemical degradation data, including multicompartment kinetics 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") ``` ![gmkin screenshot](gmkin_screenshot.png) - -## 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. diff --git a/TODO b/TODO index 944ad57..8bd0c1d 100644 --- a/TODO +++ b/TODO @@ -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 deleted file mode 100644 index 8b9dc9f..0000000 Binary files a/data/FOCUS_2006_DFOP_ref_A_to_B.rda and /dev/null differ diff --git a/data/FOCUS_2006_FOMC_ref_A_to_F.rda b/data/FOCUS_2006_FOMC_ref_A_to_F.rda deleted file mode 100644 index 246a53f..0000000 Binary files a/data/FOCUS_2006_FOMC_ref_A_to_F.rda and /dev/null differ diff --git a/data/FOCUS_2006_HS_ref_A_to_F.rda b/data/FOCUS_2006_HS_ref_A_to_F.rda deleted file mode 100644 index 1e3b386..0000000 Binary files a/data/FOCUS_2006_HS_ref_A_to_F.rda and /dev/null differ diff --git a/data/FOCUS_2006_SFO_ref_A_to_F.rda b/data/FOCUS_2006_SFO_ref_A_to_F.rda deleted file mode 100644 index 411c536..0000000 Binary files a/data/FOCUS_2006_SFO_ref_A_to_F.rda and /dev/null differ diff --git a/data/FOCUS_2006_datasets.RData b/data/FOCUS_2006_datasets.RData deleted file mode 100644 index eb9e526..0000000 Binary files a/data/FOCUS_2006_datasets.RData and /dev/null differ diff --git a/data/mccall81_245T.RData b/data/mccall81_245T.RData deleted file mode 100644 index 530d83b..0000000 Binary files a/data/mccall81_245T.RData and /dev/null differ diff --git a/data/schaefer07_complex_case.RData b/data/schaefer07_complex_case.RData deleted file mode 100644 index d58aedd..0000000 Binary files a/data/schaefer07_complex_case.RData and /dev/null differ 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 -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 - -# 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 - -# 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 -} -\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 -} -\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 @@ - - -# 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 @@ - - -# 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) -``` - -![plot of chunk unnamed-chunk-4](figure/unnamed-chunk-4.png) - -The residual plot can be easily obtained by - - -```r -mkinresplot(m.L1.SFO, ylab = "Observed", xlab = "Time") -``` - -![plot of chunk unnamed-chunk-5](figure/unnamed-chunk-5.png) - - -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) -``` - -![plot of chunk unnamed-chunk-9](figure/unnamed-chunk-9.png) - - -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) -``` - -![plot of chunk unnamed-chunk-10](figure/unnamed-chunk-10.png) - -```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) -``` - -![plot of chunk unnamed-chunk-11](figure/unnamed-chunk-11.png) - - -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) -``` - -![plot of chunk unnamed-chunk-12](figure/unnamed-chunk-12.png) - -```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) -``` - -![plot of chunk unnamed-chunk-14](figure/unnamed-chunk-14.png) - -```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) -``` - -![plot of chunk unnamed-chunk-15](figure/unnamed-chunk-15.png) - -```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) -``` - -![plot of chunk unnamed-chunk-16](figure/unnamed-chunk-16.png) - -```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) -``` - -![plot of chunk unnamed-chunk-18](figure/unnamed-chunk-18.png) - -```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) -``` - -![plot of chunk unnamed-chunk-19](figure/unnamed-chunk-19.png) - -```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} - -<>= -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. - -<>= -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). - -<>= -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: - -<>= -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}. - -<>= -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. - -<>= -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. - -<>= -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. - -<>= -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 - -<>= -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. - -<>= -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. - -<>= -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. - -<>= -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. - -<>= -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. - -<>= -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. - -<>= -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. - -<>= -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. -<>= -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. - -<>= -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 deleted file mode 100644 index e6138d0..0000000 Binary files a/vignettes/FOCUS_Z.pdf and /dev/null differ diff --git a/vignettes/figure/FOCUS_2006_Z_fits_1.pdf b/vignettes/figure/FOCUS_2006_Z_fits_1.pdf deleted file mode 100644 index 7533073..0000000 Binary files a/vignettes/figure/FOCUS_2006_Z_fits_1.pdf and /dev/null differ diff --git a/vignettes/figure/FOCUS_2006_Z_fits_10.pdf b/vignettes/figure/FOCUS_2006_Z_fits_10.pdf deleted file mode 100644 index 5ebd381..0000000 Binary files a/vignettes/figure/FOCUS_2006_Z_fits_10.pdf and /dev/null differ diff --git a/vignettes/figure/FOCUS_2006_Z_fits_11.pdf b/vignettes/figure/FOCUS_2006_Z_fits_11.pdf deleted file mode 100644 index 1ccff78..0000000 Binary files a/vignettes/figure/FOCUS_2006_Z_fits_11.pdf and /dev/null differ diff --git a/vignettes/figure/FOCUS_2006_Z_fits_11b.pdf b/vignettes/figure/FOCUS_2006_Z_fits_11b.pdf deleted file mode 100644 index 4224b93..0000000 Binary files a/vignettes/figure/FOCUS_2006_Z_fits_11b.pdf and /dev/null differ diff --git a/vignettes/figure/FOCUS_2006_Z_fits_2.pdf b/vignettes/figure/FOCUS_2006_Z_fits_2.pdf deleted file mode 100644 index fef7b27..0000000 Binary files a/vignettes/figure/FOCUS_2006_Z_fits_2.pdf and /dev/null differ diff --git a/vignettes/figure/FOCUS_2006_Z_fits_3.pdf b/vignettes/figure/FOCUS_2006_Z_fits_3.pdf deleted file mode 100644 index 4896136..0000000 Binary files a/vignettes/figure/FOCUS_2006_Z_fits_3.pdf and /dev/null differ diff --git a/vignettes/figure/FOCUS_2006_Z_fits_5.pdf b/vignettes/figure/FOCUS_2006_Z_fits_5.pdf deleted file mode 100644 index 91a06ba..0000000 Binary files a/vignettes/figure/FOCUS_2006_Z_fits_5.pdf and /dev/null differ diff --git a/vignettes/figure/FOCUS_2006_Z_fits_6.pdf b/vignettes/figure/FOCUS_2006_Z_fits_6.pdf deleted file mode 100644 index 967da63..0000000 Binary files a/vignettes/figure/FOCUS_2006_Z_fits_6.pdf and /dev/null differ diff --git a/vignettes/figure/FOCUS_2006_Z_fits_6_ff.pdf b/vignettes/figure/FOCUS_2006_Z_fits_6_ff.pdf deleted file mode 100644 index c4fe094..0000000 Binary files a/vignettes/figure/FOCUS_2006_Z_fits_6_ff.pdf and /dev/null differ diff --git a/vignettes/figure/FOCUS_2006_Z_fits_7.pdf b/vignettes/figure/FOCUS_2006_Z_fits_7.pdf deleted file mode 100644 index f40ea63..0000000 Binary files a/vignettes/figure/FOCUS_2006_Z_fits_7.pdf and /dev/null differ diff --git a/vignettes/figure/FOCUS_2006_Z_fits_8.pdf b/vignettes/figure/FOCUS_2006_Z_fits_8.pdf deleted file mode 100644 index e0f9b94..0000000 Binary files a/vignettes/figure/FOCUS_2006_Z_fits_8.pdf and /dev/null differ diff --git a/vignettes/figure/FOCUS_2006_Z_fits_9.pdf b/vignettes/figure/FOCUS_2006_Z_fits_9.pdf deleted file mode 100644 index 9b6f570..0000000 Binary files a/vignettes/figure/FOCUS_2006_Z_fits_9.pdf and /dev/null differ diff --git a/vignettes/figure/FOCUS_2006_Z_residuals_11.pdf b/vignettes/figure/FOCUS_2006_Z_residuals_11.pdf deleted file mode 100644 index 5a41904..0000000 Binary files a/vignettes/figure/FOCUS_2006_Z_residuals_11.pdf and /dev/null differ diff --git a/vignettes/figure/FOCUS_2006_Z_residuals_6.pdf b/vignettes/figure/FOCUS_2006_Z_residuals_6.pdf deleted file mode 100644 index 3e6649d..0000000 Binary files a/vignettes/figure/FOCUS_2006_Z_residuals_6.pdf and /dev/null differ diff --git a/vignettes/figure/unnamed-chunk-10.png b/vignettes/figure/unnamed-chunk-10.png deleted file mode 100644 index cd3a700..0000000 Binary files a/vignettes/figure/unnamed-chunk-10.png and /dev/null differ diff --git a/vignettes/figure/unnamed-chunk-11.png b/vignettes/figure/unnamed-chunk-11.png deleted file mode 100644 index ca488e6..0000000 Binary files a/vignettes/figure/unnamed-chunk-11.png and /dev/null differ diff --git a/vignettes/figure/unnamed-chunk-12.png b/vignettes/figure/unnamed-chunk-12.png deleted file mode 100644 index 3a64413..0000000 Binary files a/vignettes/figure/unnamed-chunk-12.png and /dev/null differ diff --git a/vignettes/figure/unnamed-chunk-14.png b/vignettes/figure/unnamed-chunk-14.png deleted file mode 100644 index 46d9c50..0000000 Binary files a/vignettes/figure/unnamed-chunk-14.png and /dev/null differ diff --git a/vignettes/figure/unnamed-chunk-15.png b/vignettes/figure/unnamed-chunk-15.png deleted file mode 100644 index 0eddbc6..0000000 Binary files a/vignettes/figure/unnamed-chunk-15.png and /dev/null differ diff --git a/vignettes/figure/unnamed-chunk-16.png b/vignettes/figure/unnamed-chunk-16.png deleted file mode 100644 index 4d1b738..0000000 Binary files a/vignettes/figure/unnamed-chunk-16.png and /dev/null differ diff --git a/vignettes/figure/unnamed-chunk-18.png b/vignettes/figure/unnamed-chunk-18.png deleted file mode 100644 index b109a11..0000000 Binary files a/vignettes/figure/unnamed-chunk-18.png and /dev/null differ diff --git a/vignettes/figure/unnamed-chunk-19.png b/vignettes/figure/unnamed-chunk-19.png deleted file mode 100644 index af84c2a..0000000 Binary files a/vignettes/figure/unnamed-chunk-19.png and /dev/null differ diff --git a/vignettes/figure/unnamed-chunk-4.png b/vignettes/figure/unnamed-chunk-4.png deleted file mode 100644 index 04187f8..0000000 Binary files a/vignettes/figure/unnamed-chunk-4.png and /dev/null differ diff --git a/vignettes/figure/unnamed-chunk-5.png b/vignettes/figure/unnamed-chunk-5.png deleted file mode 100644 index f40ba5c..0000000 Binary files a/vignettes/figure/unnamed-chunk-5.png and /dev/null differ diff --git a/vignettes/figure/unnamed-chunk-9.png b/vignettes/figure/unnamed-chunk-9.png deleted file mode 100644 index 76fd0c3..0000000 Binary files a/vignettes/figure/unnamed-chunk-9.png and /dev/null differ 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} - -<>= -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 deleted file mode 100644 index 42a44d6..0000000 Binary files a/vignettes/mkin.pdf and /dev/null differ 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/} -} -- cgit v1.2.1