Update kinfit and mkin to the latest version published on BerliOS.

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

13 files changed, 694 insertions, 109 deletions

diff --git a/R/DFOP.solution.R b/R/DFOP.solution.R
new file mode 100644
index 00000000..8531cfed
--- /dev/null
+++ b/R/DFOP.solution.R
@@ -0,0 +1,5 @@
+DFOP.solution <- function(t, parent.0, k1, k2, g)

+{

+	parent = g * parent.0 * exp(-k1 * t) +

+		 (1 - g) * parent.0 * exp(-k2 * t)

+}

diff --git a/R/FOMC.solution.R b/R/FOMC.solution.R
new file mode 100644
index 00000000..8bd13d6b
--- /dev/null
+++ b/R/FOMC.solution.R
@@ -0,0 +1,4 @@
+FOMC.solution <- function(t, parent.0, alpha, beta)

+{

+	parent = parent.0 / (t/beta + 1)^alpha

+}

diff --git a/R/HS.solution.R b/R/HS.solution.R
new file mode 100644
index 00000000..4651a6a8
--- /dev/null
+++ b/R/HS.solution.R
@@ -0,0 +1,6 @@
+HS.solution <- function(t, parent.0, k1, k2, tb)

+{

+	parent = ifelse(t <= tb, 

+		parent.0 * exp(-k1 * t),

+		parent.0 * exp(-k1 * tb) * exp(-k2 * (t - tb)))

+}

diff --git a/R/SFO.solution.R b/R/SFO.solution.R
new file mode 100644
index 00000000..3a376e48
--- /dev/null
+++ b/R/SFO.solution.R
@@ -0,0 +1,4 @@
+SFO.solution <- function(t, parent.0, k)

+{

+	parent = parent.0 * exp(-k * t)

+}

diff --git a/R/SFORB.solution.R b/R/SFORB.solution.R
new file mode 100644
index 00000000..45a4533a
--- /dev/null
+++ b/R/SFORB.solution.R
@@ -0,0 +1,9 @@
+SFORB.solution = function(t, parent.0, k_12, k_21, k_1output) {

+  sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 + k_12 * k_21 - (k_12 + k_1output) * k_21)

+  b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp

+  b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp

+  

+  parent = parent.0 * 

+        (((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) +

+        ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t))

+}

diff --git a/R/mkin_long_to_wide.R b/R/mkin_long_to_wide.R
index 45ee0a10..ca09b495 100644
--- a/R/mkin_long_to_wide.R
+++ b/R/mkin_long_to_wide.R
@@ -1,7 +1,28 @@
-mkin_long_to_wide <- function(long_data, time = "time")

+# $Id: mkin_long_to_wide.R 96 2011-04-29 11:10:40Z jranke $

+

+# Copyright (C) 2010-2011 Johannes Ranke

+# Contact: mkin-devel@lists.berlios.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/>

+

