Start of the transition to fitting transformed parameters.

Many things are broken now (see TODO list)



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

4 files changed, 158 insertions, 173 deletions

diff --git a/R/ilr.R b/R/ilr.R
new file mode 100644
index 0000000..f2a6f1b
--- /dev/null
+++ b/R/ilr.R
@@ -0,0 +1,44 @@
+# $Id: ilr.R 120 2011-09-02 14:25:35Z jranke $
+
+# Copyright (C) 2012 René Lehmann, Johannes Ranke
+# Contact: jranke@uni-bremen.de
+
+# This file is part of the R package mkin
+
+# mkin is free software: you can redistribute it and/or modify it under the
+# terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
+
+# This program is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
+# details.
+
+# You should have received a copy of the GNU General Public License along with
+# this program. If not, see <http://www.gnu.org/licenses/>
+
+ilr <- function(x) {
+  z <- vector()
+  for (i in 1:(length(x) - 1)) {
+    z[i] <- sqrt(i/(i+1)) * log((prod(x[1:i]))^(1/i) / x[i+1])
+  }
+  return(z)
+}
+
+invilr<-function(x) {
+  D <- length(x) + 1
+  z <- c(x, 0)
+  y <- rep(0, D)
+  s <- sqrt(1:D*2:(D+1))
+  q <- z/s
+  y[1] <- sum(q[1:D])
+  for (i in 2:D) {
+    y[i] <- sum(q[i:D]) - sqrt((i-1)/i) * z[i-1]
+  }
+  z <- vector()
+  for (i in 1:D) {
+    z[i] <- exp(y[i])/sum(exp(y))
+  }
+  return(z)
+}
diff --git a/R/mkinfit.R b/R/mkinfit.R
index b742742..83896bc 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -1,7 +1,7 @@
 # $Id: mkinfit.R 120 2011-09-02 14:25:35Z jranke $

 

-# Copyright (C) 2010-2011 Johannes Ranke

-# Contact: mkin-devel@lists.berlios.de

+# Copyright (C) 2010-2012 Johannes Ranke

+# Contact: jranke@uni-bremen.de

 # The summary function is an adapted and extended version of summary.modFit

 # from the FME package, v 1.1 by Soetart and Petzoldt, which was in turn

 # inspired by summary.nls.lm

@@ -22,9 +22,8 @@
 # this program. If not, see <http://www.gnu.org/licenses/>

 

 mkinfit <- function(mkinmod, observed,

-  parms.ini = rep(0.1, length(mkinmod$parms)),

+  parms.ini = "auto",

   state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)), 

-  lower = 0, upper = Inf,

   fixed_parms = NULL,

   fixed_initials = names(mkinmod$diffs)[-1],

   eigen = FALSE,

@@ -34,14 +33,24 @@ mkinfit <- function(mkinmod, observed,
   ...)

 {

   mod_vars <- names(mkinmod$diffs)

-  # Subset dataframe with mapped (modelled) variables

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

   observed <- subset(observed, name %in% names(mkinmod$map))

   # Get names of observed variables

   obs_vars = unique(as.character(observed$name))

 

-  # Name the parameters if they are not named yet

-  if(is.null(names(parms.ini))) names(parms.ini) <- mkinmod$parms

-

+  # Define initial 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) {

+    # Default values for the FOMC model

+    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

+  }

+   

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

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

 

@@ -58,14 +67,6 @@ mkinfit <- function(mkinmod, observed,
       names(state.ini.optim) <- paste(names(state.ini.optim), "0", sep="_")

   }

 

-  # Set upper limit for formation fractions to one if formation fractions are

-  # directly defined and if no user input for upper limit is given

-  if (all(upper==Inf) & any(grepl("f_", names(parms.ini)))==TRUE){

-    upper=c( rep(Inf,length(parms.optim)))

-    upper[grep("f_", names(parms.optim))]=1

-    upper=c(rep(Inf, length(state.ini.optim)), upper)

-  }

-

   # 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

@@ -80,12 +81,25 @@ mkinfit <- function(mkinmod, observed,
   # 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)),

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

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

       }

       return(list(c(diffs)))

