From 8bf2a6f8f6bd752433f06b26b5da334e958f5166 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 13 Nov 2013 16:00:01 +0100 Subject: Update of README using TortoiseGit --- README.md | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/README.md b/README.md index 5e884870..4c0e745a 100644 --- a/README.md +++ b/README.md @@ -113,11 +113,10 @@ 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, I just (2013-11-11) noticed the github repositories -[StudyKin](http://github.com/zhenglei-gao/StudyKin) and -[KineticEval](http://github.com/zhenglei-gao/KineticEval), the latter of which appears to be -actively developed, so the different tools will hopefully be able to learn -from each other in the future as well. +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 -- cgit v1.2.1 From c99b5c298713a7c14e8ab5604c68613d0b7af27a Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 13 Nov 2013 16:25:31 +0100 Subject: Update of TODO from git bash on Windows --- TODO | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/TODO b/TODO index a56ca55f..54d135a5 100644 --- a/TODO +++ b/TODO @@ -2,13 +2,13 @@ TODO for version 1.0 - Transfer (some?) unit tests to examples to make them more easily accessible for the user - they will be tested during R CMD check anyway - Let mkin call mkin_wide_to_long automatically, if observed data are suitably defined +- Convert package vignettes to knitr and Rmd - Think about what a user would expect from version 1.0 -- Check the chi2 error calculation in mkinerrmin.R. KinGUII does this without iteration +- Move the GUI to a separate package as it has special dependencies and should have it own versioning Nice to have: - Improve choice of starting parameters, in order to avoid failure of eigenvalue based solutions due to singular matrix -- Consider reverting to deSolve in case eigenvalue based solution fails due to singular matrix - 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 -- cgit v1.2.1 From 8f1bae2142b37a0ff6b8989b2d1569686937f68e Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 13 Nov 2013 20:30:11 +0100 Subject: Add initial weighting choice to GUI --- inst/GUI/mkinGUI.R | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/inst/GUI/mkinGUI.R b/inst/GUI/mkinGUI.R index fb8f22ab..d6fb82e3 100644 --- a/inst/GUI/mkinGUI.R +++ b/inst/GUI/mkinGUI.R @@ -235,7 +235,8 @@ gbutton("Configure fit for selected model and dataset", cont = dsm, stmp <<- summary(ftmp) svalue(pf) <- paste0("Dataset ", ds.i, ", Model ", m[[m.i]]$name) show_plot("Initial", default = TRUE) - svalue(f.gg.opts.st) <<- "auto" + svalue(f.gg.opts.st) <<- ftmp$solution_type + svalue(f.gg.opts.weight) <<- "manual" f.gg.parms[,] <- get_Parameters(stmp, FALSE) svalue(f.gg.summary) <- capture.output(stmp) svalue(center) <- 3 @@ -568,7 +569,6 @@ show_plot <- function(type, default = FALSE) { ftmp$ds <<- ds[[ds.i]] } - #tmp <- get_tempfile(ext=".svg") svg(tf, width = 7, height = 5) plot(ftmp, main = ftmp$ds$title, xlab = ifelse(ftmp$ds$time_unit == "", "Time", @@ -608,18 +608,26 @@ run_fit <- function() { iniparms <- Parameters.ini$Initial names(iniparms) <- sub("_0", "", Parameters.ini$Name) inifixed <- names(iniparms[Parameters.ini$Fixed]) + weight <- svalue(f.gg.opts.weight) + if (weight == "manual") { + err = "err" + } else { + err = NULL + } ftmp <<- mkinfit(ftmp$mkinmod, override(ds[[ds.i]]$data), state.ini = iniparms, fixed_initials = inifixed, parms.ini = deparms, fixed_parms = defixed, solution_type = svalue(f.gg.opts.st), - err = "err") + weight = weight, + err = err) ftmp$ds.index <<- ds.i ftmp$ds <<- ds[[ds.i]] stmp <<- summary(ftmp) show_plot("Optimised") svalue(f.gg.opts.st) <- ftmp$solution_type + svalue(f.gg.opts.weight) <- ftmp$weight.ini f.gg.parms[,] <- get_Parameters(stmp, TRUE) svalue(f.gg.summary) <- capture.output(stmp) } @@ -646,9 +654,13 @@ dev.off() plot.gs <- gsvg(tf, container = f.gg.mid, width = 490, height = 350) f.gg.opts <- gformlayout(cont = f.gg.mid) solution_types <- c("auto", "analytical", "eigen", "deSolve") +weights <- c("manual", "none", "std", "mean") f.gg.opts.st <- gcombobox(solution_types, selected = 1, label = "solution_type", width = 200, cont = f.gg.opts) +f.gg.opts.weight <- gcombobox(weights, selected = 1, + label = "weight", width = 200, + cont = f.gg.opts) # Dataframe with initial and optimised parameters {{{3 f.gg.parms <- gdf(Parameters, width = 420, height = 300, cont = pf, @@ -717,6 +729,7 @@ update_plotting_and_fitting <- function() { ", Model ", ftmp$mkinmod$name) show_plot("Optimised") svalue(f.gg.opts.st) <- ftmp$solution_type + svalue(f.gg.opts.weight) <- ftmp$weight.ini f.gg.parms[,] <- get_Parameters(stmp, TRUE) svalue(f.gg.summary) <- capture.output(stmp) } -- cgit v1.2.1 From d8dbf2ad866fb9d34cd1100000b9c116219ecef6 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 13 Nov 2013 22:09:41 +0100 Subject: Add possibility to control more mkinfit arguments from the GUI --- R/mkinfit.R | 2 ++ inst/GUI/mkinGUI.R | 45 ++++++++++++++++++++++++++++++++++++++------- 2 files changed, 40 insertions(+), 7 deletions(-) diff --git a/R/mkinfit.R b/R/mkinfit.R index d04d1ea6..3e47479f 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -281,6 +281,8 @@ mkinfit <- function(mkinmod, observed, 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 diff --git a/inst/GUI/mkinGUI.R b/inst/GUI/mkinGUI.R index d6fb82e3..07f892e0 100644 --- a/inst/GUI/mkinGUI.R +++ b/inst/GUI/mkinGUI.R @@ -236,7 +236,14 @@ gbutton("Configure fit for selected model and dataset", cont = dsm, svalue(pf) <- paste0("Dataset ", ds.i, ", Model ", m[[m.i]]$name) show_plot("Initial", default = TRUE) svalue(f.gg.opts.st) <<- ftmp$solution_type - svalue(f.gg.opts.weight) <<- "manual" + svalue(f.gg.opts.weight) <<- ftmp$weight + svalue(f.gg.opts.atol) <<- ftmp$atol + svalue(f.gg.opts.rtol) <<- ftmp$rtol + svalue(f.gg.opts.reweight.method) <<- ifelse( + is.null(ftmp$reweight.method), + "none", ftmp$reweight.method) + svalue(f.gg.opts.reweight.tol) <<- ftmp$reweight.tol + svalue(f.gg.opts.reweight.max.iter) <<- ftmp$reweight.max.iter f.gg.parms[,] <- get_Parameters(stmp, FALSE) svalue(f.gg.summary) <- capture.output(stmp) svalue(center) <- 3 @@ -614,14 +621,22 @@ run_fit <- function() { } else { err = NULL } + reweight.method <- svalue(f.gg.opts.reweight.method) + if (reweight.method == "none") reweight.method = NULL ftmp <<- mkinfit(ftmp$mkinmod, override(ds[[ds.i]]$data), state.ini = iniparms, fixed_initials = inifixed, parms.ini = deparms, fixed_parms = defixed, solution_type = svalue(f.gg.opts.st), + atol = as.numeric(svalue(f.gg.opts.atol)), + rtol = as.numeric(svalue(f.gg.opts.rtol)), weight = weight, - err = err) + err = err, + reweight.method = reweight.method, + reweight.tol = as.numeric(svalue(f.gg.opts.reweight.tol)), + reweight.max.iter = as.numeric(svalue(f.gg.opts.reweight.max.iter)) + ) ftmp$ds.index <<- ds.i ftmp$ds <<- ds[[ds.i]] stmp <<- summary(ftmp) @@ -651,16 +666,27 @@ tf <- get_tempfile(ext=".svg") svg(tf, width = 7, height = 5) plot(ftmp) dev.off() -plot.gs <- gsvg(tf, container = f.gg.mid, width = 490, height = 350) +plot.gs <- gsvg(tf, container = f.gg.mid, width = 420, height = 300) f.gg.opts <- gformlayout(cont = f.gg.mid) solution_types <- c("auto", "analytical", "eigen", "deSolve") -weights <- c("manual", "none", "std", "mean") f.gg.opts.st <- gcombobox(solution_types, selected = 1, label = "solution_type", width = 200, cont = f.gg.opts) -f.gg.opts.weight <- gcombobox(weights, selected = 1, - label = "weight", width = 200, - cont = f.gg.opts) +f.gg.opts.atol <- gedit(ftmp$atol, label = "atol", width = 20, + cont = f.gg.opts) +f.gg.opts.rtol <- gedit(ftmp$rtol, label = "rtol", width = 20, + cont = f.gg.opts) +weights <- c("manual", "none", "std", "mean") +f.gg.opts.weight <- gcombobox(weights, selected = 1, label = "weight", + width = 200, cont = f.gg.opts) +f.gg.opts.reweight.method <- gcombobox(c("none", "obs"), selected = 1, + label = "reweight.method", + width = 200, + cont = f.gg.opts) +f.gg.opts.reweight.tol <- gedit(1e-8, label = "reweight.tol", + width = 20, cont = f.gg.opts) +f.gg.opts.reweight.max.iter <- gedit(10, label = "reweight.max.iter", + width = 20, cont = f.gg.opts) # Dataframe with initial and optimised parameters {{{3 f.gg.parms <- gdf(Parameters, width = 420, height = 300, cont = pf, @@ -730,6 +756,11 @@ update_plotting_and_fitting <- function() { show_plot("Optimised") svalue(f.gg.opts.st) <- ftmp$solution_type svalue(f.gg.opts.weight) <- ftmp$weight.ini + svalue(f.gg.opts.reweight.method) <- ifelse(is.null(ftmp$reweight.method), + "none", + ftmp$reweight.method) + svalue(f.gg.opts.reweight.tol) <- ftmp$reweight.tol + svalue(f.gg.opts.reweight.max.iter) <- ftmp$reweight.max.iter f.gg.parms[,] <- get_Parameters(stmp, TRUE) svalue(f.gg.summary) <- capture.output(stmp) } -- cgit v1.2.1 From ebc6f65e4c8b865fb9207ab11dc43cf4ac122c72 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Sun, 17 Nov 2013 16:13:13 +0100 Subject: Change vignette format to knitr (see ChangeLog) --- .Rbuildignore | 13 + .gitignore | 9 + ChangeLog | 5 + DESCRIPTION | 7 +- GNUmakefile | 20 +- README.md | 42 +- TODO | 8 +- vignettes/FOCUS_L.Rmd | 243 +++++++ vignettes/FOCUS_L.md | 931 +++++++++++++++++++++++++ vignettes/FOCUS_Z.Rnw | 261 +++++++ vignettes/FOCUS_Z.pdf | Bin 0 -> 204581 bytes vignettes/examples.Rnw | 524 -------------- vignettes/figure/FOCUS_2006_Z_fits_1.pdf | Bin 0 -> 5670 bytes vignettes/figure/FOCUS_2006_Z_fits_10.pdf | Bin 0 -> 7587 bytes vignettes/figure/FOCUS_2006_Z_fits_11.pdf | Bin 0 -> 7605 bytes vignettes/figure/FOCUS_2006_Z_fits_2.pdf | Bin 0 -> 5670 bytes vignettes/figure/FOCUS_2006_Z_fits_3.pdf | Bin 0 -> 5670 bytes vignettes/figure/FOCUS_2006_Z_fits_5.pdf | Bin 0 -> 6155 bytes vignettes/figure/FOCUS_2006_Z_fits_6.pdf | Bin 0 -> 6966 bytes vignettes/figure/FOCUS_2006_Z_fits_7.pdf | Bin 0 -> 6967 bytes vignettes/figure/FOCUS_2006_Z_fits_8.pdf | Bin 0 -> 6158 bytes vignettes/figure/FOCUS_2006_Z_fits_9.pdf | Bin 0 -> 6785 bytes vignettes/figure/FOCUS_2006_Z_residuals_11.pdf | Bin 0 -> 5943 bytes vignettes/figure/FOCUS_2006_Z_residuals_6.pdf | Bin 0 -> 5959 bytes vignettes/figure/unnamed-chunk-10.png | Bin 0 -> 8140 bytes vignettes/figure/unnamed-chunk-11.png | Bin 0 -> 5089 bytes vignettes/figure/unnamed-chunk-12.png | Bin 0 -> 5518 bytes vignettes/figure/unnamed-chunk-14.png | Bin 0 -> 5690 bytes vignettes/figure/unnamed-chunk-15.png | Bin 0 -> 5742 bytes vignettes/figure/unnamed-chunk-16.png | Bin 0 -> 5691 bytes vignettes/figure/unnamed-chunk-18.png | Bin 0 -> 5373 bytes vignettes/figure/unnamed-chunk-19.png | Bin 0 -> 5333 bytes vignettes/figure/unnamed-chunk-4.png | Bin 0 -> 5950 bytes vignettes/figure/unnamed-chunk-5.png | Bin 0 -> 4665 bytes vignettes/figure/unnamed-chunk-9.png | Bin 0 -> 7851 bytes vignettes/header.tex | 8 - vignettes/mkin.Rnw | 242 ++----- vignettes/mkin.pdf | Bin 0 -> 125145 bytes vignettes/references.bib | 9 +- 39 files changed, 1595 insertions(+), 727 deletions(-) create mode 100644 .gitignore create mode 100644 vignettes/FOCUS_L.Rmd create mode 100644 vignettes/FOCUS_L.md create mode 100644 vignettes/FOCUS_Z.Rnw create mode 100644 vignettes/FOCUS_Z.pdf delete mode 100644 vignettes/examples.Rnw create mode 100644 vignettes/figure/FOCUS_2006_Z_fits_1.pdf create mode 100644 vignettes/figure/FOCUS_2006_Z_fits_10.pdf create mode 100644 vignettes/figure/FOCUS_2006_Z_fits_11.pdf create mode 100644 vignettes/figure/FOCUS_2006_Z_fits_2.pdf create mode 100644 vignettes/figure/FOCUS_2006_Z_fits_3.pdf create mode 100644 vignettes/figure/FOCUS_2006_Z_fits_5.pdf create mode 100644 vignettes/figure/FOCUS_2006_Z_fits_6.pdf create mode 100644 vignettes/figure/FOCUS_2006_Z_fits_7.pdf create mode 100644 vignettes/figure/FOCUS_2006_Z_fits_8.pdf create mode 100644 vignettes/figure/FOCUS_2006_Z_fits_9.pdf create mode 100644 vignettes/figure/FOCUS_2006_Z_residuals_11.pdf create mode 100644 vignettes/figure/FOCUS_2006_Z_residuals_6.pdf create mode 100644 vignettes/figure/unnamed-chunk-10.png create mode 100644 vignettes/figure/unnamed-chunk-11.png create mode 100644 vignettes/figure/unnamed-chunk-12.png create mode 100644 vignettes/figure/unnamed-chunk-14.png create mode 100644 vignettes/figure/unnamed-chunk-15.png create mode 100644 vignettes/figure/unnamed-chunk-16.png create mode 100644 vignettes/figure/unnamed-chunk-18.png create mode 100644 vignettes/figure/unnamed-chunk-19.png create mode 100644 vignettes/figure/unnamed-chunk-4.png create mode 100644 vignettes/figure/unnamed-chunk-5.png create mode 100644 vignettes/figure/unnamed-chunk-9.png create mode 100644 vignettes/mkin.pdf diff --git a/.Rbuildignore b/.Rbuildignore index 59a46577..363190f4 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -1 +1,14 @@ GNUmakefile +out$ +toc$ +bbl$ +blg$ +aux$ +log$ +vignettes/figure +vignettes/mkin.tex +vignettes/mkin.pdf +vignettes/FOCUS_L.md +vignettes/FOCUS_L.html +vignettes/FOCUS_Z.tex +vignettes/FOCUS_Z.pdf diff --git a/.gitignore b/.gitignore new file mode 100644 index 00000000..1cec81aa --- /dev/null +++ b/.gitignore @@ -0,0 +1,9 @@ +vignettes/*.aux +vignettes/*.bbl +vignettes/*.blg +vignettes/*.log +vignettes/*.out +vignettes/*.toc +vignettes/mkin.tex +vignettes/FOCUS_L.html +vignettes/FOCUS_Z.tex diff --git a/ChangeLog b/ChangeLog index 852ad466..74c161fd 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,8 @@ +2013-11-17 Johannes Ranke for mkin (0.9-25) + + * Change vignette format from Sweave to knitr + * Split examples vignette to FOCUS_L and FOCUS_Z + 2013-11-06 Johannes Ranke for mkin (0.9-24) * Bugfix re-enabling the fixing of any combination of initial values diff --git a/DESCRIPTION b/DESCRIPTION index c8d39fff..80f87454 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -3,7 +3,7 @@ Type: Package Title: Routines for fitting kinetic models with one or more state variables to chemical degradation data Version: 0.9-25 -Date: 2013-11-06 +Date: 2013-11-17 Author: Johannes Ranke, with contributions from Katrin Lindenberger, René Lehmann Maintainer: Johannes Ranke Description: Calculation routines based on the FOCUS Kinetics Report (2006). @@ -13,9 +13,10 @@ Description: Calculation routines based on the FOCUS Kinetics Report (2006). (default is a Levenberg-Marquardt variant). Please note that no warranty is implied for correctness of results or fitness for a particular purpose. Depends: FME, deSolve, minpack.lm -Suggests: RUnit, gWidgetsWWW2, RSVGTipsDevice +Suggests: knitr, RUnit, gWidgetsWWW2, RSVGTipsDevice License: GPL LazyLoad: yes LazyData: yes Encoding: UTF-8 -URL: http://cran.r-project.org, http://kinfit.r-forge.r-project.org +VignetteBuilder: knitr +URL: http://cran.r-project.org, http://kinfit.r-forge.r-project.org, http://github.com/jranke/mkin diff --git a/GNUmakefile b/GNUmakefile index 8cdc60c0..1c124b19 100644 --- a/GNUmakefile +++ b/GNUmakefile @@ -8,7 +8,6 @@ PKGSRC := $(shell basename $(PWD)) # containing the first instance of R on the PATH. RBIN ?= $(shell dirname "`which R`") - .PHONY: help help: @@ -17,7 +16,7 @@ help: @echo "" @echo "Development Tasks" @echo "-----------------" - @echo " build Invoke docs and then create a package" + @echo " build Create the package" @echo " check Invoke build and then check the package" @echo " install Invoke build and then install the result" @echo " test Install a new copy of the package and run it " @@ -40,6 +39,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 @@ -56,10 +59,13 @@ test: install # Packaging Tasks #------------------------------------------------------------------------------ release: - @echo "\nPull in changes from svn and merge local commits" - @git svn rebase @echo "\nHow about make test and make check?" @echo "\nIs the DESCRIPTION file up to date?" - @echo "\nPerform final changes and commit with 'git commit --amend'." - @echo "\nIf the above is taken care of, run 'git svn dcommit'" - @echo "and then 'git push origin master'" + @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'" diff --git a/README.md b/README.md index 4c0e745a..e18ec238 100644 --- a/README.md +++ b/README.md @@ -2,7 +2,8 @@ The R package **mkin** provides calculation routines for the analysis of chemical degradation data, including **m**ulticompartment **kin**etics as needed for modelling -the formation and decline of transformation products. +the formation and decline of transformation products, or if several compartments +are involved. ## Installation @@ -27,11 +28,32 @@ require(devtools) install_github("mkin", "jranke") ``` +## 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 -For a start, have a look at the examples provided in the +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 +documentation or the package vignettes referenced from the [mkin package documentation page](http://kinfit.r-forge.r-project.org/mkin_static/index.html) @@ -66,12 +88,14 @@ or the package vignettes referenced from the * 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 + 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). + 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 @@ -86,6 +110,13 @@ 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 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. @@ -118,6 +149,7 @@ Finally, there is 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 diff --git a/TODO b/TODO index 54d135a5..bbdbac9b 100644 --- a/TODO +++ b/TODO @@ -1,15 +1,11 @@ TODO for version 1.0 -- Transfer (some?) unit tests to examples to make them more easily accessible - for the user - they will be tested during R CMD check anyway - Let mkin call mkin_wide_to_long automatically, if observed data are suitably defined -- Convert package vignettes to knitr and Rmd - Think about what a user would expect from version 1.0 - Move the GUI to a separate package as it has special dependencies and should have it own versioning +- Support model definitions without pathway to sink in combination with formation fractions +- Complete the main package vignette named mkin to include a method description Nice to have: -- Improve choice of starting parameters, in order to avoid failure of eigenvalue based solutions - due to singular matrix - 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 -- Document assumptions and implications of fitting method used in FME in the package vignette diff --git a/vignettes/FOCUS_L.Rmd b/vignettes/FOCUS_L.Rmd new file mode 100644 index 00000000..957b34ab --- /dev/null +++ b/vignettes/FOCUS_L.Rmd @@ -0,0 +1,243 @@ + + +# 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 new file mode 100644 index 00000000..6c43889d --- /dev/null +++ b/vignettes/FOCUS_L.md @@ -0,0 +1,931 @@ + + +# 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 new file mode 100644 index 00000000..44cfa468 --- /dev/null +++ b/vignettes/FOCUS_Z.Rnw @@ -0,0 +1,261 @@ +%\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 +Eurofins Regulatory AG\\ +Weidenweg 15, CH--4310 Rheinfelden, Switzerland\\[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}. This model definition is not supported when formation fractions +are used. + +<>= +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) +@ + +\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") +@ + +\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"), + 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. This model appears overparameterised (no +covariance matrix returned) if the sink for Z1 is left in the model. + +<>= +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) +@ + +Looking at the confidence intervals of the SFORB model parameters of Z3, 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 new file mode 100644 index 00000000..ca67191e Binary files /dev/null and b/vignettes/FOCUS_Z.pdf differ diff --git a/vignettes/examples.Rnw b/vignettes/examples.Rnw deleted file mode 100644 index f876d14a..00000000 --- a/vignettes/examples.Rnw +++ /dev/null @@ -1,524 +0,0 @@ -% $Id: examples.Rnw 66 2010-09-03 08:50:26Z jranke $ -%%\VignetteIndexEntry{Examples for kinetic evaluations using mkin} -%%VignetteDepends{FME} -%%\usepackage{Sweave} -\documentclass[12pt,a4paper]{article} -\usepackage{a4wide} -%%\usepackage[lists,heads]{endfloat} -\input{header} -\hypersetup{ - pdftitle = {Examples for kinetic evaluations using mkin}, - pdfsubject = {Manuscript}, - pdfauthor = {Johannes Ranke}, - colorlinks = {true}, - linkcolor = {blue}, - citecolor = {blue}, - urlcolor = {red}, - hyperindex = {true}, - linktocpage = {true}, -} -\SweaveOpts{engine=R, eps=FALSE, keep.source = TRUE} -<>= -options(prompt = "R> ") -options(width = 70) -options(SweaveHooks = list( - cex = function() par(cex.lab = 1.3, cex.axis = 1.3))) -@ -\begin{document} -\title{Examples for kinetic evaluations using mkin} -\author{\textbf{Johannes Ranke} \\[0.5cm] -%EndAName -Eurofins Regulatory AG\\ -Weidenweg 15, CH--4310 Rheinfelden, Switzerland\\[0.5cm] -and\\[0.5cm] -University of Bremen\\ -} -\maketitle - -%\begin{abstract} -%\end{abstract} - - -\thispagestyle{empty} \setcounter{page}{0} - -\clearpage - -\tableofcontents - - -\textbf{Key words}: Kinetics, FOCUS, nonlinear optimisation - -\section{Kinetic evaluations for parent compounds} - -These examples are also evaluated in a parallel vignette of the -\Rpackage{kinfit} package \citep{pkg:kinfit}. The datasets are from Appendix 3, -of the FOCUS kinetics report \citep{FOCUS2006, FOCUSkinetics2011}. - -\subsection{Laboratory Data L1} - -The following code defines example dataset L1 from the FOCUS kinetics -report, p. 284 - -<>= -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 \texttt{parent}. - -<>= -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. - -<>= -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. - -<>= -plot(m.L1.SFO) -@ - -The residual plot can be easily obtained by - -<>= -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. - -<>= -m.L1.FOMC <- mkinfit(FOMC, FOCUS_2006_L1_mkin, quiet=TRUE) -summary(m.L1.FOMC) -@ - -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. - -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 \texttt{mkin}. The reason for this is not known. However, -\texttt{mkin} gives the same $\chi^2$ error levels as the \Rpackage{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. - -\subsection{Laboratory Data L2} - -The following code defines example dataset L2 from the FOCUS kinetics -report, p. 287 - -<>= -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. - -<>= -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. - -<>= -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 \textit{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. - -<>= -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. - -<>= -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. - -<>= -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. - -\subsection{Laboratory Data L3} - -The following code defines example dataset L3 from the FOCUS kinetics report, -p. 290. - -<>= -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: - -<>= -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 22\% as well as the plot suggest that the model -does not fit very well. - -The FOMC model performs better: - -<>= -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: - -<>= -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. - -\subsection{Laboratory Data L4} - -The following code defines example dataset L4 from the FOCUS kinetics -report, p. 293 - -<>= -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: - -<>= -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 - -<>= -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. - -\section{Kinetic evaluations for parent and metabolites} - -\subsection{Laboratory Data for example compound Z} - -The following code defines the example dataset from Appendix 7 to the FOCUS kinetics -report, p.350 - -<>= -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) -@ - -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}. This model definition is not supported when formation fractions -are used. - -<>= -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, parms.ini = c(k_Z0_Z1 = 0.5), - quiet = TRUE) -#m.Z.3 <- mkinfit(Z.3, FOCUS_2006_Z_mkin, solution_type = "deSolve") -plot(m.Z.3) -summary(m.Z.3, data = FALSE) -@ - -The first attempt to fit the model failed, as the default solution type chosen -by mkinfit is based on eigenvalues, and the system defined by the starting -parameters is identified as being singular to the solver. This is caused by the -fact that the rate constants for both state variables are the same using the -default starting paramters. Setting a different starting value for one of the -parameters overcomes this. Alternatively, the \Rpackage{deSolve} based model -solution can be chosen, at the cost of a bit more computing time. - -<>= -Z.4a <- mkinmod(Z0 = list(type = "SFO", to = "Z1", sink = FALSE), - Z1 = list(type = "SFO", to = "Z2"), - Z2 = list(type = "SFO")) -m.Z.4a <- mkinfit(Z.4a, FOCUS_2006_Z_mkin, parms.ini = c(k_Z0_Z1 = 0.5), - quiet = TRUE) -plot(m.Z.4a) -summary(m.Z.4a, data = FALSE) -@ - -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. Again, in order to avoid a singular system when using default starting -parameters, the starting parameter for the pathway without sink term has to be adapted. - -<>= -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, - parms.ini = c(k_Z0_Z1 = 0.5, k_Z1_Z2 = 0.2), quiet = TRUE) -plot(m.Z.5) -summary(m.Z.5, data = FALSE) -@ - -Finally, metabolite Z3 is added to the model. - -<>= -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 = c(k_Z0_Z1 = 0.5, k_Z1_Z2 = 0.2, k_Z2_Z3 = 0.3), - 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") -@ - -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"), - 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. This model appears overparameterised (no -covariance matrix returned) if the sink for Z1 is left in the model. - -<>= -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. - -Using the SFORB additionally for Z1 or Z2 did not further improve the result. - -<>= -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 = c(k_Z1_Z2 = 0.2), quiet = TRUE) -plot(m.Z.mkin.5) -summary(m.Z.mkin.5, data = FALSE) -@ - -Looking at the confidence intervals of the SFORB model parameters of Z3, 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/figure/FOCUS_2006_Z_fits_1.pdf b/vignettes/figure/FOCUS_2006_Z_fits_1.pdf new file mode 100644 index 00000000..9ee45463 Binary files /dev/null and b/vignettes/figure/FOCUS_2006_Z_fits_1.pdf differ diff --git a/vignettes/figure/FOCUS_2006_Z_fits_10.pdf b/vignettes/figure/FOCUS_2006_Z_fits_10.pdf new file mode 100644 index 00000000..c2776cab Binary files /dev/null and b/vignettes/figure/FOCUS_2006_Z_fits_10.pdf differ diff --git a/vignettes/figure/FOCUS_2006_Z_fits_11.pdf b/vignettes/figure/FOCUS_2006_Z_fits_11.pdf new file mode 100644 index 00000000..2b1c973c Binary files /dev/null and b/vignettes/figure/FOCUS_2006_Z_fits_11.pdf differ diff --git a/vignettes/figure/FOCUS_2006_Z_fits_2.pdf b/vignettes/figure/FOCUS_2006_Z_fits_2.pdf new file mode 100644 index 00000000..511fe81c Binary files /dev/null and b/vignettes/figure/FOCUS_2006_Z_fits_2.pdf differ diff --git a/vignettes/figure/FOCUS_2006_Z_fits_3.pdf b/vignettes/figure/FOCUS_2006_Z_fits_3.pdf new file mode 100644 index 00000000..e5472f9a Binary files /dev/null and b/vignettes/figure/FOCUS_2006_Z_fits_3.pdf differ diff --git a/vignettes/figure/FOCUS_2006_Z_fits_5.pdf b/vignettes/figure/FOCUS_2006_Z_fits_5.pdf new file mode 100644 index 00000000..831b7f26 Binary files /dev/null and b/vignettes/figure/FOCUS_2006_Z_fits_5.pdf differ diff --git a/vignettes/figure/FOCUS_2006_Z_fits_6.pdf b/vignettes/figure/FOCUS_2006_Z_fits_6.pdf new file mode 100644 index 00000000..6c4bcd33 Binary files /dev/null and b/vignettes/figure/FOCUS_2006_Z_fits_6.pdf differ diff --git a/vignettes/figure/FOCUS_2006_Z_fits_7.pdf b/vignettes/figure/FOCUS_2006_Z_fits_7.pdf new file mode 100644 index 00000000..e4ac1a83 Binary files /dev/null and b/vignettes/figure/FOCUS_2006_Z_fits_7.pdf differ diff --git a/vignettes/figure/FOCUS_2006_Z_fits_8.pdf b/vignettes/figure/FOCUS_2006_Z_fits_8.pdf new file mode 100644 index 00000000..11303976 Binary files /dev/null and b/vignettes/figure/FOCUS_2006_Z_fits_8.pdf differ diff --git a/vignettes/figure/FOCUS_2006_Z_fits_9.pdf b/vignettes/figure/FOCUS_2006_Z_fits_9.pdf new file mode 100644 index 00000000..920353a7 Binary files /dev/null and b/vignettes/figure/FOCUS_2006_Z_fits_9.pdf differ diff --git a/vignettes/figure/FOCUS_2006_Z_residuals_11.pdf b/vignettes/figure/FOCUS_2006_Z_residuals_11.pdf new file mode 100644 index 00000000..40a6afb1 Binary files /dev/null and b/vignettes/figure/FOCUS_2006_Z_residuals_11.pdf differ diff --git a/vignettes/figure/FOCUS_2006_Z_residuals_6.pdf b/vignettes/figure/FOCUS_2006_Z_residuals_6.pdf new file mode 100644 index 00000000..a3006d8f Binary files /dev/null and b/vignettes/figure/FOCUS_2006_Z_residuals_6.pdf differ diff --git a/vignettes/figure/unnamed-chunk-10.png b/vignettes/figure/unnamed-chunk-10.png new file mode 100644 index 00000000..cd3a700e Binary files /dev/null and b/vignettes/figure/unnamed-chunk-10.png differ diff --git a/vignettes/figure/unnamed-chunk-11.png b/vignettes/figure/unnamed-chunk-11.png new file mode 100644 index 00000000..ca488e63 Binary files /dev/null and b/vignettes/figure/unnamed-chunk-11.png differ diff --git a/vignettes/figure/unnamed-chunk-12.png b/vignettes/figure/unnamed-chunk-12.png new file mode 100644 index 00000000..3a644136 Binary files /dev/null and b/vignettes/figure/unnamed-chunk-12.png differ diff --git a/vignettes/figure/unnamed-chunk-14.png b/vignettes/figure/unnamed-chunk-14.png new file mode 100644 index 00000000..46d9c50d Binary files /dev/null and b/vignettes/figure/unnamed-chunk-14.png differ diff --git a/vignettes/figure/unnamed-chunk-15.png b/vignettes/figure/unnamed-chunk-15.png new file mode 100644 index 00000000..0eddbc63 Binary files /dev/null and b/vignettes/figure/unnamed-chunk-15.png differ diff --git a/vignettes/figure/unnamed-chunk-16.png b/vignettes/figure/unnamed-chunk-16.png new file mode 100644 index 00000000..4d1b738a Binary files /dev/null and b/vignettes/figure/unnamed-chunk-16.png differ diff --git a/vignettes/figure/unnamed-chunk-18.png b/vignettes/figure/unnamed-chunk-18.png new file mode 100644 index 00000000..b109a11f Binary files /dev/null and b/vignettes/figure/unnamed-chunk-18.png differ diff --git a/vignettes/figure/unnamed-chunk-19.png b/vignettes/figure/unnamed-chunk-19.png new file mode 100644 index 00000000..af84c2a2 Binary files /dev/null and b/vignettes/figure/unnamed-chunk-19.png differ diff --git a/vignettes/figure/unnamed-chunk-4.png b/vignettes/figure/unnamed-chunk-4.png new file mode 100644 index 00000000..04187f8f Binary files /dev/null and b/vignettes/figure/unnamed-chunk-4.png differ diff --git a/vignettes/figure/unnamed-chunk-5.png b/vignettes/figure/unnamed-chunk-5.png new file mode 100644 index 00000000..f40ba5ce Binary files /dev/null and b/vignettes/figure/unnamed-chunk-5.png differ diff --git a/vignettes/figure/unnamed-chunk-9.png b/vignettes/figure/unnamed-chunk-9.png new file mode 100644 index 00000000..76fd0c33 Binary files /dev/null and b/vignettes/figure/unnamed-chunk-9.png differ diff --git a/vignettes/header.tex b/vignettes/header.tex index 707997c4..476415e5 100644 --- a/vignettes/header.tex +++ b/vignettes/header.tex @@ -21,11 +21,3 @@ \RequirePackage{graphicx,ae,fancyvrb} \IfFileExists{upquote.sty}{\RequirePackage{upquote}}{} \usepackage{relsize} - -\DefineVerbatimEnvironment{Sinput}{Verbatim}{baselinestretch=1.05} -\DefineVerbatimEnvironment{Soutput}{Verbatim}{fontfamily=courier, - baselinestretch=1.05, - fontshape=it, - fontsize=\relsize{-1}} -\DefineVerbatimEnvironment{Scode}{Verbatim}{} -\newenvironment{Schunk}{}{} diff --git a/vignettes/mkin.Rnw b/vignettes/mkin.Rnw index 5c71e8b7..398b5418 100644 --- a/vignettes/mkin.Rnw +++ b/vignettes/mkin.Rnw @@ -1,168 +1,74 @@ -% $Id: mkin.Rnw 66 2010-09-03 08:50:26Z jranke $ -%%\VignetteIndexEntry{Routines for fitting kinetic models with one or more state variables to chemical degradation data} -%%VignetteDepends{FME} -%%\usepackage{Sweave} -\documentclass[12pt,a4paper]{article} -\usepackage{a4wide} -%%\usepackage[lists,heads]{endfloat} -\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}, - hyperindex = {true}, - linktocpage = {true}, -} -\SweaveOpts{engine=R, eps=FALSE, keep.source = TRUE} -<>= -options(prompt = "R> ") -options(SweaveHooks = list( - cex = function() par(cex.lab = 1.3, cex.axis = 1.3))) -@ -\begin{document} -\title{mkin -\\ -Routines for fitting kinetic models with one or more state variables to chemical degradation data} -\author{\textbf{Johannes Ranke} \\[0.5cm] -%EndAName -Eurofins Regulatory AG\\ -Weidenweg 15, CH--4310 Rheinfelden, Switzerland\\[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. - -\section{Example} -\label{exam} - -In the following, requirements for data formatting are explained. Then the -procedure for fitting the four kinetic models recommended by the FOCUS group -to an example dataset for parent only given in the FOCUS kinetics report is -illustrated. The explanations are kept rather verbose in order to lower the -barrier for \RR{} newcomers. - -\subsection{Data format} - -The following listing shows example dataset C from the FOCUS kinetics -report as distributed with the \Rpackage{mkin} package - -<>= -library("mkin") -FOCUS_2006_C -@ - -Note that the data needs to be in the format of a data frame containing a -variable \Robject{name} specifying the observed variable, indicating the -compound name and, if applicable, the compartment, a variable \Robject{time} -containing sampling times, and a numeric variable \Robject{value} specifying -the observed value of the variable. If a further variable \Robject{error} -is present, this will be used to give different weights to the data points -(the higher the error, the lower the weight, see the help page of the -\Robject{modCost} function of the \Rpackage{FME} package \citep{soetaert10}). -Replicate measurements are not recorded in extra columns but simply appended, -leading to multiple occurrences of the sampling times \Robject{time}. - -Small to medium size dataset can be conveniently entered directly as \RR{} code -as shown in the following listing - -<>= -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) -) -@ - -\subsection{Model definition} - -The next task is to define the model to be fitted to the data. In order to -facilitate this task, a convenience function \Robject{mkinmod} is available. - -<>= -SFO <- mkinmod(parent = list(type = "SFO")) -SFORB <- mkinmod(parent = list(type = "SFORB")) -SFO_SFO <- mkinmod( - parent = list(type = "SFO", to = "m1", sink = TRUE), - m1 = list(type = "SFO")) -SFORB_SFO <- mkinmod( - parent = list(type = "SFORB", to = "m1", sink = TRUE), - m1 = list(type = "SFO")) -@ - -The model definitions given above define sets of linear first-order ordinary -differential equations. In these cases, a coefficient matrix is also returned. - -Other models that include time on the right-hand side of the differential -equation are the first-order multi-compartment (FOMC) model and the -Hockey-Stick (HS) model. At present, these models can only be used only for the -parent compound. - -\subsection{Fitting the model} - -Then the model parameters should be fitted to the data. The function -\Robject{mkinfit} internally creates a cost function using \Robject{modCost} -from the \Rpackage{FME} package and then produces a fit using \Robject{modFit} -from the same package. In cases of linear first-order differential -equations, the solution used for calculating the cost function is based -on the fundamental system of the coefficient matrix, as proposed by -\citet{bates88}. - -<>= -SFO.fit <- mkinfit(SFO, FOCUS_2006_C) -summary(SFO.fit) -SFORB.fit <- mkinfit(SFORB, FOCUS_2006_C) -summary(SFORB.fit) -SFO_SFO.fit <- mkinfit(SFO_SFO, FOCUS_2006_D) -summary(SFO_SFO.fit, data=FALSE) -SFORB_SFO.fit <- mkinfit(SFORB_SFO, FOCUS_2006_D) -summary(SFORB_SFO.fit, data=FALSE) -@ - -\section{Acknowledgements} - -This package would 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). Parts of the package were written during my employment at Harlan. - -\bibliographystyle{plainnat} -\bibliography{references} - -\end{document} -% vim: set foldmethod=syntax: +%\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 +Eurofins Regulatory AG\\ +Weidenweg 15, CH--4310 Rheinfelden, Switzerland\\[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 new file mode 100644 index 00000000..6849fb48 Binary files /dev/null and b/vignettes/mkin.pdf differ diff --git a/vignettes/references.bib b/vignettes/references.bib index 7796000f..e069daf2 100644 --- a/vignettes/references.bib +++ b/vignettes/references.bib @@ -1,4 +1,4 @@ -% This file was created with JabRef 2.7b. +% This file was originally created with JabRef 2.7b. % Encoding: ISO8859_1 @BOOK{bates88, @@ -16,8 +16,6 @@ month = {November}, year = {2011}, file = {FOCUS kinetics 2011 Generic guidance:/home/ranke/dok/orgs/focus/FOCUSkineticsvc_1_0_Nov23.pdf:PDF}, - owner = {ranke}, - timestamp = {2012.09.20}, url = {http://focus.jrc.ec.europa.eu/dk} } @@ -46,7 +44,7 @@ degradation data}, author = {Johannes Ranke}, year = {2013}, - url = {http://CRAN.R-project.org} + url = {http://CRAN.R-project.org/package=kinfit} } @MANUAL{pkg:mkin, @@ -54,7 +52,7 @@ variables to chemical degradation data}, author = {Johannes Ranke}, year = {2013}, - url = {http://CRAN.R-project.org} + url = {http://CRAN.R-project.org/package/kinfit} } @INPROCEEDINGS{schaefer2007, @@ -79,4 +77,3 @@ number = {3}, url = {http://www.jstatsoft.org/v33/i03/} } - -- cgit v1.2.1