From 6168089ac43664c10ca2dc1281c8648fbf3b35a9 Mon Sep 17 00:00:00 2001 From: jranke Date: Sun, 14 Apr 2013 12:42:06 +0000 Subject: - Update the TODO list, setting some requirements for version 1.0 - Check that initial values specified using parms.ini are actually needed for the model, stop otherwise - List all formation fractions in the same place in the summary, also if they were fitted in the model - Include an FOMC model coupled to two metabolites in the unit tests - Some updates needed because of the above - Update of static documentation including the vignettes - Update of the mkin vignettes in the vignettes directory git-svn-id: svn+ssh://svn.r-forge.r-project.org/svnroot/kinfit/pkg/mkin@82 edb9625f-4e0d-4859-8d74-9fd3b1da38cb --- DESCRIPTION | 6 +++--- R/endpoints.R | 14 ++++++++++++++ R/mkinfit.R | 14 ++++++++++++-- R/mkinmod.R | 1 + TODO | 18 +++++++----------- inst/unitTests/runit.mkinfit.R | 39 ++++++++------------------------------ man/mkinfit.Rd | 42 +++++++++++++++++++++++++++++++++++++++++ vignettes/examples.Rnw | 6 +++--- vignettes/examples.pdf | Bin 293187 -> 293678 bytes vignettes/mkin.Rnw | 4 ++-- 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 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. <>= par(mfrow = c(2, 2)) diff --git a/vignettes/examples.pdf b/vignettes/examples.pdf index 8d32159..bd01239 100644 Binary files a/vignettes/examples.pdf and b/vignettes/examples.pdf 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 Binary files a/vignettes/mkin.pdf and b/vignettes/mkin.pdf differ -- cgit v1.2.1