- Calculate confidence intervals for parameters based on the t distribution

  Thanks to the CAKE developers at Tessella for the call to qt()
- Calculate asymetric confidence intervals for backtransformed parameters
- Overhaul and recompile vignettes



git-svn-id: svn+ssh://svn.r-forge.r-project.org/svnroot/kinfit/pkg/mkin@74 edb9625f-4e0d-4859-8d74-9fd3b1da38cb

7 files changed, 134 insertions, 91 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index c4ae84d..431df8b 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -2,8 +2,8 @@ Package: mkin
 Type: Package
 Title: Routines for fitting kinetic models with one or more state
   variables to chemical degradation data
-Version: 0.9-15
-Date: 2013-03-04
+Version: 0.9-16
+Date: 2013-03-05
 Author: Johannes Ranke, Katrin Lindenberger, René Lehmann
 Maintainer: Johannes Ranke <jranke@uni-bremen.de>
 Description: Calculation routines based on the FOCUS Kinetics Report (2006).
diff --git a/R/mkinfit.R b/R/mkinfit.R
index 1061165..12a67bc 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -213,30 +213,42 @@ mkinfit <- function(mkinmod, observed,
   return(fit)

 }

 

-summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, ...) {

+summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, ...) {

   param  <- object$par

   pnames <- names(param)

   p      <- length(param)

+  mod_vars <- names(object$mkinmod$diffs)

   covar  <- try(solve(0.5*object$hessian), silent = TRUE)   # unscaled covariance

-  if (!is.numeric(covar)) {

-    message <- "Cannot estimate covariance; system is singular"

-    warning(message)

-    covar <- matrix(data = NA, nrow = p, ncol = p)

-  } else message <- "ok"

-

-  rownames(covar) <- colnames(covar) <- pnames

   rdf    <- object$df.residual

   resvar <- object$ssr / rdf

-  se     <- sqrt(diag(covar) * resvar)

+  if (!is.numeric(covar)) {

+    covar <- NULL

+    se <- lci <- uci <- rep(NA, p)

+  } else {

+    rownames(covar) <- colnames(covar) <- pnames

+    se     <- sqrt(diag(covar) * resvar)

+    lci    <- param + qt(alpha/2, rdf) * se

+    uci    <- param + qt(1-alpha/2, rdf) * se

+

+  }

+

   names(se) <- pnames

-  tval      <- param / se

   modVariance <- object$ssr / length(object$residuals)

 

-  param <- cbind(param, se)

-  dimnames(param) <- list(pnames, c("Estimate", "Std. Error"))

-

-  bparam <- as.matrix(object$bparms.optim)

-  dimnames(bparam) <- list(pnames, c("Estimate"))

+  param <- cbind(param, se, lci, uci)

+  dimnames(param) <- list(pnames, c("Estimate", "Std. Error", "Lower", "Upper"))

+

+  blci <- buci <- numeric()

+  # Only use lower end of CI for one parameter at a time

+  for (pname in pnames) {

+    par.lower <- par.upper <- object$par

+    par.lower[pname] <- param[pname, "Lower"]

+    par.upper[pname] <- param[pname, "Upper"]

+    blci[pname] <- backtransform_odeparms(par.lower, mod_vars)[pname]

+    buci[pname] <- backtransform_odeparms(par.upper, mod_vars)[pname]

+  }

+  bparam <- cbind(object$bparms.optim, blci, buci)

+  dimnames(bparam) <- list(pnames, c("Estimate", "Lower", "Upper"))

 

   ans <- list(

     version = as.character(packageVersion("mkin")),

@@ -248,9 +260,11 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, ...) {
     residualVariance = resvar,

     sigma = sqrt(resvar),

     modVariance = modVariance,

-    df = c(p, rdf), cov.unscaled = covar,

+    df = c(p, rdf), 

+    cov.unscaled = covar,

     cov.scaled = covar * resvar,

-    info = object$info, niter = object$iterations,

+    info = object$info, 

+    niter = object$iterations,

     stopmess = message,

     par = param,

     bpar = bparam)

@@ -293,10 +307,10 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), .
   else print(x$fixed)

   

   cat("\nOptimised, transformed parameters:\n")