+mkin_long_to_wide <- function(long_data, time = "time", outtime = "time")

 {

   colnames <- unique(long_data$name)

   wide_data <- data.frame(time = subset(long_data, name == colnames[1], time))

+  names(wide_data) <- outtime

   for (var in colnames) {

     wide_data[var] <- subset(long_data, name == var, value)

   }

diff --git a/R/mkin_wide_to_long.R b/R/mkin_wide_to_long.R
index 21ab77b9..995a9be5 100644
--- a/R/mkin_wide_to_long.R
+++ b/R/mkin_wide_to_long.R
@@ -1,3 +1,23 @@
+# $Id: mkin_wide_to_long.R 59 2010-07-28 12:29:15Z jranke $

+

+# Copyright (C) 2010 Johannes Ranke

+# Contact: mkin-devel@lists.berlios.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/>

+

 mkin_wide_to_long <- function(wide_data, time = "t")

 {

   colnames <- names(wide_data)

diff --git a/R/mkinerrmin.R b/R/mkinerrmin.R
index 7a922eb4..ebed191a 100644
--- a/R/mkinerrmin.R
+++ b/R/mkinerrmin.R
@@ -1,3 +1,23 @@
+# $Id: mkinerrmin.R 59 2010-07-28 12:29:15Z jranke $

+

+# Copyright (C) 2010 Johannes Ranke

+# Contact: mkin-devel@lists.berlios.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/>

+

 mkinerrmin <- function(errdata, n.parms, alpha = 0.05)

 {

   means.mean <- mean(errdata$value_mean, na.rm=TRUE)

diff --git a/R/mkinfit.R b/R/mkinfit.R
index f8285fcd..b7427421 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -1,11 +1,36 @@
+# $Id: mkinfit.R 120 2011-09-02 14:25:35Z jranke $

+

+# Copyright (C) 2010-2011 Johannes Ranke

+# Contact: mkin-devel@lists.berlios.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

+

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

+

 mkinfit <- function(mkinmod, observed,

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

   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,

   plot = FALSE, quiet = FALSE,

   err = NULL, weight = "none", scaleVar = FALSE,

+  atol = 1e-6,

   ...)

 {

   mod_vars <- names(mkinmod$diffs)

@@ -16,18 +41,6 @@ mkinfit <- function(mkinmod, observed,
 

   # Name the parameters if they are not named yet

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

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

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

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

-  } 

 

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

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

@@ -45,6 +58,40 @@ 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

+  if (length(mkinmod$map) == 1) {

+    solution = "analytical"

+  } else {

+    if (is.matrix(mkinmod$coefmat) & eigen) solution = "eigen"

+    else solution = "deSolve"

+  }

+

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

+  # if necessary

+  if(solution == "deSolve") {

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

+      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

   out_predicted <- NA

@@ -60,19 +107,67 @@ 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)))

+    }

 

-    # Solve the ode

-    out <- ode(

-      y = odeini,

-      times = outtimes,

-      func = mkindiff, 

-      parms = odeparms)

-

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

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

         out_transformed[var] <- out[, var]

       } else {

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

@@ -90,19 +185,53 @@ 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)

-        out_plot <- ode(

-          y = odeini,

-          times = outtimes_plot,

-          func = mkindiff, 

-          parms = odeparms)

