diff options
| -rw-r--r-- | DESCRIPTION | 4 | ||||
| -rw-r--r-- | R/mkinfit.R | 59 | ||||
| -rw-r--r-- | TODO | 1 | ||||
| -rw-r--r-- | man/summary.mkinfit.Rd | 20 | ||||
| -rw-r--r-- | vignettes/examples.Rnw | 141 | ||||
| -rw-r--r-- | vignettes/examples.pdf | bin | 281843 -> 293191 bytes | |||
| -rw-r--r-- | vignettes/mkin.pdf | bin | 162843 -> 134756 bytes |
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)
@@ -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 Binary files differindex 3bd0c37..f4cfb31 100644 --- a/vignettes/examples.pdf +++ b/vignettes/examples.pdf diff --git a/vignettes/mkin.pdf b/vignettes/mkin.pdf Binary files differindex 073e361..227ea79 100644 --- a/vignettes/mkin.pdf +++ b/vignettes/mkin.pdf |