-  printCoefmat(x$par, digits = digits, ...)

+  print(signif(x$par, digits = digits))

 

   cat("\nBacktransformed parameters:\n")

-  printCoefmat(x$bpar, digits = digits, ...)

+  print(signif(x$bpar, digits = digits))

 

   cat("\nResidual standard error:",

       format(signif(x$sigma, digits)), "on", rdf, "degrees of freedom\n")

@@ -323,12 +337,13 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), .
     print(x$SFORB, digits=digits,...)

   }

 

-  printcor <- is.numeric(x$cov.unscaled)

-  if (printcor){

+  cat("\nParameter correlation:\n")

+  if (!is.null(x$cov.unscaled)){

     Corr <- cov2cor(x$cov.unscaled)

     rownames(Corr) <- colnames(Corr) <- rownames(x$par)

-    cat("\nParameter correlation:\n")

     print(Corr, digits = digits, ...)

+  } else {

+    cat("Could not estimate covariance matrix; singular system:\n")

   }

 

   printdata <- !is.null(x$data)

diff --git a/TODO b/TODO
index 0d42d78..95c1476 100644
--- a/TODO
+++ b/TODO
@@ -7,7 +7,6 @@ Must have:
   calculations
 
 Nice to have:
-- Calculate confidence intervals for parameters assuming normal distribution
 - Calculate confidence intervals for DT50 and DT90 values when only one parameter is involved
 - Add more unit tests for mkinfit
 - Document validation against fits documented in chapter 13 of FOCUS (2006)
diff --git a/man/summary.mkinfit.Rd b/man/summary.mkinfit.Rd
index 59b988a..472b5de 100644
--- a/man/summary.mkinfit.Rd
+++ b/man/summary.mkinfit.Rd
@@ -11,7 +11,7 @@
   and residual values.

 }

 \usage{

-\method{summary}{mkinfit}(object, data = TRUE, distimes = TRUE, ...)

+\method{summary}{mkinfit}(object, data = TRUE, distimes = TRUE, alpha = 0.05, ...)

 \method{print}{summary.mkinfit}(x, digits = max(3, getOption("digits") - 3), ...)

 }

 

@@ -28,6 +28,9 @@
   \item{distimes}{

   logical, indicating whether DT50 and DT90 values should be included.

 }

+  \item{alpha}{

+  error level for confidence interval estimation from t distribution

+}

   \item{digits}{

   Number of digits to use for printing

 }

