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

1 files changed, 328 insertions, 70 deletions

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)