% $Id: kinfit.Rnw 125 2011-11-10 07:19:59Z jranke $
%%\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}
<<setup, echo = FALSE, results = hide>>=
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}
\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 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

<<FOCUS_2006_C_data, echo=TRUE, eval=TRUE>>=
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

<<data_format, echo=TRUE>>=
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.

<<FOCUS_2006_C_fits, echo=TRUE>>=
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}.

<<FOCUS_2006_C_results, echo = TRUE>>=
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.

<<FOCUS_2006_C_object, echo=TRUE>>=
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.

<<FOCUS_2006_C_report, echo = TRUE>>=
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}
<<FOCUS_2006_C_figure, echo = FALSE, fig = TRUE, width = 7, height = 5, cex = TRUE>>=
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}
<<FOCUS_2006_C_res, echo = FALSE, fig = TRUE, width = 7, height = 5, cex = TRUE>>=
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.

<<kinfits_A_to_F, echo = FALSE>>=
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}.

<<SFO, echo = FALSE>>=
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}

<<FOMC, echo = FALSE>>=
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{Double First Order in Parallel Model}

<<DFOP, echo = FALSE>>=
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 convergence 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}.

<<HS, echo = FALSE>>=
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
}
@

<<HS2, echo = FALSE>>=
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.

<<KinGUI_write, echo=FALSE>>=
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.

\subsection{$R^2$ value}

In a similar manner, the coefficient of determination calculated by kinfit
according to the formula

\begin{equation}
R^2 = 1 - RSS / TSS
\end{equation}

where RSS is the sum of squares of the residuals and TSS is the sum of squares
of the deviations from the mean value is compared to the model efficiency EF
calculated by KinGUI based on the same formula (FOCUS 2006, p. 99) as shown
in Table \ref{tab:vali.R2}.  This exercise was done for FOCUS datasets A to D
only.

<<KinGUI_write, echo=FALSE>>=
R2.SFO.kinfit <- R2.FOMC.kinfit <- array(dim = 4, 
  dimnames = list(names(kinobjects[1:4])))
R2.DFOP.kinfit <- R2.HS.kinfit <- array(dim = 4,
  dimnames = list(names(kinobjects[1:4])))
for (kinobjectname in names(kinobjects[1:4]))
{
  R2.SFO.kinfit[[kinobjectname]] <- 
    kinobjects[[kinobjectname]]$results$stats[["SFO", "R2"]]
  R2.FOMC.kinfit[[kinobjectname]] <- 
    ifelse(class(kinobjects[[kinobjectname]]$fits$FOMC) == "try-error",
      NA, kinobjects[[kinobjectname]]$results$stats[["FOMC", "R2"]])
  R2.DFOP.kinfit[[kinobjectname]] <- 
    ifelse(class(kinobjects[[kinobjectname]]$fits$DFOP) == "try-error",
      NA, kinobjects[[kinobjectname]]$results$stats[["DFOP", "R2"]])
  R2.HS.kinfit[[kinobjectname]] <- 
    ifelse(class(kinobjects[[kinobjectname]]$fits$HS) == "try-error",
      NA, kinobjects[[kinobjectname]]$results$stats[["HS", "R2"]])
}

EF.SFO.KinGUI <- c(0.9845, 0.9971, 0.9714, 0.9919)
EF.FOMC.KinGUI <- c(0.9831, 0.9973, 0.9955, 0.9920)
EF.HS.KinGUI <- c(0.9972, 0.9972, 0.9980, 0.9945)
EF.DFOP.KinGUI <- c(0.9845, 0.9973, 0.9994, 0.9919)
names(EF.SFO.KinGUI) <- names(EF.FOMC.KinGUI) <- names(kinobjects[1:4])
names(EF.DFOP.KinGUI) <- names(EF.HS.KinGUI) <- names(kinobjects[1:4])

R2 <- data.frame(
  SFO.KinGUI = EF.SFO.KinGUI,
  SFO.kinfit = round(R2.SFO.kinfit, 4),
  FOMC.KinGUI = EF.FOMC.KinGUI,
  FOMC.kinfit = round(R2.FOMC.kinfit, 4),
  DFOP.KinGUI = EF.DFOP.KinGUI,
  DFOP.kinfit = round(R2.DFOP.kinfit, 4),
  HS.KinGUI = EF.HS.KinGUI,
  HS.kinfit = round(R2.HS.kinfit, 4)
)
write.table(R2,
  file = "R2_comparison.tex",
    sep=" & ", quote=FALSE, na="",
    row.names=TRUE, col.names=FALSE, eol = " \\\\ \n")
@

\begin{table}
\caption{Comparison of model efficiency (EF) values calculated 
by KinGUI and $R^2$ values calculated by \Rpackage{kinfit}. \label{tab:vali.R2}}
\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{R2_comparison}
\bottomrule
\end{tabular}
\end{table}

The comparison shows that whenever the comparison was possible, the $R^2$
value calculated by the \Rpackage{kinfit} package was equal to the model 
efficiency calculated by KinGUI, both rounded to four digits.

\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. 

\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=marker: