- Separated model prediction out into a separate function

- Created separate functions for parameter transformations (not documented yet)
- Fitting of analytical models works
- SFO_SFO model works, complex models do not at the moment
- summary.mkinfit not operational at the moment



git-svn-id: svn+ssh://svn.r-forge.r-project.org/svnroot/kinfit/pkg/mkin@21 edb9625f-4e0d-4859-8d74-9fd3b1da38cb

5 files changed, 183 insertions, 188 deletions

diff --git a/R/mkinfit.R b/R/mkinfit.R
index 83896bc..47c39cf 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -32,9 +32,12 @@ mkinfit <- function(mkinmod, observed,
   atol = 1e-6,

   ...)

 {

+  # Get the names of the state variables in the model

   mod_vars <- names(mkinmod$diffs)

-  # Subset dataframe of observed variables with mapped (modelled) variables

+

+  # 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))

 

@@ -42,26 +45,33 @@ mkinfit <- function(mkinmod, observed,
   if (parms.ini[[1]] == "auto") parms.ini = vector()

   defaultpar.names <- setdiff(mkinmod$parms, names(parms.ini))

   for (parmname in defaultpar.names) {

-    # Default values for the FOMC model

+    # Default values for rate constants, depending on the parameterisation

+    if (substr(parmname, 1, 2) == "k_") parms.ini[parmname] = 0.1 

+    # Default values for formation fractions

+    if (substr(parmname, 1, 2) == "f_") parms.ini[parmname] = 0.1

+    # Default values for the FOMC, DFOP and HS models

     if (parmname == "alpha") parms.ini[parmname] = 1

     if (parmname == "beta") parms.ini[parmname] = 10

-    # Default values for log of rate constants, depending on the parameterisation

-    if (substr(parmname, 1, 2) == "k_") parms.ini[parmname] = -2.3

-    # Default values for ilr of formation fractions

-    if (substr(parmname, 1, 2) == "f_") parms.ini[parmname] = 0.1

+    if (parmname == "k1") parms.ini[parmname] = 0.1

+    if (parmname == "k2") parms.ini[parmname] = 0.01

+    if (parmname == "tb") parms.ini[parmname] = 5

+    if (parmname == "g") parms.ini[parmname] = 0.5

   }

+  parms.ini <- transform_odeparms(parms.ini, mod_vars)

    

-  # Name the inital parameter values if they are not named yet

+  # Name the inital state variable values if they are not named yet

   if(is.null(names(state.ini))) names(state.ini) <- mod_vars

 

-  # Parameters to be optimised

+  # Parameters to be optimised:

+  # Kinetic parameters in parms.ini whose names are not in fixed_parms

   parms.fixed <- parms.ini[fixed_parms]

-  optim_parms <- setdiff(names(parms.ini), fixed_parms)

-  parms.optim <- parms.ini[optim_parms]

+  parms.optim <- parms.ini[setdiff(names(parms.ini), fixed_parms)]

 

+  # Inital state variables in state.ini whose names are not in fixed_initials

   state.ini.fixed <- state.ini[fixed_initials]

-  optim_initials <- setdiff(names(state.ini), fixed_initials)

-  state.ini.optim <- state.ini[optim_initials]

+  state.ini.optim <- state.ini[setdiff(names(state.ini), fixed_initials)]

+

+  # Preserve names of state variables before renaming initial state variable parameters

   state.ini.optim.boxnames <- names(state.ini.optim)

   if(length(state.ini.optim) > 0) {

       names(state.ini.optim) <- paste(names(state.ini.optim), "0", sep="_")

@@ -77,42 +87,17 @@ mkinfit <- function(mkinmod, observed,
     else solution = "deSolve"

   }

 

-  # Create a function calculating the differentials specified by the model

-  # if necessary

-  if(solution == "deSolve") {

-    mkindiff <- function(t, state, parms) {

-      transparms <- vector()

-      transparms <- c(transparms, exp(parms[grep("^k_", names(parms))]))

-      for (box in mod_vars) {

-        path_indices <- grep(paste("^f", box, sep = "_"), names(parms))

-        n_paths <- length(path_indices)

-        if (n_paths > 0) {

-          f <- invilr(parms[grep(paste("^f", box, sep="_"), names(parms))])

-          transparms <- c(transparms, f[1:n_paths])

-        }

-      }

-      otherparms <- parms[setdiff(names(parms), names(transparms))]

-      parms <- c(transparms, otherparms)

-

-      time <- t

-      diffs <- vector()

-      for (box in mod_vars)

-      {

-        diffname <- paste("d", box, sep="_")      

-        diffs[diffname] <- with(as.list(c(time, state, parms)),

-          eval(parse(text=mkinmod$diffs[[box]])))

-      }

-      return(list(c(diffs)))

-    } 

-  }

-

-  cost.old <- 1e100

-  calls <- 0

+  cost.old <- 1e100 # The first model cost should be smaller than this value

+  calls <- 0 # Counter for number of model solutions

   out_predicted <- NA

   # Define the model cost function

   cost <- function(P)

   {

-    assign("calls", calls+1, inherits=TRUE)

+    assign("calls", calls+1, inherits=TRUE) # Increase the model solution counter

+

+    # Time points at which observed data are available

+    outtimes = unique(observed$time)

+

     if(length(state.ini.optim) > 0) {

       odeini <- c(P[1:length(state.ini.optim)], state.ini.fixed)

       names(odeini) <- c(state.ini.optim.boxnames, names(state.ini.fixed))

@@ -120,76 +105,14 @@ mkinfit <- function(mkinmod, observed,
 

     odeparms <- c(P[(length(state.ini.optim) + 1):length(P)], parms.fixed)

 

-    outtimes = unique(observed$time)

-    evalparse <- function(string)

-    {

-      eval(parse(text=string), as.list(c(odeparms, odeini)))

-    }

+    transparms <- transform_odeparms(odeparms, mod_vars)

 

-    # Solve the system

-    if (solution == "analytical") {

-      parent.type = names(mkinmod$map[[1]])[1]  

-      parent.name = names(mkinmod$diffs)[[1]]

-      o <- switch(parent.type,

-        SFO = SFO.solution(outtimes, 

-            evalparse(parent.name),

-            evalparse(paste("k", parent.name, "sink", sep="_"))),

-        FOMC = FOMC.solution(outtimes,

-            evalparse(parent.name),

-            evalparse("alpha"), evalparse("beta")),

-        DFOP = DFOP.solution(outtimes,

-            evalparse(parent.name),

-            evalparse("k1"), evalparse("k2"),

-            evalparse("g")),

-        HS = HS.solution(outtimes,

-            evalparse(parent.name),

-            evalparse("k1"), evalparse("k2"),

-            evalparse("tb")),

-        SFORB = SFORB.solution(outtimes,

-            evalparse(parent.name),

-            evalparse(paste("k", parent.name, "bound", sep="_")),

-            evalparse(paste("k", sub("free", "bound", parent.name), "free", sep="_")),

-            evalparse(paste("k", parent.name, "sink", sep="_")))

-      )

-      out <- cbind(outtimes, o)

-      dimnames(out) <- list(outtimes, c("time", sub("_free", "", parent.name)))

-    }

-    if (solution == "eigen") {

-      coefmat.num <- matrix(sapply(as.vector(mkinmod$coefmat), evalparse), 

-        nrow = length(mod_vars))

-      e <- eigen(coefmat.num)

-      c <- solve(e$vectors, odeini)

-      f.out <- function(t) {

-        e$vectors %*% diag(exp(e$values * t), nrow=length(mod_vars)) %*% c

-      }

-      o <- matrix(mapply(f.out, outtimes), 

-        nrow = length(mod_vars), ncol = length(outtimes))

-      dimnames(o) <- list(mod_vars, outtimes)

-      out <- cbind(time = outtimes, t(o))

-    } 

-    if (solution == "deSolve")  

-    {

-      out <- ode(

-        y = odeini,

-        times = outtimes,

-        func = mkindiff, 

-        parms = odeparms,

-        atol = atol

-      )

-    }

-  

-    # Output transformation for models with unobserved compartments like SFORB

-    out_transformed <- data.frame(time = out[,"time"])

-    for (var in names(mkinmod$map)) {

-      if((length(mkinmod$map[[var]]) == 1) || solution == "analytical") {

-        out_transformed[var] <- out[, var]

-      } else {

-        out_transformed[var] <- rowSums(out[, mkinmod$map[[var]]])

-      }

-    }    

-    assign("out_predicted", out_transformed, inherits=TRUE)

+    # Solve the system with current transformed parameter values

+    out <- predict(mkinmod, transparms, odeini, outtimes, solution_type = solution)

+

+    assign("out_predicted", out, inherits=TRUE)

 

-    mC <- modCost(out_transformed, observed, y = "value",

+    mC <- modCost(out, observed, y = "value",

       err = err, weight = weight, scaleVar = scaleVar)

 

     # Report and/or plot if the model is improved

@@ -199,53 +122,9 @@ mkinfit <- function(mkinmod, observed,
       # Plot the data and current model output if requested

       if(plot) {

         outtimes_plot = seq(min(observed$time), max(observed$time), length.out=100)

-        if (solution == "analytical") {

-          o_plot <- switch(parent.type,

-            SFO = SFO.solution(outtimes_plot, 

-                evalparse(parent.name),

-                evalparse(paste("k", parent.name, "sink", sep="_"))),

-            FOMC = FOMC.solution(outtimes_plot,

-                evalparse(parent.name),

-                evalparse("alpha"), evalparse("beta")),

-            DFOP = DFOP.solution(outtimes_plot,

-                evalparse(parent.name),

-                evalparse("k1"), evalparse("k2"),

-                evalparse("g")),

-            HS = HS.solution(outtimes_plot,

-                evalparse(parent.name),

-                evalparse("k1"), evalparse("k2"),

-                evalparse("tb")),

-            SFORB = SFORB.solution(outtimes_plot,

-                evalparse(parent.name),

-                evalparse(paste("k", parent.name, "bound", sep="_")),

-                evalparse(paste("k", sub("free", "bound", parent.name), "free", sep="_")),

-                evalparse(paste("k", parent.name, "sink", sep="_")))

-          )

-          out_plot <- cbind(outtimes_plot, o_plot)

-          dimnames(out_plot) <- list(outtimes_plot, c("time", sub("_free", "", parent.name)))

-        }

-        if(solution == "eigen") {

-          o_plot <- matrix(mapply(f.out, outtimes_plot), 

-            nrow = length(mod_vars), ncol = length(outtimes_plot))

-          dimnames(o_plot) <- list(mod_vars, outtimes_plot)

-          out_plot <- cbind(time = outtimes_plot, t(o_plot))

-        } 

-        if (solution == "deSolve") {

-          out_plot <- ode(

-            y = odeini,

-            times = outtimes_plot,

-            func = mkindiff, 

-            parms = odeparms)

-        }

-        out_transformed_plot <- data.frame(time = out_plot[,"time"])

-        for (var in obs_vars) {

-          if((length(mkinmod$map[[var]]) == 1) || solution == "analytical") {

-            out_transformed_plot[var] <- out_plot[, var]

-          } else {

-            out_transformed_plot[var] <- rowSums(out_plot[, mkinmod$map[[var]]])

-          }

-        }    

-        out_transformed_plot <<- out_transformed_plot

+

+        out_plot <- predict(mkinmod, transparms, odeini, outtimes_plot, 

+          solution_type = solution)

 

         plot(0, type="n", 

           xlim = range(observed$time), ylim = range(observed$value, na.rm=TRUE),

@@ -256,7 +135,7 @@ mkinfit <- function(mkinmod, observed,
           points(subset(observed, name == obs_var, c(time, value)), 

             pch = pch_obs[obs_var], col = col_obs[obs_var])

         }

-        matlines(out_transformed_plot$time, out_transformed_plot[-1])

+        matlines(out_plot$time, out_plot[-1])

         legend("topright", inset=c(0.05, 0.05), legend=obs_vars, 

           col=col_obs, pch=pch_obs, lty=1:length(pch_obs))

       }

@@ -272,12 +151,6 @@ mkinfit <- function(mkinmod, observed,
   if (solution == "eigen") {

     fit$coefmat <- mkinmod$coefmat

   } 

-  if (solution == "deSolve") {

-    fit$mkindiff <- mkindiff

-  }

-  if (plot == TRUE) {

-    fit$out_transformed_plot = out_transformed_plot

-  }

 

   # We also need various other information for summary and plotting

   fit$map <- mkinmod$map

@@ -289,8 +162,6 @@ mkinfit <- function(mkinmod, observed,
   # Collect initial parameter values in two dataframes

   fit$start <- data.frame(initial = c(state.ini.optim, parms.optim))

   fit$start$type = c(rep("state", length(state.ini.optim)), rep("deparm", length(parms.optim)))

-  fit$start$lower <- lower

-  fit$start$upper <- upper

 

   fit$fixed <- data.frame(

     value = c(state.ini.fixed, parms.fixed))

@@ -322,11 +193,10 @@ mkinfit <- function(mkinmod, observed,
   }

   fit$errmin <- errmin

 

-  # Calculate dissipation times DT50 and DT90 and formation fractions

-  parms.all = c(fit$par, parms.fixed)

+  # Calculate dissipation times DT50 and DT90 from parameters

+  parms.all = transform_odeparms(c(fit$par, parms.fixed), mod_vars)

   fit$distimes <- data.frame(DT50 = rep(NA, length(obs_vars)), DT90 = rep(NA, length(obs_vars)), 

     row.names = obs_vars)

-  fit$ff <- vector()

   fit$SFORB <- vector()

   for (obs_var in obs_vars) {

     type = names(mkinmod$map[[obs_var]])[1]  

@@ -374,16 +244,6 @@ mkinfit <- function(mkinmod, observed,
       DT50 <- DTx(50)

       DT90 <- DTx(90)

     }

-    # Back-calculation of formation fractions in case of nonlinear parent kinetics

-    if (type %in% c("FOMC", "DFOP", "HS")) {

-      ff_names = names(mkinmod$ff)

-      for (ff_name in ff_names)

-      {

-        fit$ff[[paste(obs_var, ff_name, sep="_")]] = 

-          eval(parse(text = mkinmod$ff[ff_name]), as.list(parms.all))

-      }

-      fit$ff[[paste(obs_var, "sink", sep="_")]] = 1 - sum(fit$ff)

-    }

     if (type == "SFORB") {

       # FOCUS kinetics (2006), p. 60 f

       k_out_names = grep(paste("k", obs_var, "free", sep="_"), names(parms.all), value=TRUE)

@@ -474,9 +334,10 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, ff = TRUE, ...
   ans$diffs <- object$diffs

   if(data) ans$data <- object$data

   ans$start <- object$start

-  if (object$logk) {

-    names(ans$start) <- sub("^k_", "log k_", names(ans$start))

-  }

+

+  # Fix names of parameters to show that they were transformed

+  names(ans$start) <- sub("^k_", "log k_", names(ans$start))

+

   ans$fixed <- object$fixed

   ans$errmin <- object$errmin 

   if(distimes) ans$distimes <- object$distimes

diff --git a/R/predict.mkinmod.R b/R/predict.mkinmod.R
new file mode 100644
index 0000000..279fa52
--- /dev/null
+++ b/R/predict.mkinmod.R
@@ -0,0 +1,89 @@
+predict.mkinmod <- function(mkinmod, odeparms, odeini, outtimes, solution_type = "deSolve", map_output = TRUE, atol = 1e-6) {
+
+  # Get the names of the state variables in the model
+  mod_vars <- names(mkinmod$diffs)
+
+  # Create function for evaluation of expressions with ode parameters and initial values
+  evalparse <- function(string)
+  {
+    eval(parse(text=string), as.list(c(odeparms, odeini)))
+  }
+
+  # Create a function calculating the differentials specified by the model
+  # if necessary
+  if (solution_type == "analytical") {
+    parent.type = names(mkinmod$map[[1]])[1]  
+    parent.name = names(mkinmod$diffs)[[1]]
+    o <- switch(parent.type,
+      SFO = SFO.solution(outtimes, 
+          evalparse(parent.name),
+          evalparse(paste("k", parent.name, sep="_"))),
+      FOMC = FOMC.solution(outtimes,
+          evalparse(parent.name),
+          evalparse("alpha"), evalparse("beta")),
+      DFOP = DFOP.solution(outtimes,
+          evalparse(parent.name),
+          evalparse("k1"), evalparse("k2"),
+          evalparse("g")),
+      HS = HS.solution(outtimes,
+          evalparse(parent.name),
+          evalparse("k1"), evalparse("k2"),
+          evalparse("tb")),
+      SFORB = SFORB.solution(outtimes,
+          evalparse(parent.name),
+          evalparse(paste("k", parent.name, "bound", sep="_")),
+          evalparse(paste("k", sub("free", "bound", parent.name), "free", sep="_")),
+          evalparse(paste("k", parent.name, sep="_")))
+    )
+    out <- cbind(outtimes, o)
+    dimnames(out) <- list(outtimes, c("time", sub("_free", "", parent.name)))
+  }
+  if (solution_type == "eigen") {
+    coefmat.num <- matrix(sapply(as.vector(mkinmod$coefmat), evalparse), 
+      nrow = length(mod_vars))
+    e <- eigen(coefmat.num)
+    c <- solve(e$vectors, odeini)
+    f.out <- function(t) {
+      e$vectors %*% diag(exp(e$values * t), nrow=length(mod_vars)) %*% c
+    }
+    o <- matrix(mapply(f.out, outtimes), 
+      nrow = length(mod_vars), ncol = length(outtimes))
+    dimnames(o) <- list(mod_vars, outtimes)
+    out <- cbind(time = outtimes, t(o))
+  } 
+  if (solution_type == "deSolve") {
+    mkindiff <- function(t, state, parms) {
+
+      time <- t
+      diffs <- vector()
+      for (box in names(mkinmod$diffs))
+      {
+        diffname <- paste("d", box, sep="_")      
+        diffs[diffname] <- with(as.list(c(time, state, parms)),
+          eval(parse(text=mkinmod$diffs[[box]])))
+      }
+      return(list(c(diffs)))
+    } 
+    out <- ode(
+      y = odeini,
+      times = outtimes,
+      func = mkindiff, 
+      parms = odeparms,
+      atol = atol
+    )
+  }
+  if (map_output) {
+    # Output transformation for models with unobserved compartments like SFORB
+    out_mapped <- data.frame(time = out[,"time"])
+    for (var in names(mkinmod$map)) {
+      if((length(mkinmod$map[[var]]) == 1) || solution == "analytical") {
+        out_mapped[var] <- out[, var]
+      } else {
+        out_mapped[var] <- rowSums(out[, mkinmod$map[[var]]])
+      }
+    }
+    return(out_mapped) 
+  } else {
+    return(out)
+  }
+}
diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R
new file mode 100644
index 0000000..4b8bfd1
--- /dev/null
+++ b/R/transform_odeparms.R
@@ -0,0 +1,46 @@
+transform_odeparms <- function(odeparms, mod_vars) {

+  # Transform rate constants and formation fractions

+  transparms <- odeparms

+  # Exponential transformation for rate constants

+  index_k <- grep("^k_", names(odeparms))

+  if (length(index_k) > 0) {

+    transparms[index_k] <- exp(odeparms[index_k])

+  }

+

+  # Go through state variables and apply inverse isotropic logratio transformation

+  for (box in mod_vars) {

+    indices_f <- grep(paste("^f", box, sep = "_"), names(odeparms))

+    f_names <- grep(paste("^f", box, sep = "_"), names(odeparms), value = TRUE)

+    n_paths <- length(indices_f)

+    if (n_paths > 0) {

+      f <- invilr(odeparms[indices_f])[1:n_paths] # We do not need the last component

+      names(f) <- f_names

+      transparms[indices_f] <- f

+    }

+  }

+  return(transparms)

+}

+

+backtransform_odeparms <- function(transparms, mod_vars) {

+  # Transform rate constants and formation fractions

+  odeparms <- transparms

+  # Log transformation for rate constants

+  index_k <- grep("^k_", names(transparms))

+  if (length(index_k) > 0) {

+    odeparms[index_k] <- log(transparms[index_k])

+  }

+

+  # Go through state variables and apply isotropic logratio transformation

+  for (box in mod_vars) {

+    indices_f <- grep(paste("^f", box, sep = "_"), names(transparms))

+    f_names <- grep(paste("^f", box, sep = "_"), names(transparms), value = TRUE)

+    n_paths <- length(indices_f)

+    if (n_paths > 0) {

+      trans_f <- transparms[indices_f]

+      f <- ilr(c(trans_f, 1 - sum(trans_f)))

+      names(f) <- f_names

+      odeparms[indices_f] <- f

+    }

+  }

+  return(odeparms)

+}

diff --git a/TODO b/TODO
index 3da3339..8f89002 100644
--- a/TODO
+++ b/TODO
@@ -1,5 +1,6 @@
 - Fix analytical and Eivenvalue based solutions after transition to fitting
   transformed k and f values
+- Fix FOMC model with FOCUS_2006_C
 - Transfer calculation of DT50 and DT90 calculations from mkinfit to summary.mkinfit
 - Calculate back-transformed parameters in summary.mkinfit
 - Report back-transformed parameters in print.summary.mkinfit
diff --git a/man/mkinfit.Rd b/man/mkinfit.Rd
index 24dcb83..17e0a11 100644
--- a/man/mkinfit.Rd
+++ b/man/mkinfit.Rd
@@ -34,11 +34,9 @@ mkinfit(mkinmod, observed,
     \code{\link{modFit}}.
   }
   \item{parms.ini}{
-    A named vector if initial values for the parameters, including parameters to
+    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}.
-    The default is to set the initial values to 0.1. The setting of the initial values for 
-    the parameters has a strong impart on performance and it lies in the responsibilty of the 
-    user to set sensible initial values.   
+    If set to "auto", initial values for rate constants are set to default values.
   }
   \item{state.ini}{
     A named vector of initial values for the state variables of the model. In case the