+        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 names(mkinmod$map)) {

-          if(length(mkinmod$map[[var]]) == 1) {

+        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

 

         plot(0, type="n", 

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

@@ -124,8 +253,17 @@ mkinfit <- function(mkinmod, observed,
   }

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

 

-  # We need the function for plotting

-  fit$mkindiff <- mkindiff

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

+  fit$solution <- solution

+  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

@@ -161,54 +299,109 @@ mkinfit <- function(mkinmod, observed,
     n.optim <- n.k.optim + n.initials.optim

     if ("alpha" %in% names(parms.optim)) n.optim <- n.optim + 1

     if ("beta" %in% names(parms.optim)) n.optim <- n.optim + 1

+    if ("k1" %in% names(parms.optim)) n.optim <- n.optim + 1

+    if ("k2" %in% names(parms.optim)) n.optim <- n.optim + 1

+    if ("g" %in% names(parms.optim)) n.optim <- n.optim + 1

+    if ("tb" %in% names(parms.optim)) n.optim <- n.optim + 1

     errmin.tmp <- mkinerrmin(errdata.var, n.optim)

     errmin[obs_var, c("err.min", "n.optim", "df")] <- errmin.tmp

   }

   fit$errmin <- errmin

 

-  # Calculate dissipation times DT50 and DT90

+  # Calculate dissipation times DT50 and DT90 and formation fractions

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

   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]  

-      if (type == "SFO") {

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

-        k_tot = sum(parms.all[k_names])

-        DT50 = log(2)/k_tot

-        DT90 = log(10)/k_tot

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

+    if (type == "SFO") {

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

+      k_tot = sum(parms.all[k_names])

+      DT50 = log(2)/k_tot

+      DT90 = log(10)/k_tot

+      for (k_name in k_names)

+      {

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

       }

-      if (type == "FOMC") {

-        alpha = parms.all["alpha"]

-        beta = parms.all["beta"]

-        DT50 = beta * (2^(1/alpha) - 1)

-        DT90 = beta * (10^(1/alpha) - 1)

+    }

+    if (type == "FOMC") {

+      alpha = parms.all["alpha"]

+      beta = parms.all["beta"]

+      DT50 = beta * (2^(1/alpha) - 1)

+      DT90 = beta * (10^(1/alpha) - 1)

+    }

+    if (type == "DFOP") {

+      k1 = parms.all["k1"]

+      k2 = parms.all["k2"]

+      g = parms.all["g"]

+      f <- function(t, x) {

+        fraction <- g * exp( - k1 * t) + (1 - g) * exp( - k2 * t)

+        (fraction - (1 - x/100))^2

       }

-      if (type == "SFORB") {

-        # FOCUS kinetics (2006), p. 60 f

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

-        k_out_names = setdiff(k_out_names, paste("k", obs_var, "free", "bound", sep="_"))

-        k_1output = sum(parms.all[k_out_names])

-        k_12 = parms.all[paste("k", obs_var, "free", "bound", sep="_")]

-        k_21 = parms.all[paste("k", obs_var, "bound", "free", sep="_")]

-

-        sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 + k_12 * k_21 - (k_12 + k_1output) * k_21)

-        b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp

-        b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp

-

-        SFORB_fraction = function(t) {

-          ((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) +

-          ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t)

-        }

-        f_50 <- function(t) (SFORB_fraction(t) - 0.5)^2

-        max_DT <- 1000

-        DT50.o <- optimize(f_50, c(0.01, max_DT))$minimum

-        if (abs(DT50.o - max_DT) < 0.01) DT50 = NA else DT50 = DT50.o

-        f_90 <- function(t) (SFORB_fraction(t) - 0.1)^2

-        DT90.o <- optimize(f_90, c(0.01, 1000))$minimum

-        if (abs(DT90.o - max_DT) < 0.01) DT90 = NA else DT90 = DT90.o

+      DTmax <- 1000

+      DT50.o <- optimize(f, c(0.001, DTmax), x=50)$minimum

+      DT50 = ifelse(DTmax - DT50.o < 0.1, NA, DT50.o)

+      DT90.o <- optimize(f, c(0.001, DTmax), x=90)$minimum

+      DT90 = ifelse(DTmax - DT90.o < 0.1, NA, DT90.o)

+    }

+    if (type == "HS") {

+      k1 = parms.all["k1"]

+      k2 = parms.all["k2"]

+      tb = parms.all["tb"]

+      DTx <- function(x) {

+        DTx.a <- (log(100/(100 - x)))/k1

+        DTx.b <- tb + (log(100/(100 - x)) - k1 * tb)/k2

+        if (DTx.a < tb) DTx <- DTx.a

+        else DTx <- DTx.b

+        return(DTx)

       }

-      fit$distimes[obs_var, ] = c(DT50, DT90)

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

+      k_out_names = setdiff(k_out_names, paste("k", obs_var, "free", "bound", sep="_"))

+      k_1output = sum(parms.all[k_out_names])

+      k_12 = parms.all[paste("k", obs_var, "free", "bound", sep="_")]

+      k_21 = parms.all[paste("k", obs_var, "bound", "free", sep="_")]

+

+      sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 + k_12 * k_21 - (k_12 + k_1output) * k_21)

+      b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp

+      b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp

+

+      SFORB_fraction = function(t) {

+        ((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) +

+        ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t)

+      }

+      f_50 <- function(t) (SFORB_fraction(t) - 0.5)^2

+      max_DT <- 1000

+      DT50.o <- optimize(f_50, c(0.01, max_DT))$minimum

+      if (abs(DT50.o - max_DT) < 0.01) DT50 = NA else DT50 = DT50.o

+      f_90 <- function(t) (SFORB_fraction(t) - 0.1)^2

+      DT90.o <- optimize(f_90, c(0.01, max_DT))$minimum

+      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

+      }

+      # Return the eigenvalues for comparison with DFOP rate constants

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

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

+    }

+    fit$distimes[obs_var, ] = c(DT50, DT90)

   }

 

   # Collect observed, predicted and residuals

@@ -217,25 +410,74 @@ mkinfit <- function(mkinmod, observed,
   data$residual <- data$observed - data$predicted

   data$variable <- ordered(data$variable, levels = obs_vars)

   fit$data <- data[order(data$variable, data$time), ]

+  fit$atol <- atol

 

   class(fit) <- c("mkinfit", "modFit") 

   return(fit)

 }

 

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

-  ans <- FME:::summary.modFit(object, cov = cov)

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

+  param  <- object$par

+  pnames <- names(param)

+  p      <- length(param)

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

+  if (!is.numeric(covar)) {

+    message <- "Cannot estimate covariance; system is singular"

+    warning(message)

+    covar <- matrix(data = NA, nrow = p, ncol = p)

+  } else message <- "ok"

+

+  rownames(covar) <- colnames(covar) <-pnames

+  rdf    <- object$df.residual

+  resvar <- object$ssr / rdf

+  se     <- sqrt(diag(covar) * resvar)

+  names(se) <- pnames

+  tval      <- param / se

+  modVariance <- object$ssr / length(object$residuals)

+

+  if (!all(object$start$lower >=0)) {

+    message <- "Note that the one-sided t-test may not be appropriate if

+      parameter values below zero are possible."

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

+

   ans$diffs <- object$diffs

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

   ans$start <- object$start

   ans$fixed <- object$fixed

   ans$errmin <- object$errmin 

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

+  if(ff) ans$ff <- object$ff

+  if(length(object$SFORB) != 0) ans$SFORB <- object$SFORB

   class(ans) <- c("summary.mkinfit", "summary.modFit") 

   return(ans)  

 }

 

 # Expanded from print.summary.modFit

-print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), ...) {

+print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), tval = TRUE, ...) {

   cat("\nEquations:\n")

   print(noquote(as.character(x[["diffs"]])))

   df  <- x$df

@@ -249,7 +491,11 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), .
   else print(x$fixed)

   

   cat("\nOptimised parameters:\n")