@@ -251,7 +265,7 @@ mkinfit <- function(mkinmod, observed,
     }

     return(mC)

   }

-  fit <- modFit(cost, c(state.ini.optim, parms.optim), lower = lower, upper = upper, ...)

+  fit <- modFit(cost, c(state.ini.optim, parms.optim), ...)

 

   # We need to return some more data for summary and plotting

   fit$solution <- solution

@@ -323,7 +337,7 @@ mkinfit <- function(mkinmod, observed,
       DT90 = log(10)/k_tot

       for (k_name in k_names)

       {

-        fit$ff[[sub("k_", "", k_name)]] = parms.all[[k_name]] / k_tot

+        fit$ff[[sub("^k_", "", k_name)]] = parms.all[[k_name]] / k_tot

       }

     }

     if (type == "FOMC") {

@@ -395,7 +409,7 @@ mkinfit <- function(mkinmod, observed,
       if (abs(DT90.o - max_DT) < 0.01) DT90 = NA else DT90 = DT90.o

       for (k_out_name in k_out_names)

       {

-        fit$ff[[sub("k_", "", k_out_name)]] = parms.all[[k_out_name]] / k_1output

+        fit$ff[[sub("^k_", "", k_out_name)]] = parms.all[[k_out_name]] / k_1output

       }

       # Return the eigenvalues for comparison with DFOP rate constants

       fit$SFORB[[paste(obs_var, "b1", sep="_")]] = b1

@@ -416,8 +430,11 @@ mkinfit <- function(mkinmod, observed,
   return(fit)

 }

 

-summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, ff = TRUE, cov = FALSE,...) {

+summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, ff = TRUE, ...) {

   param  <- object$par

+  if (object$logk) {

+    names(param) <- sub("^k_", "log k_", names(param))

+  }

   pnames <- names(param)

   p      <- length(param)

   covar  <- try(solve(0.5*object$hessian), silent = TRUE)   # unscaled covariance

@@ -441,32 +458,25 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, ff = TRUE, cov
     warning(message)

   } else message <- "ok"

 

-  param <- cbind(param, se, tval, pt(tval, rdf, lower.tail = FALSE))

-  dimnames(param) <- list(pnames, c("Estimate", "Std. Error",

-                                    "t value", "Pr(>t)"))

-  if(cov)

-    ans <- list(residuals = object$residuals,

-                residualVariance = resvar,

-                sigma = sqrt(resvar),

-                modVariance = modVariance,

-                df = c(p, rdf), cov.unscaled = covar,

-                cov.scaled = covar * resvar,

-                info = object$info, niter = object$iterations,

-                stopmess = message,

-                par = param)

-  else

-    ans <- list(residuals = object$residuals,

-                residualVariance = resvar,

-                sigma = sqrt(resvar),

-                modVariance = modVariance,

-                df = c(p, rdf),

-                info = object$info, niter = object$iterations,

-                stopmess = message,

-                par = param)

+  param <- cbind(param, se)

+  dimnames(param) <- list(pnames, c("Estimate", "Std. Error"))

+

+  ans <- list(residuals = object$residuals,

+          residualVariance = resvar,

+          sigma = sqrt(resvar),

+          modVariance = modVariance,

+          df = c(p, rdf), cov.unscaled = covar,

+          cov.scaled = covar * resvar,

+          info = object$info, niter = object$iterations,

+          stopmess = message,

+          par = param)

 

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

+  }

   ans$fixed <- object$fixed

   ans$errmin <- object$errmin 

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

diff --git a/R/mkinmod.R b/R/mkinmod.R
index 31c778c..b49c813 100644
--- a/R/mkinmod.R
+++ b/R/mkinmod.R
@@ -1,7 +1,7 @@
 # $Id: mkinmod.R 71 2010-09-12 01:13:36Z jranke $

 

-# Copyright (C) 2010 Johannes Ranke

-# Contact: mkin-devel@lists.berlios.de

+# Copyright (C) 2010-2012 Johannes Ranke

+# Contact: jranke@uni-bremen.de

 

 # This file is part of the R package mkin

 