@@ -36,15 +39,24 @@
 }

 }

 \value{

-  The summary function returns a list with the same components as

-  \code{\link{summary.modFit}}, and the additional components

+  The summary function returns a list derived from

+  \code{\link{summary.modFit}}, with components, among others

+  \item{version, Rversion}{The mkin and R versions used}

+  \item{date.fit, date.summary}{The dates where the fit and the summary were produced}

+  \item{use_of_ff}{Was maximum or minimum use made of formation fractions}

+  \item{residuals, residualVariance, sigma, modVariance, df}{As in summary.modFit}

+  \item{cov.unscaled, cov.scaled, info, niter, stopmess, par}{As in summary.modFit}

+  \item{bpar}{Optimised and backtransformed parameters}

   \item{diffs }{The differential equations used in the model}

   \item{data }{The data (see Description above).}

   \item{start }{The starting values and bounds, if applicable, for optimised parameters.}

   \item{fixed }{The values of fixed parameters.}

   \item{errmin }{The chi2 error levels for each observed variable.}

-  \item{disstimes }{The DT50 and DT90 values for each observed variable.}

+  \item{bparms.ode }{All backtransformed ODE parameters, for use as starting parameters for

+    related models.}

   \item{ff }{The estimated formation fractions derived from the fitted model.}

+  \item{distimes }{The DT50 and DT90 values for each observed variable.}

+  \item{SFORB}{If applicable, eigenvalues of SFORB components of the model.}

   The print method is called for its side effect, i.e. printing the summary.

 }

 \references{

diff --git a/vignettes/examples.Rnw b/vignettes/examples.Rnw
index 6f3cfc9..fa35c8f 100644
--- a/vignettes/examples.Rnw
+++ b/vignettes/examples.Rnw
@@ -93,13 +93,13 @@ summary(m.L1.SFO)
 

 A plot of the fit is obtained with the plot function for mkinfit objects.

 

-<<L1_SFO_plot, fig=TRUE, echo=TRUE>>=

+<<L1_SFO_plot, fig=TRUE, echo=TRUE, height=4>>=

 plot(m.L1.SFO)

 @

 

 The residual plot can be easily obtained by

 

-<<L1_SFO_residuals, fig=TRUE, echo=TRUE>>=

+<<L1_SFO_residuals, fig=TRUE, echo=TRUE, height=4>>=

 mkinresplot(m.L1.SFO, ylab = "Observed", xlab = "Time")

 @

 

@@ -113,7 +113,8 @@ 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\%).

+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

@@ -167,22 +168,22 @@ models generally only implement SFO kinetics.
 For comparison, the FOMC model is fitted as well, and the $\chi^2$ error level

 is checked.

 

-<<L2_FOMC, echo=TRUE, fig=TRUE>>=

-m.L2.FOMC <- mkinfit(FOMC, FOCUS_2006_L2_mkin, quiet=TRUE)

+<<L2_FOMC, echo=TRUE, fig=TRUE, height=8>>=

+m.L2.FOMC <- mkinfit(FOMC, FOCUS_2006_L2_mkin, quiet = TRUE)

+par(mfrow = c(2, 1))

 plot(m.L2.FOMC)

-s.m.L2.FOMC <- summary(m.L2.FOMC)

-s.m.L2.FOMC$errmin

+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 does not further reduce the 

-$\chi^2$ error level. 

+Fitting the four parameter DFOP model further reduces the $\chi^2$ error level. 

 

-<<L2_DFOP, echo=TRUE, fig=TRUE>>=

-m.L2.DFOP <- mkinfit(DFOP, FOCUS_2006_L2_mkin, quiet=TRUE)

+<<L2_DFOP, echo=TRUE, fig=TRUE, height=4>>=

+m.L2.DFOP <- mkinfit(DFOP, FOCUS_2006_L2_mkin, quiet = TRUE)

 plot(m.L2.DFOP)

 @

 

@@ -190,17 +191,18 @@ 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.

 

-<<L2_DFOP_2, echo=TRUE, fig=TRUE>>=

+<<L2_DFOP_2, echo=TRUE, fig=TRUE, height=4>>=

 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)

-s.m.L2.DFOP <- summary(m.L2.DFOP)

-s.m.L2.DFOP$errmin

+summary(m.L2.DFOP, data = FALSE)

 @

 

-Therefore, the FOMC model is clearly the best-fit model for dataset L1 based on the 

-$\chi^2$ error level criterion.

+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}

 

@@ -216,10 +218,10 @@ FOCUS_2006_L3_mkin <- mkin_wide_to_long(FOCUS_2006_L3)
 

 SFO model, summary and plot:

 

-<<L3_SFO, echo=TRUE, fig=TRUE>>=

+<<L3_SFO, echo=TRUE, fig=TRUE, height=4>>=

 m.L3.SFO <- mkinfit(SFO, FOCUS_2006_L3_mkin, quiet = TRUE)

-summary(m.L3.SFO)

 plot(m.L3.SFO)

+summary(m.L3.SFO)

 @

 

 The $\chi^2$ error level of 22\% as well as the plot suggest that the model