-  printCoefmat(x$par, digits = digits, ...)

+  if (tval) printCoefmat(x$par, digits = digits, ...)

+  else {

+    printCoefmat(x$par[,c(1,2,4)], cs.in = c(1,2), tst.ind = integer(0), 

+      P.values = TRUE, has.Pvalue = TRUE, digits = digits, ...)

+  }

 

   cat("\nResidual standard error:",

       format(signif(x$sigma, digits)), "on", rdf, "degrees of freedom\n")

@@ -260,10 +506,22 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), .
 

   printdistimes <- !is.null(x$distimes)

   if(printdistimes){

-    cat("\nEstimated disappearance times\n")

+    cat("\nEstimated disappearance times:\n")

     print(x$distimes, digits=digits,...)

   }    

 

+  printff <- !is.null(x$ff)

+  if(printff){

+    cat("\nEstimated formation fractions:\n")

+    print(data.frame(ff = x$ff), digits=digits,...)

+  }    

+

+  printSFORB <- !is.null(x$SFORB)

+  if(printSFORB){

+    cat("\nEstimated Eigenvalues of SFORB model(s):\n")

+    print(x$SFORB, digits=digits,...)

+  }    

+

   printcor <- !is.null(x$cov.unscaled)

   if (printcor){

     Corr <- cov2cor(x$cov.unscaled)

diff --git a/R/mkinmod.R b/R/mkinmod.R
index 643e855e..31c778cb 100644
--- a/R/mkinmod.R
+++ b/R/mkinmod.R
@@ -1,3 +1,23 @@
+# $Id: mkinmod.R 71 2010-09-12 01:13:36Z jranke $

+

+# Copyright (C) 2010 Johannes Ranke

+# Contact: mkin-devel@lists.berlios.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/>

+

 mkinmod <- function(...)

 {

   spec <- list(...)

@@ -5,22 +25,36 @@ mkinmod <- function(...)
 

   # The returned model will be a list of character vectors, containing

   # differential equations, parameter names and a mapping from model variables

-  # to observed variables

+  # to observed variables. If possible, a matrix representation of the 

+  # differential equations is included

   parms <- vector()

   diffs <- vector()

   map <- list()

+  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",

+        sep="\n")

+      warning(message)

+    } else message <- "ok"

+  } else mat = TRUE

 

   # Establish list of 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", "SFORB")) stop("Available types are SFO, FOMC and SFORB only")

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

 

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

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

       SFO = varname,

       FOMC = varname,

+      DFOP = varname,