@@ -33,21 +33,27 @@ mkinmod <- function(...)
   if(spec[[1]]$type %in% c("FOMC", "DFOP", "HS")) {

     mat = FALSE 

     if(!is.null(spec[[1]]$to)) {

-      message <- paste("Only constant formation fractions over time are implemented.",

-        "Depending on the reason for the time dependence of degradation, this may be unrealistic",

+      message <- paste(

+        "Only constant formation fractions over time are implemented.",

+        "Depending on the reason for the time dependence of degradation",

+        "this may be unrealistic",

         sep="\n")

       warning(message)

     } else message <- "ok"

   } else mat = TRUE

 

   # Establish list of differential equations

+  # as well as map from observed compartments to differential equations

   for (varname in obs_vars)

   {

     if(is.null(spec[[varname]]$type)) stop(

       "Every argument to mkinmod must be a list containing a type component")

     if(!spec[[varname]]$type %in% c("SFO", "FOMC", "DFOP", "HS", "SFORB")) stop(

       "Available types are SFO, FOMC, DFOP, HS and SFORB only")

-    new_parms <- vector()

+    if(spec[[varname]]$type %in% c("FOMC", "DFOP", "HS") & match(varname, obs_vars) != 1) {

+        stop(paste("Types FOMC, DFOP and HS are only implemented for the first compartment,", 

+                   "which is assumed to be the source compartment"))

+    }

 

     # New (sub)compartments (boxes) needed for the model type

     new_boxes <- switch(spec[[varname]]$type,

@@ -62,124 +68,75 @@ mkinmod <- function(...)
 

     # Start a new differential equation for each new box

     new_diffs <- paste("d_", new_boxes, " =", sep="")

+    names(new_diffs) <- new_boxes

+    diffs <- c(diffs, new_diffs)

+  }

+

+  # Create content of differential equations and build parameter list

+  for (varname in obs_vars)

+  {

+    # Add first-order term to first (or only) box for SFO and SFORB

+    box_1 = map[[varname]][[1]] # This is the only box unless type is SFORB

+    if(spec[[varname]]$type %in% c("SFO", "SFORB")) {

+      k_compound <- paste("k", box_1, sep="_")

+      origin_term <- paste(k_compound, "*", box_1)

+      parms <- c(parms, k_compound)

+    }

 

-    # Turn on sink if not specified otherwise

-    if(is.null(spec[[varname]]$sink)) spec[[varname]]$sink <- TRUE

 

     # Construct and add FOMC term and add FOMC parameters if needed

     if(spec[[varname]]$type == "FOMC") {

-      if(match(varname, obs_vars) != 1) {

-        stop("Type FOMC is only possible for the first compartment, which is assumed to be the source compartment")

-      }

-      if(spec[[varname]]$sink == FALSE) {

-        stop("Turning off the sink for the FOMC model is not implemented")

-      }

       # From p. 53 of the FOCUS kinetics report

-      nonlinear_term <- paste("(alpha/beta) * ((time/beta) + 1)^-1 *", new_boxes[[1]])

-      new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term)

-      new_parms <- c("alpha", "beta")

-      ff <- vector()

+      origin_term <- paste("(alpha/beta) * ((time/beta) + 1)^-1 *", box_1)

+      parms <- c(parms, "alpha", "beta")

     } 

 

     # Construct and add DFOP term and add DFOP parameters if needed

     if(spec[[varname]]$type == "DFOP") {

-      if(match(varname, obs_vars) != 1) {

-        stop("Type DFOP is only possible for the first compartment, which is assumed to be the source compartment")

-      }

-      if(spec[[varname]]$sink == FALSE) {

-        stop("Turning off the sink for the DFOP model is not implemented")

-      }

       # From p. 57 of the FOCUS kinetics report

-      nonlinear_term <- paste("((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) *", new_boxes[[1]])

-      new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term)

-      new_parms <- c("k1", "k2", "g")

-      ff <- vector()

+      origin_term <- paste("((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) *", box_1)

+      parms <- c(parms, "k1", "k2", "g")

     } 

 

     # Construct and add HS term and add HS parameters if needed

     if(spec[[varname]]$type == "HS") {

-      if(match(varname, obs_vars) != 1) {

-        stop("Type HS is only possible for the first compartment, which is assumed to be the source compartment")

-      }

-      if(spec[[varname]]$sink == FALSE) {

-        stop("Turning off the sink for the HS model is not implemented")

-      }

       # From p. 55 of the FOCUS kinetics report

-      nonlinear_term <- paste("ifelse(time <= tb, k1, k2)", "*", new_boxes[[1]])

-      new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term)

-      new_parms <- c("k1", "k2", "tb")

-      ff <- vector()

+      origin_term <- paste("ifelse(time <= tb, k1, k2)", "*", box_1)

+      parms <- c(parms, "k1", "k2", "tb")

     } 

 

