From 4110808d1953b8704ef48d0dc5b5e445cbde9af0 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Tue, 20 May 2014 08:36:42 +0200 Subject: Step back regarding confidence intervals of formation fractions - Don't give confidence intervals for formation fractions that belong to a set, only for single formation fractions - Update the vignette with history and background - Correct 'isotropic' into 'isometric' in the ilr documentation, the README and the comments in transform_odeparms - Candidate for release to CRAN --- ChangeLog | 15 +++++++- DESCRIPTION | 2 +- R/mkinfit.R | 18 +++++++-- R/transform_odeparms.R | 4 +- README.md | 2 +- man/ilr.Rd | 6 +-- vignettes/header.tex | 1 - vignettes/mkin.Rnw | 95 +++++++++++++++++++++++++++++++++++++++++++++-- vignettes/mkin.pdf | Bin 124442 -> 160326 bytes vignettes/references.bib | 35 ++++++++++------- 10 files changed, 147 insertions(+), 31 deletions(-) diff --git a/ChangeLog b/ChangeLog index c60caf5e..2ce73e5c 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,12 +1,23 @@ -2014-05-07 Johannes Ranke for mkin (0.9-26) +2014-05-20 Johannes Ranke for mkin (0.9-28) + + * Do not backtransform confidence intervals for formation fractions if more + than one compound is formed, as such parameters only define the pathways + as a set + * Add historical remarks and some background to the main package vignette + * Correct 'isotropic' into 'isometric' for the ilr transformation + +2014-05-09 Johannes Ranke for mkin (0.9-27) * Fork the GUI into a separate package gmkin * DESCRIPTION, NAMESPACE, TODO: Adapt and add copyright information * Remove files belonging to the GUI - * Possibility to fit without parameter transformation + * Possibility to fit without parameter transformations, using bounds + as implemented in FME * Add McCall 2,4,5-T dataset * Enable selection of observed variables in plotting * Add possibility to show residual plot in plot.mkinfit + * R/mkinparplot.R, man/mkinparplot.Rd: plot parameters with confidence + intervals 2014-02-27 Johannes Ranke for mkin (0.9-25) diff --git a/DESCRIPTION b/DESCRIPTION index fc92b968..702d3604 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-28 -Date: 2014-05-15 +Date: 2014-05-20 Authors@R: c(person("Johannes", "Ranke", role = c("aut", "cre", "cph"), email = "jranke@uni-bremen.de"), person("Katrin", "Lindenberger", role = "ctb"), diff --git a/R/mkinfit.R b/R/mkinfit.R index ab1977d1..612e43d2 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -364,15 +364,27 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, "t value", "Pr(>|t|)", "Pr(>t)")) blci <- buci <- numeric() - # Only transform boundaries of CI for one parameter at a time + # Transform boundaries of CI for one parameter at a time, + # with the exception of sets of formation fractions (single fractions are OK). + f_names_skip <- character(0) + for (box in mod_vars) { # Figure out sets of fractions to skip + f_names <- grep(paste("^f", box, sep = "_"), pnames, value = TRUE) + n_paths <- length(f_names) + if (n_paths > 1) f_names_skip <- c(f_names_skip, f_names) + } + for (pname in pnames) { par.lower <- par.upper <- object$par par.lower[pname] <- param[pname, "Lower"] par.upper[pname] <- param[pname, "Upper"] - blci[pname] <- backtransform_odeparms(par.lower, mod_vars, + if (!pname %in% f_names_skip) { + blci[pname] <- backtransform_odeparms(par.lower, mod_vars, object$transform_rates, object$transform_fractions)[pname] - buci[pname] <- backtransform_odeparms(par.upper, mod_vars, + buci[pname] <- backtransform_odeparms(par.upper, mod_vars, object$transform_rates, object$transform_fractions)[pname] + } else { + blci[pname] <- buci[pname] <- NA + } } bparam <- cbind(object$bparms.optim, blci, buci) dimnames(bparam) <- list(pnames, c("Estimate", "Lower", "Upper")) diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R index a0c302f5..2438151b 100644 --- a/R/transform_odeparms.R +++ b/R/transform_odeparms.R @@ -30,7 +30,7 @@ transform_odeparms <- function(parms, mod_vars, else transparms[index_k] <- parms[index_k] } - # Go through state variables and apply isotropic logratio transformation if requested + # Go through state variables and apply isometric logratio transformation if requested for (box in mod_vars) { indices_f <- grep(paste("^f", box, sep = "_"), names(parms)) f_names <- grep(paste("^f", box, sep = "_"), names(parms), value = TRUE) @@ -73,7 +73,7 @@ backtransform_odeparms <- function(transparms, mod_vars, else parms[index_k] <- transparms[index_k] } - # Go through state variables and apply inverse isotropic logratio transformation + # Go through state variables and apply inverse isometric logratio transformation for (box in mod_vars) { indices_f <- grep(paste("^f", box, sep = "_"), names(transparms)) f_names <- grep(paste("^f", box, sep = "_"), names(transparms), value = TRUE) diff --git a/README.md b/README.md index 6b2fbede..9b510b0a 100644 --- a/README.md +++ b/README.md @@ -79,7 +79,7 @@ or the package vignettes referenced from the so their estimators can more reasonably be expected to follow a normal distribution. This has the side effect that no constraints are needed in the optimisation. Thanks to René Lehmann for the nice - cooperation on this, especially the isotropic logration transformation + cooperation on this, especially the isometric logration transformation that is now used for the formation fractions. * A side effect of this is that when parameter estimates are backtransformed to match the model definition, confidence intervals calculated from diff --git a/man/ilr.Rd b/man/ilr.Rd index cedb49ca..9d1483f3 100644 --- a/man/ilr.Rd +++ b/man/ilr.Rd @@ -2,10 +2,10 @@ \alias{ilr} \alias{invilr} \title{ - Function to perform isotropic log-ratio transformation + Function to perform isometric log-ratio transformation } \description{ - This implementation is a special case of the class of isotropic log-ratio transformations. + This implementation is a special case of the class of isometric log-ratio transformations. } \usage{ ilr(x) @@ -28,7 +28,7 @@ René Lehmann and Johannes Ranke } \seealso{ - Other implementations are in R packages \code{compositions} and \code{robCompositions}. + Another implementation can be found in R package \code{robCompositions}. } \examples{ # Order matters diff --git a/vignettes/header.tex b/vignettes/header.tex index 476415e5..e0775d34 100644 --- a/vignettes/header.tex +++ b/vignettes/header.tex @@ -8,7 +8,6 @@ \usepackage[round]{natbib} \usepackage{amstext} \usepackage{hyperref} -\usepackage[utf8]{inputenc} \newcommand{\Rpackage}[1]{{\normalfont\fontseries{b}\selectfont #1}} \newcommand{\Robject}[1]{\texttt{#1}} diff --git a/vignettes/mkin.Rnw b/vignettes/mkin.Rnw index 0ac114ed..1befe009 100644 --- a/vignettes/mkin.Rnw +++ b/vignettes/mkin.Rnw @@ -2,6 +2,7 @@ %\VignetteEngine{knitr::knitr} \documentclass[12pt,a4paper]{article} \usepackage{a4wide} +\usepackage[utf8]{inputenc} \input{header} \hypersetup{ pdftitle = {mkin - Routines for fitting kinetic models with one or more state variables to chemical degradation data}, @@ -29,7 +30,7 @@ chemical degradation data} Wissenschaftlicher Berater\\ Kronacher Str. 8, 79639 Grenzach-Wyhlen, Germany\\[0.5cm] and\\[0.5cm] -University of Bremen\\ +Privatdozent at the University of Bremen\\ } \maketitle @@ -58,14 +59,100 @@ metabolites. 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 +\citep{rcore2014} 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. +with transformation products, commonly termed metabolites, and to more than one +compartment. It is also possible to include back reactions, so equilibrium reactions +and equilibrium partitioning can be specified, although this oftentimes leads +to an overparameterisation of the model. + +When mkin was first published in May 2010, the most commonly used tools +for fitting more complex kinetic degradation models to experimental data were KinGUI +\citep{schaefer2007}, a MATLAB$^\circledR$ based tool with a graphical user +interface that was specifically tailored to the task and included some output +as proposed by the FOCUS Kinetics Workgroup, and ModelMaker, a general purpose +compartment based tool providing infrastructure for fitting dynamic simulation +models based on differential equations to data. + +At that time, the R package \Rpackage{FME} (Flexible Modelling Environment) +\citep{soetaert2010} was already available, and provided a good basis for +developing a package specifically tailored to the task. The remaining challenge +was to make it as easy as possible for the users (including the author of this +vignette) to specify the system of differential equations and to include the +output requested by the FOCUS guidance, such as the relative standard deviation +that has to be assumed for the residuals, such that the $\chi^2$ +goodness-of-fit test as defined by the FOCUS kinetics workgroup would pass +using an significance level $\alpha$ of 0.05. + +Also, mkin introduced using analytical solutions for parent only kinetics for +improved optimization speed. Later, Eigenvalue based solutions were +introduced to mkin for the case of linear differential equations (\textit{i.e.} +where the FOMC or DFOP models were not used for the parent compound), greatly +improving the optimization speed for these cases. + +Soon after the publication of mkin, two derived tools were published, namely +KinGUII (available from Bayer Crop Science) and CAKE (commissioned to Tessella +by Syngenta), which added a graphical user interface (GUI), and added fitting by +iteratively reweighted least squares (IRLS) and characterisation of likely +parameter distributions by Markov Chain Monte Carlo (MCMC) sampling. + +CAKE focuses on a smooth use experience, sacrificing some flexibility in the model +definition, allowing only two primary metabolites in parallel. KinGUI offers +quite a flexible widget for specifying complex kinetic models. Back-reactions +(non-instanteneous equilibria) were not present in the first version of +KinGUII, and only simple first-order models could be specified for +transformation products. As of May 2014 (KinGUII version 2.1), back-reactions +and biphasic modelling of metabolites are also available in KinGUII. + +Currently (May 2014), the main feature available in \Rpackage{mkin} which is +not present in KinGUII or CAKE, is the estimation of parameter confidence +intervals based on transformed parameters. For rate constants, the log +transformation is used, as proposed by Bates and Watts \citep[p. 77, p. +149]{bates1988}. Approximate intervals are constructed for the transformed rate +constants \citep[compare][p. 153]{bates1988}, \textit{i.e.} for their logarithms. +Confidence intervals for the rate constants are then obtained using the +appropriate backtransformation using the exponential function. + +In the first version of \Rpackage{mkin} allowing for specifying models using +formation fractions, a home-made reparameterisation was used in order to ensure +that the sum of formation fractions would not exceed unity. + +This method is still used in the current version of KinGUII (v2.1), with a +modification that allows for fixing the pathway to sink to zero. CAKE uses +penalties in the objective function in order to enforce this constraint. + +In 2012, an alternative reparameterisation of the formation fractions was +proposed together with René Lehmann \citep{ranke2012}, based on isometric +logratio transformation (ILR). The aim was to improve the validity of the +linear approximation of the objective function during the parameter +estimation procedure as well as in the subsequent calculation of parameter +confidence intervals. + +In the first attempt at providing improved parameter confidence intervals +introduced to \Rpackage{mkin} in 2013, confidence intervals obtained from +FME on the transformed parameters were simply all backtransformed one by one +to yield asymetric confidence intervals for the backtransformed parameters. + +However, while there is a 1:1 relation between the rate constants in the model +and the transformed parameters fitted in the model, the parameters obtained by the +isometric logratio transformation are calculated from the set of formation +fractions that quantify the paths to each of the compounds formed from a +specific parent compound, and no such 1:1 relation exists. + +Therefore, parameter confidence intervals for formation fractions obtained with +this method only appear valid for the case of a single transformation product, where +only one formation fraction is to be estimated, directly corresponding to one +component of the ilr transformed parameter. + +The confidence intervals obtained by backtransformation for the cases where a +1:1 relation between transformed and original parameter exist are considered by +the author of this vignette to be more accurate than those obtained using a +re-estimation of the Hessian matrix after backtransformation, as implemented +in the FME package. \bibliographystyle{plainnat} \bibliography{references} diff --git a/vignettes/mkin.pdf b/vignettes/mkin.pdf index 0e4d1be8..12293506 100644 Binary files a/vignettes/mkin.pdf and b/vignettes/mkin.pdf differ diff --git a/vignettes/references.bib b/vignettes/references.bib index e069daf2..88caea60 100644 --- a/vignettes/references.bib +++ b/vignettes/references.bib @@ -1,7 +1,7 @@ % This file was originally created with JabRef 2.7b. % Encoding: ISO8859_1 -@BOOK{bates88, +@BOOK{bates1988, title = {Nonlinear regression and its applications}, publisher = {Wiley-Interscience}, year = {1988}, @@ -29,12 +29,12 @@ url = {http://focus.jrc.ec.europa.eu/dk} } -@MANUAL{rcore2013, +@MANUAL{rcore2014, title = {\textsf{R}: A Language and Environment for Statistical Computing}, author = {{R Development Core Team}}, organization = {R Foundation for Statistical Computing}, address = {Vienna, Austria}, - year = {2013}, + year = {2014}, note = {{ISBN} 3-900051-07-0}, url = {http://www.R-project.org} } @@ -52,21 +52,20 @@ variables to chemical degradation data}, author = {Johannes Ranke}, year = {2013}, - url = {http://CRAN.R-project.org/package/kinfit} + url = {http://CRAN.R-project.org/package=mkin} } -@INPROCEEDINGS{schaefer2007, - author = {D. Sch\"{a}fer and M. Mikolasch and P. Rainbird and B. Harvey}, - title = {{KinGUI}: a new kinetic software tool for evaluations according to - FOCUS degradation kinetics}, - booktitle = {Proceedings of the XIII Symposium Pesticide Chemistry}, - year = {2007}, - editor = {Del Re A. A. M. and Capri E. and Fragoulis G. and Trevisan M.}, - pages = {916--923}, - address = {Piacenza} +@Inproceedings{ schaefer2007, + title = {{KinGUI}: a new kinetic software tool for evaluations according to {FOCUS} degradation kinetics}, + author = {D. Sch\"{a}fer and B. Mikolasch and P. Rainbird and B. Harvey}, + booktitle = {Proceedings of the XIII Symposium Pesticide Chemistry}, + editor = {Del Re A. A. M. and Capri E. and Fragoulis G. and Trevisan M.}, + year = {2007}, + address = {Piacenza}, + pages = {916--923} } -@ARTICLE{soetaert10, +@ARTICLE{soetaert2010, author = {Karline Soetaert and Thomas Petzoldt}, title = {Inverse Modelling, Sensitivity and Monte Carlo Analysis in {R} Using Package {FME}}, @@ -77,3 +76,11 @@ number = {3}, url = {http://www.jstatsoft.org/v33/i03/} } + +@Inproceedings{ ranke2012, + title = {Parameter reliability in kinetic evaluation of environmental metabolism data - Assessment and the influence of model specification}, + author = {J. Ranke and R. Lehmann}, + booktitle = {SETAC World 20-24 May}, + year = {2012}, + address = {Berlin} +} -- cgit v1.2.1