+      HS = varname,

       SFORB = paste(varname, c("free", "bound"), sep="_")

     )

     map[[varname]] <- new_boxes

@@ -29,25 +63,64 @@ mkinmod <- function(...)
     # Start a new differential equation for each new box

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

 

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

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

+    } 

+

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

+    } 

+

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

+    } 

+

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

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

       }

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

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

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

-        }

-        # From p. 53 of the FOCUS kinetics report

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

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

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

-      }

     }

    

     # Add reversible binding if appropriate

@@ -66,7 +139,6 @@ mkinmod <- function(...)
     names(new_diffs) <- new_boxes

     diffs <- c(diffs, new_diffs)

   }

-

   # Transfer between compartments

   for (varname in obs_vars) {

     to <- spec[[varname]]$to

@@ -74,6 +146,8 @@ mkinmod <- function(...)
       origin_box <- switch(spec[[varname]]$type,

         SFO = varname,

         FOMC = varname,

+        DFOP = varname,

+        HS = varname,

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

       fraction_left <- NULL

       for (target in to) {

@@ -82,27 +156,63 @@ mkinmod <- function(...)
           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)

+          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 == "FOMC") {

+        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="")

+          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="")

+            fraction_really_to_target = paste(fraction_left, " * ", 

+              fraction_to_target, sep="")

+            fraction_left = paste(fraction_left, " * ", 

+              fraction_not_to_target, sep="")

           }

-          diffs[[target_box]] <- paste(diffs[[target_box]], "+", fraction_really_to_target, "*", fomc_term)

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

         }

       }

     }

   }

   model <- list(diffs = diffs, parms = parms, map = map)

+

+  # Create coefficient matrix if appropriate

+  if (mat) {

+    boxes <- names(diffs)

+    n <- length(boxes)

+    m <- matrix(nrow=n, ncol=n, dimnames=list(boxes, boxes))

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

+        } else {

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

+        }

+      }

+    }

+    model$coefmat <- m

+  }

+

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

   class(model) <- "mkinmod"

   invisible(model)

 }

diff --git a/R/mkinplot.R b/R/mkinplot.R
index d2880c0e..d665bc23 100644
--- a/R/mkinplot.R
+++ b/R/mkinplot.R
@@ -1,5 +1,6 @@
-mkinplot <- function(fit, xlab = "Time", ylab = "Observed", xlim = range(fit$data$time), ylim = range(fit$data$observed, na.rm = TRUE), ...)
+mkinplot <- function(fit, xlab = "Time", ylab = "Observed", xlim = range(fit$data$time), ylim = range(fit$data$observed, na.rm = TRUE), legend = TRUE, ...)
 {
+  solution = fit$solution
   fixed <- fit$fixed$value
   names(fixed) <- rownames(fit$fixed)
   parms.all <- c(fit$par, fixed)
@@ -16,12 +17,60 @@ mkinplot <- function(fit, xlab = "Time", ylab = "Observed", xlim = range(fit$dat
     rownames(subset(fit$fixed, type == "deparm")))
   odeparms <- parms.all[odenames]
 
-  # Solve the ode
-  out <- ode(
-    y = odeini,
-    times = outtimes,
-    func = fit$mkindiff, 
-    parms = odeparms)
+  # Solve the system
+  evalparse <- function(string)
+  {
+    eval(parse(text=string), as.list(c(odeparms, odeini)))
+  }
+  if (solution == "analytical") {
+    parent.type = names(fit$map[[1]])[1]  
+    parent.name = names(fit$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, "free_bound", sep="_")),
+          evalparse(paste("k", parent.name, "bound_free", sep="_")),
+          evalparse(paste("k", parent.name, "free_sink", sep="_")))
+    )
+    out <- cbind(outtimes, o)
+    dimnames(out) <- list(outtimes, c("time", parent.name))
+  }
+  if (solution == "eigen") {
+    coefmat.num <- matrix(sapply(as.vector(fit$coefmat), evalparse), 
+      nrow = length(odeini))
+    e <- eigen(coefmat.num)
+    c <- solve(e$vectors, odeini)
+    f.out <- function(t) {
+      e$vectors %*% diag(exp(e$values * t), nrow=length(odeini)) %*% c
+    }
+    o <- matrix(mapply(f.out, outtimes), 
+      nrow = length(odeini), ncol = length(outtimes))
+    dimnames(o) <- list(names(odeini), NULL)
+    out <- cbind(time = outtimes, t(o))
+  } 
+  if (solution == "deSolve") {
+    out <- ode(
+      y = odeini,
+      times = outtimes,
+      func = fit$mkindiff, 
+      parms = odeparms,
+      atol = fit$atol
+    )
+  }
     
   # Output transformation for models with unobserved compartments like SFORB
   out_transformed <- data.frame(time = out[,"time"])