@@ -227,12 +229,10 @@ does not fit very well.
 

 The FOMC model performs better:

 

-<<L3_FOMC, echo=TRUE, fig=TRUE>>=

+<<L3_FOMC, echo=TRUE, fig=TRUE, height=4>>=

 m.L3.FOMC <- mkinfit(FOMC, FOCUS_2006_L3_mkin, quiet = TRUE)

 plot(m.L3.FOMC)

-s.m.L3.FOMC <- summary(m.L3.FOMC)

-s.m.L3.FOMC$errmin

-endpoints(m.L3.FOMC)

+summary(m.L3.FOMC, data = FALSE)

 @

 

 The error level at which the $\chi^2$ test passes is 7\% in this case.

@@ -240,14 +240,15 @@ 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:

 

-<<L3_DFOP, echo=TRUE, fig=TRUE>>=

+<<L3_DFOP, echo=TRUE, fig=TRUE, height=4>>=

 m.L3.DFOP <- mkinfit(DFOP, FOCUS_2006_L3_mkin, quiet = TRUE)

 plot(m.L3.DFOP)

-s.m.L3.DFOP <- summary(m.L3.DFOP)

-s.m.L3.DFOP$errmin

+summary(m.L3.DFOP, data = FALSE)

 @

 

-Therefore, the DFOP model is the best-fit model based on the $\chi^2$ error

+Here, a look to the model plot, the confidence intervals of the parameters 

+and the correlation matrix suggest that the paramter 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}

@@ -264,10 +265,10 @@ FOCUS_2006_L4_mkin <- mkin_wide_to_long(FOCUS_2006_L4)
 

 SFO model, summary and plot:

 

-<<L4_SFO, echo=TRUE, fig=TRUE>>=

+<<L4_SFO, echo=TRUE, fig=TRUE, height=4>>=

 m.L4.SFO <- mkinfit(SFO, FOCUS_2006_L4_mkin, quiet = TRUE)

-summary(m.L4.SFO)

 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

@@ -275,11 +276,10 @@ fits very well.
 

 The FOMC model for comparison

 

-<<L4_FOMC, echo=TRUE, fig=TRUE>>=

+<<L4_FOMC, echo=TRUE, fig=TRUE, height=4>>=

 m.L4.FOMC <- mkinfit(FOMC, FOCUS_2006_L4_mkin, quiet = TRUE)

 plot(m.L4.FOMC)

-s.m.L4.FOMC <- summary(m.L4.FOMC)

-s.m.L4.FOMC$errmin

+summary(m.L4.FOMC, data = FALSE)

 @

 

 The error level at which the $\chi^2$ test passes is slightly lower for the FOMC 

@@ -315,29 +315,29 @@ 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).

 

-<<FOCUS_2006_Z_fits_1, echo=TRUE, fig=TRUE>>=

+<<FOCUS_2006_Z_fits_1, echo=TRUE, fig=TRUE, height=4>>=

 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)

-summary(m.Z.2a, data = FALSE)

 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 to the analysis in the FOCUS kinetics report, to simplify 

+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:

 

-<<FOCUS_2006_Z_fits_2, echo=TRUE, fig=TRUE>>=

+<<FOCUS_2006_Z_fits_2, echo=TRUE, fig=TRUE, height=4>>=

 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)

-summary(m.Z.2a.ff, data = FALSE)

 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, 

@@ -348,18 +348,17 @@ 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.

 

-<<FOCUS_2006_Z_fits_3, echo=TRUE, fig=TRUE>>=

+<<FOCUS_2006_Z_fits_3, echo=TRUE, fig=TRUE, height=4>>=

 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", 

-                 quiet = TRUE)

-summary(m.Z.3, data = FALSE)

+#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 fails, as the default solution type chosen

+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

@@ -367,36 +366,36 @@ 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.

 

-<<FOCUS_2006_Z_fits_4, echo=TRUE, fig=TRUE>>=

