diff options
-rw-r--r-- | DESCRIPTION | 6 | ||||
-rw-r--r-- | R/endpoints.R | 14 | ||||
-rw-r--r-- | R/mkinfit.R | 14 | ||||
-rw-r--r-- | R/mkinmod.R | 1 | ||||
-rw-r--r-- | TODO | 18 | ||||
-rw-r--r-- | inst/unitTests/runit.mkinfit.R | 39 | ||||
-rw-r--r-- | man/mkinfit.Rd | 42 | ||||
-rw-r--r-- | vignettes/examples.Rnw | 6 | ||||
-rw-r--r-- | vignettes/examples.pdf | bin | 293187 -> 293678 bytes | |||
-rw-r--r-- | vignettes/mkin.Rnw | 4 | ||||
-rw-r--r-- | vignettes/mkin.pdf | bin | 163652 -> 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"))
@@ -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 Binary files differindex 8d32159..bd01239 100644 --- a/vignettes/examples.pdf +++ b/vignettes/examples.pdf 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 Binary files differindex 212feed..3eebf14 100644 --- a/vignettes/mkin.pdf +++ b/vignettes/mkin.pdf |