From e3162e617bc268d9de92640311e2fbe650aa636a Mon Sep 17 00:00:00 2001
From: jranke <jranke@edb9625f-4e0d-4859-8d74-9fd3b1da38cb>
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