+<<FOCUS_2006_Z_fits_4, echo=TRUE, fig=TRUE, height=4>>=

 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)

-summary(m.Z.4a, data = FALSE)

 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. While this step appears questionable on the basis of the above results, it 

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

 

-<<FOCUS_2006_Z_fits_5, echo=TRUE, fig=TRUE>>=

+<<FOCUS_2006_Z_fits_5, echo=TRUE, fig=TRUE, height=4>>=

 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)

-summary(m.Z.5, data = FALSE)

 plot(m.Z.5)

+summary(m.Z.5, data = FALSE)

 @

 

 Finally, metabolite Z3 is added to the model.

 

-<<FOCUS_2006_Z_fits_6, echo=TRUE, fig=TRUE>>=

+<<FOCUS_2006_Z_fits_6, echo=TRUE, fig=TRUE, height=4>>=

 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"),

@@ -404,8 +403,8 @@ Z.FOCUS <- mkinmod(Z0 = list(type = "SFO", to = "Z1", sink = FALSE),
 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)

-summary(m.Z.FOCUS, data = FALSE)

 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 

@@ -428,7 +427,7 @@ 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.

 

-<<FOCUS_2006_Z_fits_7, echo=TRUE, fig=TRUE>>=

+<<FOCUS_2006_Z_fits_7, echo=TRUE, fig=TRUE, height=4>>=

 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"),

@@ -436,46 +435,49 @@ Z.mkin.1 <- mkinmod(Z0 = list(type = "SFO", to = "Z1", sink = FALSE),
 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, k_Z2_Z3 = 0.2), 

                   quiet = TRUE)

-summary(m.Z.mkin.1, data = FALSE)

 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.

 

-<<FOCUS_2006_Z_fits_8, echo=TRUE, fig=TRUE>>=

+<<FOCUS_2006_Z_fits_8, echo=TRUE, fig=TRUE, height=4>>=

 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)

-summary(m.Z.mkin.2, data = FALSE)

 plot(m.Z.mkin.2)

+summary(m.Z.mkin.2, data = FALSE)

 @

 

-The sink is for Z1 is turned off again, for the same reasons as in the original analysis.

-Then, metabolite Z2 is added.

+When metabolite Z2 is added, the additional sink for Z1 is turned off again,

+for the same reasons as in the original analysis.

 

-<<FOCUS_2006_Z_fits_9, echo=TRUE, fig=TRUE>>=

+<<FOCUS_2006_Z_fits_9, echo=TRUE, fig=TRUE, height=4>>=

 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)

-summary(m.Z.mkin.3, data = FALSE)

 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.

 

-<<FOCUS_2006_Z_fits_10, echo=TRUE, fig=TRUE>>=

+<<FOCUS_2006_Z_fits_10, echo=TRUE, fig=TRUE, height=4>>=

 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)

-summary(m.Z.mkin.4, data = FALSE)

 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

@@ -483,20 +485,35 @@ 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.

-Therefore, the model \texttt{Z.mkin.5} is proposed as the best-fit model

-for the dataset from Appendix 7 of the FOCUS report.

 

-<<FOCUS_2006_Z_fits_11, echo=TRUE, fig=TRUE>>=

+<<FOCUS_2006_Z_fits_11, echo=TRUE, fig=TRUE, height=4>>=

 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)

-summary(m.Z.mkin.5, data = FALSE)

 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 a degradation rate of Z3. However, this

+appears to be a feature of the data.

+

+<<FOCUS_2006_Z_residuals_11, fig=TRUE>>=

+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}

 

diff --git a/vignettes/examples.pdf b/vignettes/examples.pdf
index 3bd0c37..f4cfb31 100644
--- a/vignettes/examples.pdf
+++ b/vignettes/examples.pdf
diff --git a/vignettes/mkin.pdf b/vignettes/mkin.pdf
index 073e361..227ea79 100644
--- a/vignettes/mkin.pdf
+++ b/vignettes/mkin.pdf