From 12934520d7ec3218ce1505787b6066334a24a562 Mon Sep 17 00:00:00 2001 From: jranke Date: Tue, 30 Mar 2010 19:49:44 +0000 Subject: Initial import of the kinfit package developed from 2008-07 to 2010-03 at Harlan Laboratories Ltd (former RCC Ltd). Supports fitting of parent data with the usual FOCUS kinetic models. git-svn-id: svn+ssh://svn.r-forge.r-project.org/svnroot/kinfit/pkg/kinfit@2 edb9625f-4e0d-4859-8d74-9fd3b1da38cb --- inst/doc/kinfit.Rnw | 933 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 933 insertions(+) create mode 100644 inst/doc/kinfit.Rnw (limited to 'inst/doc/kinfit.Rnw') diff --git a/inst/doc/kinfit.Rnw b/inst/doc/kinfit.Rnw new file mode 100644 index 0000000..7a590d5 --- /dev/null +++ b/inst/doc/kinfit.Rnw @@ -0,0 +1,933 @@ +%%\VignetteIndexEntry{Routines for fitting kinetic models to chemical degradation data} +%%VignetteDepends{nls} +%%\usepackage{Sweave} +\documentclass[12pt,a4paper]{article} +\usepackage{a4wide} +%%\usepackage[lists,heads]{endfloat} +\input{header} +\hypersetup{ + pdftitle = {kinfit - Routines for fitting kinetic models 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{kinfit -\\ +Routines for fitting kinetic models to chemical degradation +data\footnote{This is a preprint of an article published in +The R Journal, Volume 1, Number ?, ?--?. available online +\url{http://journal.r-project.org}.}} +\author{\textbf{Johannes Ranke} \\ +%EndAName +Product Safety \\ +Harlan Laboratories Ltd. \\ +Zelgliweg 1, CH--4452 Itingen, Switzerland} +\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 regression. +The \RR{} add-on package \Rpackage{kinfit} implements fitting the models +recommended in this guidance from within R and calculates the +recommended statistical measures for data series within one compartment +without metabolite data. +\end{abstract} + + +\thispagestyle{empty} \setcounter{page}{0} + +\clearpage + +\tableofcontents + +\textbf{Key words}: Kinetics, FOCUS, nonlinear fitting + +\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{rcore2009} implements the approach recommended in the kinetics report +provided by the FOrum for Co-ordination of pesticide fate models and their +USe \citep{FOCUS2006} for simple data series for one parent compound in one +compartment. + +\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 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{kinfit} package + +<>= +library("kinfit") +data("FOCUS_2006_C", package = "kinfit") +print(FOCUS_2006_C) +@ + +Note that the data needs to be in the format of a data frame containing +a variable \Robject{t} containing sampling times and a variable +\Robject{parent} containing the measured data. Replicate measurements are +not recorded in extra columns but simply appended, leading to multiple +occurrences of the sampling times \Robject{t}. + +Small to medium size dataset can be conveniently entered directly as \RR{} code +as shown in the following listing + +<>= +kindata_example <- data.frame( + t = c(0, 1, 3, 7, 14, 28, 63, 91, 119), + parent = c(85.1, 57.9, 29.9, 14.6, 9.7, 6.6, 4, 3.9, 0.6) +) +@ + + +\subsection{Fitting the kinetic models} + +The user can choose for which kinetic models the \Robject{kinfit} function +will try to find optimised parameters. This is achieved by the argument +\Robject{kinmodels} to the function, as shown below. The models currently +implemented are abbreviated \Robject{SFO} (Single First-Order), \Robject{FOMC} +(First-Order Multi-Compartment), \Robject{DFOP} (Double First-Order in Parallel) +and \Robject{HS} (Hockey-Stick) as defined by the \cite{FOCUS2006}. From the +DFOP model, corresponding parameters in the notation of the SFORB model (Single +First-Order Reversible Binding) are additionally calculated. + +<>= +kinfits.C <- kinfit(FOCUS_2006_C, kinmodels = c("SFO", "FOMC", "DFOP", "HS")) +@ + +The results of the fitting procedure are returned by the function, and can +then be inspected by the function \Robject{kinresults}. + +<>= +kinresults(kinfits.C) +@ + +The higher level functions \Robject{kinplot} and \Robject{kinreport} +work on lists called \Robject{kinobject}. They contain the fitted models, +optionally the data used for fitting the models, and the name of the parent +compound as well as the test system type used for generating the data, +as well as some more optional entries. The construction of such an +object is shown below. + +<>= +kinobject.C <- kinobject <- list( + parent = "Compound XY", + type = "Degradation in the environment", + system = "System 1", + source = "Synthetic example data from FOCUS kinetics", + data = FOCUS_2006_C, + fits = kinfits.C, + results = kinresults(kinfits.C)) +@ + +The plotting and reporting functions then work on this object. The example +below outputs the report to the console, because no \Robject{file} argument +is specified. If a filename is specified, the report will be written to +a text file. + +<>= +kinreport(kinobject.C) +@ + +Plotting is done on an on-screen device. Graphics files in vector based +formats can be obtained using the \RR{} devices \Robject{pdf}, \Robject{eps}, +or, subject to platform restrictions, \Robject{windows.metafile}. + +\setkeys{Gin}{width=0.6\textwidth} +\begin{figure}[t] +\begin{center} +<>= +kinfits.C <- kinfit(FOCUS_2006_C, kinmodels = c("SFO", "FOMC", "DFOP", "HS")) +kinplot(kinobject.C) +title("FOCUS dataset C") +@ +\caption{Fits of standard models to FOCUS dataset C.} +\label{fig:FOCUS_2006_C} +\end{center} +\end{figure} + +A residual plot can be obtained with the function \Robject{kinresplot} as +shown in Figure \ref{fig:FOCUS_2006_C_res}. + +\setkeys{Gin}{width=0.6\textwidth} +\begin{figure}[t] +\begin{center} +<>= +kinresplot(kinobject.C, "SFO") +@ +\caption{Residual plot for fitting the SFO model to FOCUS dataset C.} +\label{fig:FOCUS_2006_C_res} +\end{center} +\end{figure} + +\section{Validation} +\label{vali} + +In the following comparisons, the results for fitting the four recommended +kinetic models to FOCUS datasets A to F with \Rpackage{kinfit} were +obtained. + +<>= +datasets <- LETTERS[1:4] +data(list=paste("FOCUS_2006_", datasets, sep=""), package = "kinfit") +kinobjects <- list() +for (dataset in datasets) +{ + kinobjects[[dataset]] <- list() + kinobjects[[dataset]]$data <- get(paste("FOCUS_2006_", dataset, sep="")) + kinobjects[[dataset]]$fits <- + kinfit(kinobjects[[dataset]]$data, + kinmodels = c("SFO", "FOMC", "DFOP", "HS")) + kinobjects[[dataset]]$results <- + kinresults(kinobjects[[dataset]]$fits) +} + +data(FOCUS_2006_F, package = "kinfit") +# Set the initial concentration in the sediment to zero +FOCUS_2006_F[1, "parent.sediment"] <- 0 +# Calculate total system values for the parent compound +FOCUS_2006_F = transform(FOCUS_2006_F, + parent.system = parent.water + parent.sediment) + +subsets <- c("system", "water") +for (subset in subsets) +{ + index <- paste("F", subset, sep=" ") + kinobjects[[index]] <- list() + kinobjects[[index]]$data <- data.frame( + t = FOCUS_2006_F$t, + parent = FOCUS_2006_F[[paste("parent", subset, sep=".")]]) + kinobjects[[index]]$fits <- + kinfit(kinobjects[[index]]$data, + kinmodels = c("SFO", "FOMC", "DFOP", "HS")) + kinobjects[[index]]$results <- + kinresults(kinobjects[[index]]$fits) +} +@ + +\subsection{Single First Order Model} + +In Tables \ref{tab:vali.SFO.A} to \ref{tab:vali.SFO.F_water}, +the results from fitting the SFO model to FOCUS example datasets with +various software packages as given in the report by the \cite{FOCUS2006} are +compared with the results obtained with \Rpackage{kinfit}. + +<>= +data("FOCUS_2006_SFO_ref_A_to_F", package = "kinfit") +kinmodel = "SFO" +refs <- list() +for (kinobjectname in names(kinobjects)) +{ + ref <- subset(FOCUS_2006_SFO_ref_A_to_F, dataset == kinobjectname) + ref <- ref[-6] + texfile <- paste("FOCUS_2006_", kinmodel, "_", + gsub(" ", "_", kinobjectname), "_ref.tex", sep="") + write.table(format(ref, nsmall=2), + file = texfile, + sep=" & ", quote=FALSE, + row.names=FALSE, col.names=FALSE, eol = " \\\\ \n") + refs[[kinobjectname]] <- ref +} +@ + +\begin{table} +\caption{Results of fitting the SFO model to the example dataset A +\citep{FOCUS2006}, as given in the report, in comparison to the results +obtained by \Rpackage{kinfit}. \label{tab:vali.SFO.A}} +\begin{center} +\vspace{0.5cm} +\begin{tabular}{lcccc} +\toprule +Package & $M_0$ & $k$ & DT$_{50}$ & DT$_{90}$ \\ +\midrule +\input{FOCUS_2006_SFO_A_ref} +\midrule +Median & +\Sexpr{format(median(refs$A$M0), nsmall=2)} & +\Sexpr{format(median(refs$A$k), nsmall=4)} & +\Sexpr{format(median(refs$A$DT50), nsmall=2)} & +\Sexpr{format(median(refs$A$DT90), nsmall=2)} \\ +\midrule +\Rpackage{kinfit} & +\Sexpr{round(kinobjects$A$results$parms$SFO$parent.0, 2)} & +\Sexpr{round(kinobjects$A$results$parms$SFO$k, 4)} & +\Sexpr{round(kinobjects$A$results$results[["SFO", "DT50"]], 2)} & +\Sexpr{round(kinobjects$A$results$results[["SFO", "DT90"]], 2)} \\ +\bottomrule +\end{tabular} +\end{center} +\end{table} + +\begin{table} +\caption{Results of fitting the SFO model to the example dataset B +\citep{FOCUS2006}, as given in the report, in comparison to the results +obtained by \Rpackage{kinfit}. \label{tab:vali.SFO.B}} +\begin{center} +\vspace{0.5cm} +\begin{tabular}{lcccc} +\toprule +Package & $M_0$ & $k$ & DT$_{50}$ & DT$_{90}$ \\ +\midrule +\input{FOCUS_2006_SFO_B_ref} +\midrule +Median & +\Sexpr{format(median(refs$B$M0), nsmall=2)} & +\Sexpr{format(median(refs$B$k), nsmall=4)} & +\Sexpr{format(median(refs$B$DT50), nsmall=2)} & +\Sexpr{format(median(refs$B$DT90), nsmall=2)} \\ +\midrule +\Rpackage{kinfit} & +\Sexpr{round(kinobjects$B$results$parms$SFO$parent.0, 2)} & +\Sexpr{round(kinobjects$B$results$parms$SFO$k, 4)} & +\Sexpr{round(kinobjects$B$results$results[["SFO", "DT50"]], 2)} & +\Sexpr{round(kinobjects$B$results$results[["SFO", "DT90"]], 2)} \\ +\bottomrule +\end{tabular} +\end{center} +\end{table} + +\begin{table} +\caption{Results of fitting the SFO model to the example dataset C +\citep{FOCUS2006}, as given in the report, in comparison to the results +obtained by \Rpackage{kinfit}. \label{tab:vali.SFO.C}} +\begin{center} +\vspace{0.5cm} +\begin{tabular}{lcccc} +\toprule +Package & $M_0$ & $k$ & DT$_{50}$ & DT$_{90}$ \\ +\midrule +\input{FOCUS_2006_SFO_C_ref} +\midrule +Median & +\Sexpr{format(median(refs$C$M0), nsmall=2)} & +\Sexpr{format(median(refs$C$k), nsmall=4)} & +\Sexpr{format(median(refs$C$DT50), nsmall=2)} & +\Sexpr{format(median(refs$C$DT90), nsmall=2)} \\ +\midrule +\Rpackage{kinfit} & +\Sexpr{round(kinobjects$C$results$parms$SFO$parent.0, 2)} & +\Sexpr{round(kinobjects$C$results$parms$SFO$k, 4)} & +\Sexpr{round(kinobjects$C$results$results[["SFO", "DT50"]], 2)} & +\Sexpr{round(kinobjects$C$results$results[["SFO", "DT90"]], 2)} \\ +\bottomrule +\end{tabular} +\end{center} +\end{table} + +\begin{table} +\caption{Results of fitting the SFO model to the example dataset D +\citep{FOCUS2006}, as given in the report, in comparison to the results +obtained by \Rpackage{kinfit}. \label{tab:vali.SFO.D}} +\begin{center} +\vspace{0.5cm} +\begin{tabular}{lcccc} +\toprule +Package & $M_0$ & $k$ & DT$_{50}$ & DT$_{90}$ \\ +\midrule +\input{FOCUS_2006_SFO_D_ref} +\midrule +Median & +\Sexpr{format(median(refs$D$M0), nsmall=2)} & +\Sexpr{format(median(refs$D$k), nsmall=4)} & +\Sexpr{format(median(refs$D$DT50), nsmall=2)} & +\Sexpr{format(median(refs$D$DT90), nsmall=2)} \\ +\midrule +\Rpackage{kinfit} & +\Sexpr{round(kinobjects$D$results$parms$SFO$parent.0, 2)} & +\Sexpr{round(kinobjects$D$results$parms$SFO$k, 4)} & +\Sexpr{round(kinobjects$D$results$results[["SFO", "DT50"]], 2)} & +\Sexpr{round(kinobjects$D$results$results[["SFO", "DT90"]], 2)} \\ +\bottomrule +\end{tabular} +\end{center} +\end{table} + +\begin{table} +\caption{Results of fitting the SFO model to the total system data from +example dataset F \citep{FOCUS2006}, as given in the report, in comparison to +the results obtained by \Rpackage{kinfit}. \label{tab:vali.SFO.F_system}} +\begin{center} +\vspace{0.5cm} +\begin{tabular}{lcccc} +\toprule +Package & $M_0$ & $k$ & DT$_{50}$ & DT$_{90}$ \\ +\midrule +\input{FOCUS_2006_SFO_F_system_ref} +\midrule +Median & +\Sexpr{format(median(refs[["F system"]]$M0), nsmall=2)} & +\Sexpr{format(median(refs[["F system"]]$k), nsmall=4)} & +\Sexpr{format(median(refs[["F system"]]$DT50), nsmall=2)} & +\Sexpr{format(median(refs[["F system"]]$DT90), nsmall=2)} \\ +\midrule +\Rpackage{kinfit} & +\Sexpr{round(kinobjects[["F system"]]$results$parms$SFO$parent.0, 2)} & +\Sexpr{round(kinobjects[["F system"]]$results$parms$SFO$k, 4)} & +\Sexpr{round(kinobjects[["F system"]]$results$results[["SFO", "DT50"]], 2)} & +\Sexpr{round(kinobjects[["F system"]]$results$results[["SFO", "DT90"]], 2)} \\ +\bottomrule +\end{tabular} +\end{center} +\end{table} + +\begin{table} +\caption{Results of fitting the SFO model to the water phase data from +example dataset F \citep{FOCUS2006}, as given in the report, in comparison to +the results obtained by \Rpackage{kinfit}. \label{tab:vali.SFO.F_water}} +\begin{center} +\vspace{0.5cm} +\begin{tabular}{lcccc} +\toprule +Package & $M_0$ & $k$ & DT$_{50}$ & DT$_{90}$ \\ +\midrule +\input{FOCUS_2006_SFO_F_water_ref} +\midrule +Median & +\Sexpr{format(median(refs[["F water"]]$M0), nsmall=2)} & +\Sexpr{format(median(refs[["F water"]]$k), nsmall=4)} & +\Sexpr{format(median(refs[["F water"]]$DT50), nsmall=2)} & +\Sexpr{format(median(refs[["F water"]]$DT90), nsmall=2)} \\ +\midrule +\Rpackage{kinfit} & +\Sexpr{round(kinobjects[["F water"]]$results$parms$SFO$parent.0, 2)} & +\Sexpr{round(kinobjects[["F water"]]$results$parms$SFO$k, 4)} & +\Sexpr{round(kinobjects[["F water"]]$results$results[["SFO", "DT50"]], 2)} & +\Sexpr{format(kinobjects[["F water"]]$results$results[["SFO", "DT90"]], digits=4, nsmall=2)} \\ +\bottomrule +\end{tabular} +\end{center} +\end{table} + +The comparisons show that all packages evaluated in the FOCUS report give +very similar results for the SFO model. The results obtained with +\Rpackage{kinfit} are very close to the median of the results reported for +the other packages. + +\subsection{First Order Multi Compartment Model} + +<>= +data("FOCUS_2006_FOMC_ref_A_to_F", package = "kinfit") +kinmodel = "FOMC" +refs <- list() +for (kinobjectname in names(kinobjects)[c(1:3, 5:6)]) +{ + ref <- subset(FOCUS_2006_FOMC_ref_A_to_F, dataset == kinobjectname) + ref$package <- gsub("#", "$^a$", ref$package) + ref <- ref[-7] + texfile <- paste("FOCUS_2006_", kinmodel, "_", + gsub(" ", "_", kinobjectname), "_ref.tex", sep="") + write.table(format(ref, nsmall=2), + file = texfile, + sep=" & ", quote=FALSE, + row.names=FALSE, col.names=FALSE, eol = " \\\\ \n") + refs[[kinobjectname]] <- ref +} +@ + +\begin{table} +\caption{Results of fitting the FOMC model to the example dataset A +\citep{FOCUS2006}, as given in the report, in comparison to the results +obtained by \Rpackage{kinfit}. \label{tab:vali.FOMC.A}} +\begin{center} +\vspace{0.5cm} +\begin{tabular}{lccccc} +\toprule +Package & $M_0$ & $\alpha$ & $\beta$ & DT$_{50}$ & DT$_{90}$ \\ +\midrule +\input{FOCUS_2006_FOMC_A_ref} +\midrule +Median & +\Sexpr{format(median(refs$A$M0), nsmall=2)} & +\Sexpr{format(median(refs$A$alpha), scientific=TRUE)} & +\Sexpr{format(median(as.numeric(refs$A$beta)), nsmall=0)} & +\Sexpr{format(median(refs$A$DT50), nsmall=2)} & +\Sexpr{format(median(refs$A$DT90), nsmall=2)} \\ +\midrule +\Rpackage{kinfit} & +no fit \\ +\bottomrule +\end{tabular} +\end{center} +\end{table} + +\begin{table} +\caption{Results of fitting the FOMC model to the example dataset B +\citep{FOCUS2006}, as given in the report, in comparison to the results +obtained by \Rpackage{kinfit}. \label{tab:vali.FOMC.B}} +\begin{center} +\vspace{0.5cm} +\begin{tabular}{lccccc} +\toprule +Package & $M_0$ & $\alpha$ & $\beta$ & DT$_{50}$ & DT$_{90}$ \\ +\midrule +\input{FOCUS_2006_FOMC_B_ref} +\midrule +Median & +\Sexpr{format(median(refs$B$M0), digits=4, nsmall=2)} & +\Sexpr{format(median(refs$B$alpha), scientific=TRUE)} & +\Sexpr{format(median(as.numeric(refs$B$beta)), nsmall=0)} & +\Sexpr{format(median(refs$B$DT50), nsmall=2)} & +\Sexpr{format(median(refs$B$DT90), digits=4, nsmall=2)} \\ +\midrule +\Rpackage{kinfit} & +\Sexpr{round(kinobjects$B$results$parms$FOMC$parent.0, 2)} & +\Sexpr{round(kinobjects$B$results$parms$FOMC$alpha, 4)} & +\Sexpr{round(kinobjects$B$results$parms$FOMC$beta, 4)} & +\Sexpr{round(kinobjects$B$results$results[["FOMC", "DT50"]], 2)} & +\Sexpr{format(kinobjects$B$results$results[["FOMC", "DT90"]], digits=4, nsmall=2)} \\ +\bottomrule +\end{tabular} +\end{center} +\end{table} + + +\begin{table} +\caption{Results of fitting the FOMC model to the example dataset C +\citep{FOCUS2006}, as given in the report, in comparison to the results +obtained by \Rpackage{kinfit}. \label{tab:vali.FOMC.C}} +\begin{center} +\vspace{0.5cm} +\begin{tabular}{lccccc} +\toprule +Package & $M_0$ & $\alpha$ & $\beta$ & DT$_{50}$ & DT$_{90}$ \\ +\midrule +\input{FOCUS_2006_FOMC_C_ref} +\midrule +Median & +\Sexpr{format(median(refs$C$M0), nsmall=2)} & +\Sexpr{format(median(refs$C$alpha), nsmall=2)} & +\Sexpr{format(median(as.numeric(refs$C$beta)), nsmall=2)} & +\Sexpr{format(median(refs$C$DT50), nsmall=2)} & +\Sexpr{format(median(refs$C$DT90), nsmall=2)} \\ +\midrule +\Rpackage{kinfit} & +\Sexpr{round(kinobjects$C$results$parms$FOMC$parent.0, 2)} & +\Sexpr{round(kinobjects$C$results$parms$FOMC$alpha, 4)} & +\Sexpr{round(kinobjects$C$results$parms$FOMC$beta, 4)} & +\Sexpr{round(kinobjects$C$results$results[["FOMC", "DT50"]], 2)} & +\Sexpr{format(kinobjects$C$results$results[["FOMC", "DT90"]], digits=4, nsmall=2)} \\ +\bottomrule +\end{tabular} +\end{center} +\end{table} + +\begin{table} +\caption{Results of fitting the FOMC model to the total system data from +example dataset F \citep{FOCUS2006}, as given in the report, in comparison to +the results obtained by \Rpackage{kinfit}. \label{tab:vali.FOMC.F_system}} +\begin{center} +\vspace{0.5cm} +\begin{tabular}{lccccc} +\toprule +Package & $M_0$ & $\alpha$ & $\beta$ & DT$_{50}$ & DT$_{90}$ \\ +\midrule +\input{FOCUS_2006_FOMC_F_system_ref} +\midrule +Median & +\Sexpr{format(median(refs[["F system"]]$M0), nsmall=2)} & +\Sexpr{format(median(refs[["F system"]]$alpha), nsmall=4)} & +\Sexpr{format(median(refs[["F system"]]$beta), nsmall=4)} & +\Sexpr{format(median(refs[["F system"]]$DT50), nsmall=2)} & +\Sexpr{format(median(refs[["F system"]]$DT90), nsmall=2)} \\ +\midrule +\Rpackage{kinfit} & +no fit \\ +\bottomrule +\end{tabular} +\end{center} +\end{table} + +\begin{table} +\caption{Results of fitting the FOMC model to the water phase data from +example dataset F \citep{FOCUS2006}, as given in the report, in comparison to +the results obtained by \Rpackage{kinfit}. \label{tab:vali.FOMC.F_water}} +\begin{center} +\vspace{0.5cm} +\begin{tabular}{lccccc} +\toprule +Package & $M_0$ & $\alpha$ & $\beta$ & DT$_{50}$ & DT$_{90}$ \\ +\midrule +\input{FOCUS_2006_FOMC_F_water_ref} +\midrule +Median & +\Sexpr{format(median(refs[["F water"]]$M0), nsmall=2)} & +\Sexpr{format(median(refs[["F water"]]$alpha), nsmall=4)} & +\Sexpr{format(median(refs[["F water"]]$beta), nsmall=4)} & +\Sexpr{format(median(refs[["F water"]]$DT50), nsmall=2)} & +\Sexpr{format(median(refs[["F water"]]$DT90), nsmall=2)} \\ +\midrule +\Rpackage{kinfit} & +no fit \\ +\bottomrule +\end{tabular} +\end{center} +\end{table} + +The comparison of the results obtained for the FOMC model show much more +variability between software packages. For dataset A, results for the +\Robject{alpha} and \Robject{beta} parameters differ over several +orders of magnitude between the different packages. The method used +by the \Robject{kinfit} routine does not converge for this dataset. +The same applies to the total system and water phase only data for +example dataset F and the FOMC model. + +For datasets B and C, the \Rpackage{kinfit} function produces results +which are very close to the median of the results obtained by the +other packages. + +\subsection{Dual First Order in Parallel Model} + +<>= +data("FOCUS_2006_DFOP_ref_A_to_B", package = "kinfit") +kinmodel = "DFOP" +refs <- list() +for (kinobjectname in names(kinobjects)[1:2]) +{ + ref <- subset(FOCUS_2006_DFOP_ref_A_to_B, dataset == kinobjectname) + ref <- ref[-8] + texfile <- paste("FOCUS_2006_", kinmodel, "_", + gsub(" ", "_", kinobjectname), "_ref.tex", sep="") + write.table(format(ref, nsmall=2), + file = texfile, + sep=" & ", quote=FALSE, + row.names=FALSE, col.names=FALSE, eol = " \\\\ \n") + refs[[kinobjectname]] <- ref +} +@ + +\begin{table} +\caption{Results of fitting the DFOP model to the example dataset A +\citep{FOCUS2006}, as given in the report, in comparison to the results +obtained by \Rpackage{kinfit}. \label{tab:vali.DFOP.A}} +\begin{center} +\vspace{0.5cm} +\begin{tabular}{lcccccc} +\toprule +Package & $M_0$ & $f$ & $k_1$ & $k_2$ & DT$_{50}$ & DT$_{90}$ \\ +\midrule +\input{FOCUS_2006_DFOP_A_ref} +\midrule +Median & +\Sexpr{format(median(refs$A$M0), nsmall=2)} & +\Sexpr{format(median(refs$A$f), nsmall=2)} & +\Sexpr{format(median(refs$A$k1), nsmall=4)} & +\Sexpr{format(median(refs$A$k2), nsmall=4)} & +\Sexpr{format(median(refs$A$DT50), nsmall=2)} & +\Sexpr{format(median(refs$A$DT90), nsmall=2)} \\ +\midrule +\Rpackage{kinfit} & +no fit \\ +\bottomrule +\end{tabular} +\end{center} +\end{table} + +\begin{table} +\caption{Results of fitting the DFOP model to the example dataset B +\citep{FOCUS2006}, as given in the report, in comparison to the results +obtained by \Rpackage{kinfit}. \label{tab:vali.DFOP.B}} +\begin{center} +\vspace{0.5cm} +\begin{tabular}{lcccccc} +\toprule +Package & $M_0$ & $f$ & $k_1$ & $k_2$ & DT$_{50}$ & DT$_{90}$ \\ +\midrule +\input{FOCUS_2006_DFOP_B_ref} +\midrule +Median & +\Sexpr{format(median(refs$B$M0), nsmall=2)} & +\Sexpr{format(median(refs$B$f), nsmall=2)} & +\Sexpr{format(median(refs$B$k1), nsmall=4)} & +\Sexpr{format(median(refs$B$k2), nsmall=4)} & +\Sexpr{format(median(refs$B$DT50), nsmall=2)} & +\Sexpr{format(median(refs$B$DT90), digits=4, nsmall=2)} \\ +\midrule +\Rpackage{kinfit} & +\Sexpr{round(kinobjects$B$results$parms$DFOP$parent.0, 2)} & +\Sexpr{round(kinobjects$B$results$parms$DFOP$g, 2)} & +\Sexpr{round(kinobjects$B$results$parms$DFOP$k1, 4)} & +\Sexpr{round(kinobjects$B$results$parms$DFOP$k2, 4)} & +\Sexpr{round(kinobjects$B$results$results[["DFOP", "DT50"]], 2)} & +\Sexpr{format(kinobjects$B$results$results[["DFOP", "DT90"]], digits=4, nsmall=2)} \\ +\bottomrule +\end{tabular} +\end{center} +\end{table} + +Regarding fitting the DFOP model to FOCUS example dataset A, it is already +indicated in the report that it is not a good example dataset for fitting +this particular model, as the two kinetic constants postulated by the DFOP +model are hardly distinguishable. As a consequence, the software packages +strongly disagree especially on the model parameter $f$ specifying the +distribution between the kinetic domains that are characterised by the +two kinetic constants. Again, the \Rpackage{kinfit} routine does not +show conversion for this model and this dataset (Table \ref{tab:vali.DFOP.A}). + +Fitting the DFOP model with \Rpackage{kinfit} to dataset B yields results +that are very close to the median of the results obtained by other packages, +as illustrated in Table \ref{tab:vali.DFOP.B}. + +\subsection{Hockey Stick Model} + +Analysis of dataset A shows basically two different parameter sets +generated by the 8 packages reported in the FOCUS report \citep{FOCUS2006}. +The \Rpackage{kinfit} package does not show conversion with the standard +paramater defaults, but can reproduce the two parameter sets when given +the respective paramter values as starting values, as shown in the last +two lines in Table \ref{tab:vali.HS.A}. + +<>= +data("FOCUS_2006_HS_ref_A_to_F", package = "kinfit") +kinmodel = "HS" +refs <- list() +for (kinobjectname in names(kinobjects)[c(1:3, 5:6)]) +{ + ref <- subset(FOCUS_2006_HS_ref_A_to_F, dataset == kinobjectname) + ref$package <- gsub("\\*", "$^a$", ref$package) + ref <- ref[-8] + texfile <- paste("FOCUS_2006_", kinmodel, "_", + gsub(" ", "_", kinobjectname), "_ref.tex", sep="") + write.table(format(ref, nsmall=2), + file = texfile, + sep=" & ", quote=FALSE, + row.names=FALSE, col.names=FALSE, eol = " \\\\ \n") + refs[[kinobjectname]] <- ref +} +@ + +<>= +kinobjects$A$fits.2 <- kinfit(kinobjects$A$data, + kinmodels = c("HS"), + start.HS = list(parent.0 = 100, tb = 5, k1 = 0.017, k2 = 0.05)) +kinobjects$A$results.2 <- kinresults(kinobjects$A$fits.2) + +kinobjects$A$fits.3 <- kinfit(kinobjects$A$data, + kinmodels = c("HS"), + start.HS = list(parent.0 = 100, tb = 11, k1 = 0.017, k2 = 0.05)) +kinobjects$A$results.3 <- kinresults(kinobjects$A$fits.3) +@ + +\begin{table} +\caption{Results of fitting the HS model to the example dataset A +\citep{FOCUS2006}, as given in the report, in comparison to the results +obtained by \Rpackage{kinfit}. \label{tab:vali.HS.A}} +\begin{center} +\vspace{0.5cm} +\begin{tabular}{lcccccc} +\toprule +Package & $M_0$ & $t_b$ & $k_1$ & $k_2$ & DT$_{50}$ & DT$_{90}$ \\ +\midrule +\input{FOCUS_2006_HS_A_ref} +\midrule +Median & +\Sexpr{format(median(refs$A$M0), nsmall=2)} & +\Sexpr{format(median(refs$A$tb), nsmall=2)} & +\Sexpr{format(median(refs$A$k1), nsmall=4)} & +\Sexpr{format(median(refs$A$k2), nsmall=4)} & +\Sexpr{format(median(refs$A$DT50), nsmall=2)} & +\Sexpr{format(median(refs$A$DT90), digits=4, nsmall=2)} \\ +\midrule +\Rpackage{kinfit} & no fit \\ +\Rpackage{kinfit} & +\Sexpr{round(kinobjects$A$results.2$parms$HS$parent.0, 2)} & +\Sexpr{round(kinobjects$A$results.2$parms$HS$tb, 2)} & +\Sexpr{round(kinobjects$A$results.2$parms$HS$k1, 4)} & +\Sexpr{round(kinobjects$A$results.2$parms$HS$k2, 4)} & +\Sexpr{round(kinobjects$A$results.2$results[["HS", "DT50"]], 2)} & +\Sexpr{format(kinobjects$A$results.2$results[["HS", "DT90"]], digits=4, nsmall=2)} \\ +\Rpackage{kinfit} & +\Sexpr{round(kinobjects$A$results.3$parms$HS$parent.0, 2)} & +\Sexpr{round(kinobjects$A$results.3$parms$HS$tb, 2)} & +\Sexpr{round(kinobjects$A$results.3$parms$HS$k1, 4)} & +\Sexpr{round(kinobjects$A$results.3$parms$HS$k2, 4)} & +\Sexpr{round(kinobjects$A$results.3$results[["HS", "DT50"]], 2)} & +\Sexpr{format(kinobjects$A$results.3$results[["HS", "DT90"]], digits=4, nsmall=2)} \\ +\bottomrule +\end{tabular} +\end{center} +\end{table} + +\begin{table} +\caption{Results of fitting the HS model to the example dataset B +\citep{FOCUS2006}, as given in the report, in comparison to the results +obtained by \Rpackage{kinfit}. \label{tab:vali.HS.B}} +\begin{center} +\vspace{0.5cm} +\begin{tabular}{lcccccc} +\toprule +Package & $M_0$ & $t_b$ & $k_1$ & $k_2$ & DT$_{50}$ & DT$_{90}$ \\ +\midrule +\input{FOCUS_2006_HS_B_ref} +\midrule +Median & +\Sexpr{format(median(refs$B$M0), nsmall=2)} & +\Sexpr{format(median(refs$B$tb), nsmall=2)} & +\Sexpr{format(median(refs$B$k1), nsmall=4)} & +\Sexpr{format(median(refs$B$k2), nsmall=4)} & +\Sexpr{format(median(refs$B$DT50), nsmall=2)} & +\Sexpr{format(median(refs$B$DT90), digits=4, nsmall=2)} \\ +\midrule +\Rpackage{kinfit} & no fit \\ +\bottomrule +\end{tabular} +\end{center} +\end{table} + +\begin{table} +\caption{Results of fitting the HS model to the example dataset C +\citep{FOCUS2006}, as given in the report, in comparison to the results +obtained by \Rpackage{kinfit}. \label{tab:vali.HS.C}} +\begin{center} +\vspace{0.5cm} +\begin{tabular}{lcccccc} +\toprule +Package & $M_0$ & $t_b$ & $k_1$ & $k_2$ & DT$_{50}$ & DT$_{90}$ \\ +\midrule +\input{FOCUS_2006_HS_C_ref} +\midrule +Median & +\Sexpr{format(median(refs$C$M0), nsmall=2)} & +\Sexpr{format(median(refs$C$tb), nsmall=2)} & +\Sexpr{format(median(refs$C$k1), nsmall=4)} & +\Sexpr{format(median(as.numeric(refs$C$k2)), nsmall=4)} & +\Sexpr{format(median(refs$C$DT50), nsmall=2)} & +\Sexpr{format(median(refs$C$DT90), digits=4, nsmall=2)} \\ +\midrule +\Rpackage{kinfit} & +\Sexpr{round(kinobjects$C$results$parms$HS$parent.0, 2)} & +\Sexpr{round(kinobjects$C$results$parms$HS$tb, 2)} & +\Sexpr{round(kinobjects$C$results$parms$HS$k1, 4)} & +\Sexpr{round(kinobjects$C$results$parms$HS$k2, 4)} & +\Sexpr{round(kinobjects$C$results$results[["HS", "DT50"]], 2)} & +\Sexpr{format(kinobjects$C$results$results[["HS", "DT90"]], digits=4, nsmall=2)} \\ +\bottomrule +\end{tabular} +\end{center} +\end{table} + +The HS fit did not converge for dataset B with \Rpackage{kinfit}. Again, this +should be viewed in the light of the vastly differing results produced by +the other software packages as listed in Table \ref{tab:vali.HS.B}. + +The results from fitting the HS model to dataset C with \Rpackage{kinfit} +agree nicely with the median of the results obtained with the other packages, +as shown in Table \ref{tab:vali.HS.C}. + +\subsection{$\chi^2$ statistics} + +As no values for the minimum error rate that has to be assumed for the +model to agree with the data ($\chi^2$ statistics) are reported for +the FOCUS datasets A to F, the respective values calculated by +\Rpackage{kinfit} are compared to the $\chi^2$ values calculated by +the KinGUI package \citep{schaefer2007} as shown in Table +\ref{tab:vali.chi2}. + +For this, the possibility to write KinGUI input files using the function +\Robject{kinwrite.KinGUI} from \Rpackage{kinfit} was used. + +<>= +chi2.SFO.kinfit <- chi2.FOMC.kinfit <- array(dim = length(kinobjects), + dimnames = list(names(kinobjects))) +chi2.DFOP.kinfit <- chi2.HS.kinfit <- array(dim = length(kinobjects), + dimnames = list(names(kinobjects))) +for (kinobjectname in names(kinobjects)) +{ + outname <- paste("KinGUI/", gsub(" ", "_", kinobjectname), "_KinGUI.txt", + sep="") + kinwrite.KinGUI(kinobjects[[kinobjectname]], outname) + chi2.SFO.kinfit[[kinobjectname]] <- + kinobjects[[kinobjectname]]$results$stats[["SFO", "err.min"]] + chi2.FOMC.kinfit[[kinobjectname]] <- + ifelse(class(kinobjects[[kinobjectname]]$fits$FOMC) == "try-error", + NA, kinobjects[[kinobjectname]]$results$stats[["FOMC", "err.min"]]) + chi2.DFOP.kinfit[[kinobjectname]] <- + ifelse(class(kinobjects[[kinobjectname]]$fits$DFOP) == "try-error", + NA, kinobjects[[kinobjectname]]$results$stats[["DFOP", "err.min"]]) + chi2.HS.kinfit[[kinobjectname]] <- + ifelse(class(kinobjects[[kinobjectname]]$fits$HS) == "try-error", + NA, kinobjects[[kinobjectname]]$results$stats[["HS", "err.min"]]) +} + +chi2.SFO.KinGUI <- c(8.3852, 4.4562, 15.8456, 6.4539, 12.5386, 10.8069) +chi2.FOMC.KinGUI <- c(9.3116, 4.6641, 6.6574, 6.8080, 13.4533, 11.6682) +chi2.DFOP.KinGUI <- c(9.6600, 4.9562, 2.6613, 7.2751, 14.1524, 12.1821) +chi2.HS.KinGUI <- c(4.1106, 4.4535, 4.6963, 5.8196, 3.2178, 1.6558) +names(chi2.SFO.KinGUI) <- names(chi2.FOMC.KinGUI) <- names(kinobjects) +names(chi2.DFOP.KinGUI) <- names(chi2.HS.KinGUI) <- names(kinobjects) + +chi2 <- data.frame( + SFO.KinGUI = chi2.SFO.KinGUI, + SFO.kinfit = round(100 * chi2.SFO.kinfit, 4), + FOMC.KinGUI = chi2.FOMC.KinGUI, + FOMC.kinfit = round(100 * chi2.FOMC.kinfit, 4), + DFOP.KinGUI = chi2.DFOP.KinGUI, + DFOP.kinfit = round(100 * chi2.DFOP.kinfit, 4), + HS.KinGUI = chi2.HS.KinGUI, + HS.kinfit = round(100 * chi2.HS.kinfit, 4) +) +write.table(chi2, + file = "chi2_comparison.tex", + sep=" & ", quote=FALSE, na="", + row.names=TRUE, col.names=FALSE, eol = " \\\\ \n") +@ + +\begin{table} +\caption{Comparison of $\chi^2$ error levels in percent calculated for model +fits by the KinGUI and \Rpackage{kinfit} packages. \label{tab:vali.chi2}} +\vspace{0.5cm} +\begin{tabular}{lcccccccc} +\toprule + & \multicolumn{2}{c}{SFO} & +\multicolumn{2}{c}{FOMC} & +\multicolumn{2}{c}{DFOP} & +\multicolumn{2}{c}{HS} \\ +Dataset & KinGUI & \Rpackage{kinfit} & KinGUI & \Rpackage{kinfit} & +KinGUI & \Rpackage{kinfit} & KinGUI & \Rpackage{kinfit} \\ +\midrule +\input{chi2_comparison} +\bottomrule +\end{tabular} +\end{table} + +The comparison shows that whenever a minimum error level $\chi^2$ was +calculated using the \Rpackage{kinfit} package, it was very close to +the value generated by KinGUI. + +\section{Conclusion} + +The \Rpackage{kinfit} package for \RR{} gives access to the possibility to +fit the kinetic models recommended by the FOCUS group \citep{FOCUS2006} +from within \RR. Comparison with the results obtained with other +software packages shows that \Rpackage{kinfit} produces kinetic endpoints that +are within the variability and even very close to the median of results obtained +with other packages, except for some cases where \Rpackage{kinfit} does not +produce results and the results obtained with other software packages are +strongly divergent. + +\bibliographystyle{plainnat} +\bibliography{references} + +\end{document} +% vim: set foldmethod=marker: -- cgit v1.2.1