-    # Construct terms for transfer to sink and add if appropriate

-

-    if(spec[[varname]]$sink) {

-      # Add first-order sink term to first (or only) box for SFO and SFORB

-      if(spec[[varname]]$type %in% c("SFO", "SFORB")) {

-        k_compound_sink <- paste("k", new_boxes[[1]], "sink", sep="_")

-        sink_term <- paste("-", k_compound_sink, "*", new_boxes[[1]])

-        new_diffs[[1]] <- paste(new_diffs[[1]], sink_term)

-        new_parms <- k_compound_sink

-      }

-    }

+    # Add origin decline term to box 1 (usually the only box, unless type is SFORB)

+    diffs[[box_1]] <- paste(diffs[[box_1]], "-", origin_term)

    

-    # Add reversible binding if appropriate

+    # Add reversible binding for SFORB models

     if(spec[[varname]]$type == "SFORB") {

-      k_free_bound <- paste("k", varname, "free", "bound", sep="_")      

-      k_bound_free <- paste("k", varname, "bound", "free", sep="_")      

-      reversible_binding_terms <- c(

-        paste("-", k_free_bound, "*", new_boxes[[1]], "+", k_bound_free, "*", new_boxes[[2]]),

-        paste("+", k_free_bound, "*", new_boxes[[1]], "-", k_bound_free, "*", new_boxes[[2]]))

-      new_diffs <- paste(new_diffs, reversible_binding_terms)

-      new_parms <- c(new_parms, k_free_bound, k_bound_free)

+      box_2 = map[[varname]][[2]]

+      k_free_bound <- paste("k", varname, "free", "bound", sep="_")

+      k_bound_free <- paste("k", varname, "bound", "free", sep="_")

+      reversible_binding_term_1 <- paste("-", k_free_bound, "*", box_1, "+",

+        k_bound_free, "*", box_2)

+      diffs[[box_1]] <- paste(diffs[[box_1]], reversible_binding_term_1)

+      reversible_binding_term_2 <- paste("+", k_free_bound, "*", box_1, "-",

+        k_bound_free, "*", box_2)

+      diffs[[box_2]] <- paste(diffs[[box_2]], reversible_binding_term_2)

+      parms <- c(parms, k_free_bound, k_bound_free)

     } 

 

-    # Add observed variable to model

-    parms <- c(parms, new_parms)

-    names(new_diffs) <- new_boxes

-    diffs <- c(diffs, new_diffs)

-  }

-  # Transfer between compartments

