- 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

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}}. 
   }