From e3162e617bc268d9de92640311e2fbe650aa636a Mon Sep 17 00:00:00 2001 From: jranke Date: Mon, 23 Apr 2012 06:36:24 +0000 Subject: - Some fixes to mkinfit, to account for the changes in mkinmod - See TODO for what remains to be done (not complete) git-svn-id: svn+ssh://svn.r-forge.r-project.org/svnroot/kinfit/pkg/mkin@29 edb9625f-4e0d-4859-8d74-9fd3b1da38cb --- R/mkinfit.R | 43 +++++++++++++++++++++++++------------------ TODO | 2 +- man/mkinfit.Rd | 23 +++++++++++++++++------ 3 files changed, 43 insertions(+), 25 deletions(-) diff --git a/R/mkinfit.R b/R/mkinfit.R index 64edaba..37eee33 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -26,25 +26,22 @@ mkinfit <- function(mkinmod, observed, state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)), fixed_parms = NULL, fixed_initials = names(mkinmod$diffs)[-1], - eigen = FALSE, + solution_type = "auto", plot = FALSE, quiet = FALSE, err = NULL, weight = "none", scaleVar = FALSE, - atol = 1e-6, + atol = 1e-6, n.outtimes = 100, ...) { # Get the names of the state variables in the model mod_vars <- names(mkinmod$diffs) - # See which variant of the model specification was used - use_of_ff <- mkinmod$use_of_ff - # Subset observed data with names of observed data in the model observed <- subset(observed, name %in% names(mkinmod$map)) # Get names of observed variables obs_vars = unique(as.character(observed$name)) - # Define initial values for parameters where not specified by the user + # Define starting values for parameters where not specified by the user if (parms.ini[[1]] == "auto") parms.ini = vector() defaultpar.names <- setdiff(mkinmod$parms, names(parms.ini)) for (parmname in defaultpar.names) { @@ -54,7 +51,7 @@ mkinfit <- function(mkinmod, observed, if (grepl("free_bound$", parmname)) parms.ini[parmname] = 0.1 if (grepl("bound_free$", parmname)) parms.ini[parmname] = 0.02 # Default values for formation fractions - if (substr(parmname, 1, 2) == "f_") parms.ini[parmname] = 0.1 + if (substr(parmname, 1, 2) == "f_") parms.ini[parmname] = 0.2 # Default values for the FOMC, DFOP and HS models if (parmname == "alpha") parms.ini[parmname] = 1 if (parmname == "beta") parms.ini[parmname] = 10 @@ -88,13 +85,21 @@ mkinfit <- function(mkinmod, observed, # Decide if the solution of the model can be based on a simple analytical # formula, the spectral decomposition of the matrix (fundamental system) # or a numeric ode solver from the deSolve package - if (length(mkinmod$map) == 1) { - solution = "analytical" - } else { - if (is.matrix(mkinmod$coefmat) && eigen) { - solution = "eigen" + if (!solution_type %in% c("auto", "analytical", "eigen", "deSolve")) + stop("solution_type must be auto, analytical, eigen or de Solve") + if (solution_type == "analytical" && length(mkinmod$map) > 1) + stop("Analytical solution not implemented for models with metabolites.") + if (solution_type == "eigen" && !is.matrix(mkinmod$coefmat)) + stop("Eigenvalue based solution not possible, coefficient matrix not present.") + if (solution_type == "auto") { + if (length(mkinmod$map) == 1) { + solution_type = "analytical" } else { - solution = "deSolve" + if (is.matrix(mkinmod$coefmat)) { + solution_type = "eigen" + } else { + solution_type = "deSolve" + } } } @@ -107,7 +112,9 @@ mkinfit <- function(mkinmod, observed, assign("calls", calls+1, inherits=TRUE) # Increase the model solution counter # Time points at which observed data are available - outtimes = unique(observed$time) + # Make sure we include time 0, so initial values for state variables are for time 0 + outtimes = sort(unique(c(observed$time, + seq(min(observed$time), max(observed$time), length.out=n.outtimes)))) if(length(state.ini.optim) > 0) { odeini <- c(P[1:length(state.ini.optim)], state.ini.fixed) @@ -119,7 +126,7 @@ mkinfit <- function(mkinmod, observed, parms <- backtransform_odeparms(odeparms, mod_vars) # Solve the system with current transformed parameter values - out <- mkinpredict(mkinmod, parms, odeini, outtimes, solution_type = solution, ...) + out <- mkinpredict(mkinmod, parms, odeini, outtimes, solution_type = solution_type, ...) assign("out_predicted", out, inherits=TRUE) @@ -135,7 +142,7 @@ mkinfit <- function(mkinmod, observed, outtimes_plot = seq(min(observed$time), max(observed$time), length.out=100) out_plot <- mkinpredict(mkinmod, parms, odeini, outtimes_plot, - solution_type = solution, ...) + solution_type = solution_type, ...) plot(0, type="n", xlim = range(observed$time), ylim = range(observed$value, na.rm=TRUE), @@ -158,8 +165,8 @@ mkinfit <- function(mkinmod, observed, fit <- modFit(cost, c(state.ini.optim, parms.optim), ...) # We need to return some more data for summary and plotting - fit$solution <- solution - if (solution == "eigen") { + fit$solution_type <- solution_type + if (solution_type == "eigen") { fit$coefmat <- mkinmod$coefmat } diff --git a/TODO b/TODO index 134a4a1..10ccf91 100644 --- a/TODO +++ b/TODO @@ -1,4 +1,4 @@ -- Debug mkinmod +- Fix transformation and backtransformation of formation fractions - Adapt mkinplot function - Calculate confidence intervals for parameters assuming normal distribution - Calculate confidence intervals for DT50 and DT90 values when only one parameter is involved diff --git a/man/mkinfit.Rd b/man/mkinfit.Rd index 17e0a11..0c8e48f 100644 --- a/man/mkinfit.Rd +++ b/man/mkinfit.Rd @@ -14,10 +14,10 @@ mkinfit(mkinmod, observed, parms.ini = "auto", state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)), fixed_parms = NULL, fixed_initials = names(mkinmod$diffs)[-1], - eigen = FALSE, + solution_type = "auto", plot = FALSE, quiet = FALSE, err = NULL, weight = "none", scaleVar = FALSE, - atol = 1e-6, ...) + atol = 1e-6, n.outtimes, ...) } \arguments{ \item{mkinmod}{ @@ -53,11 +53,16 @@ mkinfit(mkinmod, observed, The names of model variables for which the initial state at time 0 should be excluded from the optimisation. Defaults to all state variables except for the first one. } - \item{eigen}{ - Should the solution of the system of differential equations be based on the + \item{solution_type}{ + If set to "eigen", the solution of the system of differential equations is based on the spectral decomposition of the coefficient matrix in cases that this is - possible? Be aware that the results may differ from the results returned using - the ode solver. + possible. If set to "deSolve", a numerical ode solver from package + \code{\link{deSolve}} is used. If set to "analytical", an analytical solution + of the model is used. This is only implemented for simple degradation experiments + with only one state variable, i.e. with no metabolites. The default is "auto", + which uses "analytical" if possible, otherwise "eigen" if the model can be expressed + using eigenvalues and eigenvectors, and finally "deSolve" for the remaining + models (time dependence of degradation rates and metabolites). } \item{plot}{ Should the observed values and the numerical solutions be plotted at each stage @@ -81,6 +86,12 @@ mkinfit(mkinmod, observed, Absolute error tolerance, passed to \code{\link{ode}}. Default is 1e-6 as in \code{\link{lsoda}}. } + \item{n.outtimes}{ + The length of the dataseries that is produced by the model prediction + function \code{\link{mkinpredict}}. This impacts the accuracy of + the numerical solver if that is used (see \code{solution} argument. + The default value is 100. + } \item{\dots}{ Further arguments that will be passed to \code{\link{modFit}}. } -- cgit v1.2.1