aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--DESCRIPTION6
-rw-r--r--R/endpoints.R14
-rw-r--r--R/mkinfit.R14
-rw-r--r--R/mkinmod.R1
-rw-r--r--TODO18
-rw-r--r--inst/unitTests/runit.mkinfit.R39
-rw-r--r--man/mkinfit.Rd42
-rw-r--r--vignettes/examples.Rnw6
-rw-r--r--vignettes/examples.pdfbin293187 -> 293678 bytes
-rw-r--r--vignettes/mkin.Rnw4
-rw-r--r--vignettes/mkin.pdfbin163652 -> 163548 bytes
11 files changed, 92 insertions, 52 deletions
diff --git a/DESCRIPTION b/DESCRIPTION
index ed5abdb..868e79e 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-18
-Date: 2013-03-15
+Version: 0.9-19
+Date: 2013-04-14
Author: Johannes Ranke, Katrin Lindenberger, René Lehmann
Maintainer: Johannes Ranke <jranke@uni-bremen.de>
Description: Calculation routines based on the FOCUS Kinetics Report (2006).
@@ -12,7 +12,7 @@ Description: Calculation routines based on the FOCUS Kinetics Report (2006).
and a choice of the optimisation methods made available by the FME package
(default is a Levenberg-Marquardt variant). Please note that no warranty is
implied for correctness of results or fitness for a particular purpose.
-Depends: FME, deSolve
+Depends: FME, deSolve, minpack.lm
Suggests: RUnit
License: GPL
LazyLoad: yes
diff --git a/R/endpoints.R b/R/endpoints.R
index 163dc8d..c9a6a51 100644
--- a/R/endpoints.R
+++ b/R/endpoints.R
@@ -11,6 +11,20 @@ endpoints <- function(fit, pseudoDT50 = FALSE) {
ep$SFORB <- vector()
for (obs_var in obs_vars) {
type = names(fit$mkinmod$map[[obs_var]])[1]
+
+ # Get formation fractions if directly fitted, and calculate remaining fraction to sink
+ f_names = grep(paste("f", obs_var, sep = "_"), names(parms.all), value=TRUE)
+ f_values = parms.all[f_names]
+ f_to_sink = 1 - sum(f_values)
+ names(f_to_sink) = ifelse(type == "SFORB",
+ paste(obs_var, "free", "sink", sep = "_"),
+ paste(obs_var, "sink", sep = "_"))
+ for (f_name in f_names) {
+ ep$ff[[sub("f_", "", sub("_to_", "_", f_name))]] = f_values[[f_name]]
+ }
+ ep$ff = append(ep$ff, f_to_sink)
+
+ # Get the rest
if (type == "SFO") {
k_names = grep(paste("k", obs_var, sep="_"), names(parms.all), value=TRUE)
k_tot = sum(parms.all[k_names])
diff --git a/R/mkinfit.R b/R/mkinfit.R
index b80be97..cf018f0 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -28,6 +28,8 @@ mkinfit <- function(mkinmod, observed,
fixed_parms = NULL,
fixed_initials = names(mkinmod$diffs)[-1],
solution_type = "auto",
+ method.modFit = "Marq",
+ control.modFit = list(),
plot = FALSE, quiet = FALSE,
err = NULL, weight = "none", scaleVar = FALSE,
atol = 1e-8, rtol = 1e-10, n.outtimes = 100,
@@ -45,6 +47,13 @@ mkinfit <- function(mkinmod, observed,
# Define starting values for parameters where not specified by the user
if (parms.ini[[1]] == "auto") parms.ini = vector()
+
+ # Prevent inital parameter specifications that are not in the model
+ wrongpar.names <- setdiff(names(parms.ini), mkinmod$parms)
+ if (length(wrongpar.names) > 0) {
+ stop("Initial parameter(s) ", paste(wrongpar.names, collapse = ", "), " not used in the model")
+ }
+
defaultpar.names <- setdiff(mkinmod$parms, names(parms.ini))
for (parmname in defaultpar.names) {
# Default values for rate constants, depending on the parameterisation
@@ -168,7 +177,7 @@ mkinfit <- function(mkinmod, observed,
}
return(mC)
}
- fit <- modFit(cost, c(state.ini.optim, parms.optim), ...)
+ fit <- modFit(cost, c(state.ini.optim, parms.optim), method = method.modFit, control = control.modFit, ...)
# We need to return some more data for summary and plotting
fit$solution_type <- solution_type
@@ -203,6 +212,7 @@ mkinfit <- function(mkinmod, observed,
fit$data <- data[order(data$variable, data$time), ]
fit$atol <- atol
fit$rtol <- rtol
+
# Return all backtransformed parameters for summary
fit$bparms.optim <- bparms.optim
fit$bparms.fixed <- bparms.fixed
@@ -239,7 +249,7 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05,
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
+ # Only transform boundaries of CI for one parameter at a time
for (pname in pnames) {
par.lower <- par.upper <- object$par
par.lower[pname] <- param[pname, "Lower"]
diff --git a/R/mkinmod.R b/R/mkinmod.R
index adfd9ea..09deea4 100644
--- a/R/mkinmod.R
+++ b/R/mkinmod.R
@@ -26,6 +26,7 @@ mkinmod <- function(..., use_of_ff = "min")
# Check if any of the names of the observed variables contains any other
for (obs_var in obs_vars) {
if (length(grep(obs_var, obs_vars)) > 1) stop("Sorry, variable names can not contain each other")
+ if (grepl("_to_", obs_var)) stop("Sorry, names of observed variables can not contain _to_")
}
if (!use_of_ff %in% c("min", "max"))
diff --git a/TODO b/TODO
index c081ee7..25f5478 100644
--- a/TODO
+++ b/TODO
@@ -1,18 +1,14 @@
-Must have:
-- Some more tests with coupled models with FOMC, DFOP, HS and SFORB for parent
-- Some tests with SFORB for metabolites
-- Make the "method" argument work again - currently it is being passed both to ode() and modFit()
-- Check the contents of the parms.ini argument to mkinfit for arguments that
- are not used in the model in order to avoid erroneous chi2 error level
- calculations
+TODO for version 1.0
+- Transfer (some?) unit tests to examples to make them more easily accessible
+ for the user - they will be tested during R CMD check anyway
+- Let mkin call mkin_wide_to_long automatically, if observed data are suitably defined
+- Think about what a user would expect from version 1.0
Nice to have:
-- Calculate confidence intervals for DT50 and DT90 values when only one parameter is involved
-- Add more unit tests for mkinfit
- Improve choice of starting parameters, in order to avoid failure of eigenvalue based solutions
due to singular matrix
- Consider reverting to deSolve in case eigenvalue based solution fails due to singular matrix
-- Reproduce example anaylses (L1, L2, ...) in FOCUS (2006)
+- Calculate confidence intervals for DT50 and DT90 values when only one parameter is involved
- Calculate transformation only DT50 values (exclude pathways to sink) as additional information
-- Generate confidence intervals for DT50 and DT90 values for other cases
+- Calculate DT50 values from smaller rate constant/eigenvalue of DFOP, HS and SFORB models
- Document assumptions and implications of fitting method used in FME in the package vignette
diff --git a/inst/unitTests/runit.mkinfit.R b/inst/unitTests/runit.mkinfit.R
index 784054b..278560d 100644
--- a/inst/unitTests/runit.mkinfit.R
+++ b/inst/unitTests/runit.mkinfit.R
@@ -295,42 +295,19 @@ test.mkinfit.schaefer07_complex_example <- function()
checkIdentical(r$mkin.deviation < 10, rep(TRUE, length(r$mkin.deviation)))
} # }}}
-# Test deSolve based fit to Schaefer 2007 data against solution from conference paper {{{
-test.mkinfit.schaefer07_complex_example <- function()
+# Test deSolve based fit to Schaefer 2007 using FOMC for parent (eigenvalue based solution not possible)
+# and skipping B1 because of its high scatter {{{
+test.mkinfit.schaefer07_complex_example_2 <- function()
{
- schaefer07_complex_model <- mkinmod(
- parent = list(type = "SFO", to = c("A1", "B1", "C1"), sink = FALSE),
+ schaefer07_complex_model_2 <- mkinmod(
+ parent = list(type = "FOMC", to = c("A1", "C1")),
A1 = list(type = "SFO", to = "A2"),
- B1 = list(type = "SFO"),
C1 = list(type = "SFO"),
A2 = list(type = "SFO"))
- # Works fine with n.outtimes = 1000 but takes too much time
- # fit <- mkinfit(schaefer07_complex_model,
- # mkin_wide_to_long(schaefer07_complex_case, time = "time"),
- # n.outtimes = 1000, solution_type = "deSolve")
- # s <- summary(fit)
- # r <- schaefer07_complex_results
- # attach(as.list(fit$parms.all))
- # k_parent <- sum(k_parent_A1, k_parent_B1, k_parent_C1)
- # r$mkin <- c(
- # k_parent,
- # s$distimes["parent", "DT50"],
- # s$ff["parent_A1"],
- # sum(k_A1_sink, k_A1_A2),
- # s$distimes["A1", "DT50"],
- # s$ff["parent_B1"],
- # k_B1_sink,
- # s$distimes["B1", "DT50"],
- # s$ff["parent_C1"],
- # k_C1_sink,
- # s$distimes["C1", "DT50"],
- # s$ff["A1_A2"],
- # k_A2_sink,
- # s$distimes["A2", "DT50"])
- # r$means <- (r$KinGUI + r$ModelMaker)/2
- # r$mkin.deviation <- abs(round(100 * ((r$mkin - r$means)/r$means), digits=1))
- # checkIdentical(r$mkin.deviation < 10, rep(TRUE, length(r$mkin.deviation)))
+ fit <- mkinfit(schaefer07_complex_model_2,
+ mkin_wide_to_long(schaefer07_complex_case, time = "time"))
+ checkTrue(fit$ssr < 122)
} # }}}
# vim: set foldmethod=marker ts=2 sw=2 expandtab:
diff --git a/man/mkinfit.Rd b/man/mkinfit.Rd
index 319b1b3..445bce2 100644
--- a/man/mkinfit.Rd
+++ b/man/mkinfit.Rd
@@ -15,6 +15,8 @@ mkinfit(mkinmod, observed,
state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)),
fixed_parms = NULL, fixed_initials = names(mkinmod$diffs)[-1],
solution_type = "auto",
+ method.modFit = "Marq",
+ control.modFit = list(),
plot = FALSE, quiet = FALSE, err = NULL, weight = "none",
scaleVar = FALSE,
atol = 1e-8, rtol = 1e-10, n.outtimes = 100,
@@ -38,6 +40,13 @@ mkinfit(mkinmod, observed,
A named vector of initial values for the parameters, including parameters to
be optimised and potentially also fixed parameters as indicated by \code{fixed_parms}.
If set to "auto", initial values for rate constants are set to default values.
+ Using parameter names that are not in the model gives an error.
+
+ It is possible to only specify a subset of the parameters that the model
+ needs. You can use the parameter lists "bparms.ode" from a previously
+ fitted model, which contains the differential equation parameters from this
+ model. This works nicely if the models are nested. An example is given
+ below.
}
\item{state.ini}{
A named vector of initial values for the state variables of the model. In case the
@@ -65,6 +74,15 @@ mkinfit(mkinmod, observed,
using eigenvalues and eigenvectors, and finally "deSolve" for the remaining
models (time dependence of degradation rates and metabolites).
}
+ \item{method.modFit}{
+ The optimisation method passed to \code{\link{modFit}}. The default "Marq" is the Levenberg Marquardt
+ algorithm \code{\link{nls.lm}} from the package \code{minpack.lm}. Often other methods need
+ more iterations to find the same result. When using "Pseudo", "upper" and "lower" need to be
+ specified for the transformed parameters.
+ }
+ \item{control.modFit}{
+ Additional arguments passed to the optimisation method used by \code{\link{modFit}}.
+ }
\item{plot}{
Should the observed values and the numerical solutions be plotted at each stage
of the optimisation?
@@ -120,6 +138,30 @@ SFO_SFO <- mkinmod(
fit <- mkinfit(SFO_SFO, FOCUS_2006_D)
str(fit)
summary(fit)
+
+# Use stepwise fitting, using optimised parameters from parent only fit, FOMC
+\dontrun{
+FOMC <- mkinmod(parent = list(type = "FOMC"))
+FOMC_SFO <- mkinmod(
+ parent = list(type = "FOMC", to = "m1", sink = TRUE),
+ m1 = list(type = "SFO"))
+# Fit the model to the FOCUS example dataset D using defaults
+fit.FOMC_SFO <- mkinfit(FOMC_SFO, FOCUS_2006_D)
+# Use starting parameters from parent only FOMC fit (not really needed in this case)
+fit.FOMC = mkinfit(FOMC, FOCUS_2006_D)
+fit.FOMC_SFO <- mkinfit(FOMC_SFO, FOCUS_2006_D, parms.ini = fit.FOMC$bparms.ode, plot=TRUE)
+}
+
+# Use stepwise fitting, using optimised parameters from parent only fit, SFORB
+SFORB <- mkinmod(parent = list(type = "SFORB"))
+SFORB_SFO <- mkinmod(
+ parent = list(type = "SFORB", to = "m1", sink = TRUE),
+ m1 = list(type = "SFO"))
+# Fit the model to the FOCUS example dataset D using defaults
+fit.SFORB_SFO <- mkinfit(SFORB_SFO, FOCUS_2006_D)
+# Use starting parameters from parent only SFORB fit (not really needed in this case)
+fit.SFORB = mkinfit(SFORB, FOCUS_2006_D)
+fit.SFORB_SFO <- mkinfit(SFORB_SFO, FOCUS_2006_D, parms.ini = fit.SFORB$bparms.ode, plot=TRUE)
}
\keyword{ models }
\keyword{ optimize }
diff --git a/vignettes/examples.Rnw b/vignettes/examples.Rnw
index fa35c8f..4897f7c 100644
--- a/vignettes/examples.Rnw
+++ b/vignettes/examples.Rnw
@@ -433,7 +433,7 @@ Z.mkin.1 <- mkinmod(Z0 = list(type = "SFO", to = "Z1", sink = FALSE),
Z2 = list(type = "SFO", to = "Z3"),
Z3 = list(type = "SFORB"))
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),
+ parms.ini = c(k_Z0_Z1 = 0.5, k_Z1_Z2 = 0.3),
quiet = TRUE)
plot(m.Z.mkin.1)
summary(m.Z.mkin.1, data = FALSE)
@@ -498,8 +498,8 @@ 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.
+clear that nothing can be said about the degradation rate of Z3 towards the end
+of the experiment. However, this appears to be a feature of the data.
<<FOCUS_2006_Z_residuals_11, fig=TRUE>>=
par(mfrow = c(2, 2))
diff --git a/vignettes/examples.pdf b/vignettes/examples.pdf
index 8d32159..bd01239 100644
--- a/vignettes/examples.pdf
+++ b/vignettes/examples.pdf
Binary files differ
diff --git a/vignettes/mkin.Rnw b/vignettes/mkin.Rnw
index 2f79cc5..22c9c15 100644
--- a/vignettes/mkin.Rnw
+++ b/vignettes/mkin.Rnw
@@ -149,9 +149,9 @@ SFO.fit <- mkinfit(SFO, FOCUS_2006_C)
summary(SFO.fit)
SFORB.fit <- mkinfit(SFORB, FOCUS_2006_C)
summary(SFORB.fit)
-SFO_SFO.fit <- mkinfit(SFO_SFO, FOCUS_2006_D, plot=TRUE)
+SFO_SFO.fit <- mkinfit(SFO_SFO, FOCUS_2006_D)
summary(SFO_SFO.fit, data=FALSE)
-SFORB_SFO.fit <- mkinfit(SFORB_SFO, FOCUS_2006_D, plot=TRUE)
+SFORB_SFO.fit <- mkinfit(SFORB_SFO, FOCUS_2006_D)
summary(SFORB_SFO.fit, data=FALSE)
@
diff --git a/vignettes/mkin.pdf b/vignettes/mkin.pdf
index 212feed..3eebf14 100644
--- a/vignettes/mkin.pdf
+++ b/vignettes/mkin.pdf
Binary files differ

Contact - Imprint