@@ -36,14 +85,16 @@ mkinplot <- function(fit, xlab = "Time", ylab = "Observed", xlim = range(fit$dat
   # Plot the data and model output
   plot(0, type="n", 
     xlim = xlim, ylim = ylim,
-      xlab = xlab, ylab = ylab, ...)
-    col_obs <- pch_obs <- 1:length(fit$map)
-    names(col_obs) <- names(pch_obs) <- names(fit$map)
-    for (obs_var in names(fit$map)) {
-      points(subset(fit$data, variable == obs_var, c(time, observed)), 
-        pch = pch_obs[obs_var], col = col_obs[obs_var])
-    }
-    matlines(out_transformed$time, out_transformed[-1])
-      legend("topright", inset=c(0.05, 0.05), legend=names(fit$map),
+    xlab = xlab, ylab = ylab, ...)
+  col_obs <- pch_obs <- 1:length(fit$map)
+  names(col_obs) <- names(pch_obs) <- names(fit$map)
+  for (obs_var in names(fit$map)) {
+    points(subset(fit$data, variable == obs_var, c(time, observed)), 
+    pch = pch_obs[obs_var], col = col_obs[obs_var])
+  }
+  matlines(out_transformed$time, out_transformed[-1])
+  if (legend == TRUE) {
+    legend("topright", inset=c(0.05, 0.05), legend=names(fit$map),
       col=col_obs, pch=pch_obs, lty=1:length(pch_obs))
+  }
 }
diff --git a/R/mkinresplot.R b/R/mkinresplot.R
new file mode 100644
index 00000000..bede5f78
--- /dev/null
+++ b/R/mkinresplot.R
@@ -0,0 +1,54 @@
+# $Id:  $

+

+# Copyright (C) 2008-2011 Katrin Lindenberger and Johannes Ranke

+# Contact: mkin-devel@lists.berlios.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/>

+

+mkinresplot <- function (object, obs_vars = vector(), 

+  xlab = "Time [days]", ylab = "Residual [% of applied radioactivity]",

+	maxabs = "auto", legend= TRUE, lpos = "topright", ...) 

+{

+	obs_vars_all <- as.character(unique(object$data$variable))

+

+  if (length(obs_vars) > 0){

+      vars <- intersect(obs_vars_all, obs_vars)	

+  } else vars <- obs_vars_all

+

+  residuals <- subset(object$data, variable %in% vars, residual)

+

+  if (maxabs == "auto") maxabs = max(abs(residuals), na.rm = TRUE)

+

+	col_obs <- pch_obs <- 1:length(vars)

+ 	names(col_obs) <- names(pch_obs) <- vars

+

+  plot(0,  xlab = xlab, ylab = ylab, 

+       xlim = c(0, 1.1 * max(object$data$time)), 

+       ylim = c(-1.2 * maxabs, 1.2 * maxabs), ...)

+

+	for(var in vars){

+		residuals_plot <- subset(object$data, variable == var, c("time", "residual"))

+		points(residuals_plot, pch = pch_obs[var], col = col_obs[var])

+	}

+

+  abline(h = 0, lty = 2)

+  title(paste("Residuals of mkin fit"), font.main = 1)

+

+ 	if (legend == TRUE) {

+    legend(lpos, inset = c(0.05, 0.05), legend = vars, 

+    col = col_obs, pch = pch_obs)

+  }

+}

diff --git a/R/mkinstart.R b/R/mkinstart.R
new file mode 100644
index 00000000..46e5bed7
--- /dev/null
+++ b/R/mkinstart.R
@@ -0,0 +1,23 @@
+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)

+  }

+}