-  for (varname in obs_vars) {

+    # Transfer between compartments

     to <- spec[[varname]]$to

     if(!is.null(to)) {

-      origin_box <- switch(spec[[varname]]$type,

-        SFO = varname,

-        FOMC = varname,

-        DFOP = varname,

-        HS = varname,

-        SFORB = paste(varname, "free", sep="_"))

-      fraction_left <- NULL

+      # Name of box from which transfer takes place

+      origin_box <- box_1

+

+      # Add transfer terms to listed compartments

       for (target in to) {

         target_box <- switch(spec[[target]]$type,

           SFO = target,

           SFORB = paste(target, "free", sep="_"))

-        if(spec[[varname]]$type %in% c("SFO", "SFORB")) {

-          k_from_to <- paste("k", origin_box, target_box, sep="_")

-          diffs[[origin_box]] <- paste(diffs[[origin_box]], "-", 

-            k_from_to, "*", origin_box)

-          diffs[[target_box]] <- paste(diffs[[target_box]], "+", 

-            k_from_to, "*", origin_box)

-          parms <- c(parms, k_from_to)

-        }

-        if(spec[[varname]]$type %in% c("FOMC", "DFOP", "HS")) {

-          fraction_to_target = paste("f_to", target, sep="_")

-          fraction_not_to_target = paste("(1 - ", fraction_to_target, ")", 

-            sep="")

-          if(is.null(fraction_left)) {

-            fraction_really_to_target = fraction_to_target

-            fraction_left = fraction_not_to_target

-          } else {

-            fraction_really_to_target = paste(fraction_left, " * ", 

-              fraction_to_target, sep="")

-            fraction_left = paste(fraction_left, " * ", 

-              fraction_not_to_target, sep="")

-          }

-          ff[target_box] = fraction_really_to_target

-          diffs[[target_box]] <- paste(diffs[[target_box]], "+", 

-            ff[target_box], "*", nonlinear_term)

-          parms <- c(parms, fraction_to_target)

-        }

+        fraction_to_target = paste("f", origin_box, "to", target, sep="_")

+        diffs[[target_box]] <- paste(diffs[[target_box]], "+", 

+            fraction_to_target, "*", origin_term)

+        parms <- c(parms, fraction_to_target)

       }

     }

   }

@@ -193,26 +150,23 @@ mkinmod <- function(...)
     for (from in boxes) {

       for (to in boxes) {

         if (from == to) {

-          k.candidate = paste("k", from, c(boxes, "sink"), sep="_")

-          k.candidate = sub("free.*bound", "free_bound", k.candidate)

-          k.candidate = sub("bound.*free", "bound_free", k.candidate)

-          k.effective = intersect(model$parms, k.candidate)

-          m[from,to] = ifelse(length(k.effective) > 0,

-              paste("-", k.effective, collapse = " "), "0")

+          k.candidate = paste("k", from, sep="_")

+          m[from,to] = ifelse(k.candidate %in% model$parms,

+              paste("-", k.candidate), "0")

         } else {

+          f.candidate = paste("f", from, "to", to, sep="_")

           k.candidate = paste("k", from, to, sep="_")

           k.candidate = sub("free.*bound", "free_bound", k.candidate)

           k.candidate = sub("bound.*free", "bound_free", k.candidate)

-          k.effective = intersect(model$parms, k.candidate)

-          m[to, from] = ifelse(length(k.effective) > 0,

-              k.effective, "0")

+          m[to, from] = ifelse(f.candidate %in% model$parms,

+            paste(f.candidate, " * k_", from, sep=""), 

+            ifelse(k.candidate %in% model$parms, k.candidate, "0"))

         }

       }

     }

     model$coefmat <- m

   }

 

-  if (exists("ff")) model$ff = ff

   class(model) <- "mkinmod"

-  invisible(model)

+  return(model)

 }

diff --git a/R/mkinstart.R b/R/mkinstart.R
deleted file mode 100644
index 46e5bed..0000000
--- a/R/mkinstart.R
+++ /dev/null
@@ -1,23 +0,0 @@
-mkinstart <- function(model, data, mode = "auto")

-{

-  if (class(model) != "mkinmod") stop("The first argument must be a model of class mkinmod")

-  names <- model$parms

-  observed <- names(model$map)

-  if(!all(observed %in% levels(data$name))) stop("The data must contain the observed variables used in the model")

-  for (obs in observed)

-  {

-    tmp <- subset(data, name == obs)

-    max <- tmp[which.max(tmp$value), ]

-    type = names(model$map[[obs]])[[1]]

-    kinmodel <- ifelse(type == "SFORB", "DFOP", type)

-    tmp.longdata <- subset(data, name == obs & time >= max$time)

-    tmp.widedata <- mkin_long_to_wide(tmp.longdata, outtime = "t")

-    names(tmp.widedata) <- c("t", "parent")

-    tmp.fit <- kinfit(

-      kindata = tmp.widedata,

-      kinmodels = kinmodel,

-      parent.0.user = max$value)

-    if(class(tmp.fit[[kinmodel]]) == "try-error") stop(paste("Automatic generation of starting parameters failed\nkinfit failed to find a", kinmodel, "fit for", obs))

-    tmp.results <- kinresults(tmp.fit)

-  }

-}