Put the calculation of dissipation times, formation fractions and SFORB eigenvalues into a separate function

Some testing of this functionality in test.R



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

3 files changed, 100 insertions, 94 deletions

diff --git a/R/endpoints.R b/R/endpoints.R
new file mode 100644
index 0000000..98e290d
--- /dev/null
+++ b/R/endpoints.R
@@ -0,0 +1,94 @@
+endpoints <- function(fit, pseudoDT50 = FALSE) {

+  # Calculate dissipation times DT50 and DT90 and, if necessary, formation

+  # fractions and SFORB eigenvalues from optimised parameters

+  ep <- list()

+  obs_vars <- fit$obs_vars

+  parms.all <- fit$parms.all

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

+			    DT90 = rep(NA, length(obs_vars)), 

+    row.names = obs_vars)

+  ep$ff <- vector()

+  ep$SFORB <- vector()

+  for (obs_var in obs_vars) {

+    type = names(fit$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)

+      {

+        ep$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 == "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

+      }

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

+      }

+      DT50 <- DTx(50)

+      DT90 <- DTx(90)

+    }

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

+      {

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

+      }

+

+      # Return the eigenvalues for comparison with DFOP rate constants

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

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

+    }

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

+  }

+  return(ep)

+}

diff --git a/R/mkinerrmin.R b/R/mkinerrmin.R
index 5415b41..c278647 100644
--- a/R/mkinerrmin.R
+++ b/R/mkinerrmin.R
@@ -76,5 +76,5 @@ mkinerrmin <- function(fit, alpha = 0.05)
     errmin[obs_var, c("err.min", "n.optim", "df")] <- errmin.tmp

   }

 

-	return(errmin)

+  return(errmin)

 }

diff --git a/R/mkinfit.R b/R/mkinfit.R
index b2641e6..d3ee5e7 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -189,97 +189,8 @@ mkinfit <- function(mkinmod, observed,
     value = c(state.ini.fixed, parms.fixed))

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

 

-

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

 

-  # Calculate dissipation times DT50 and DT90 and, if necessary, formation

-  # fractions and SFORB eigenvalues from optimised parameters

-  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

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

-      }

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

-      }

-      DT50 <- DTx(50)

-      DT90 <- DTx(90)

-    }

-    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

   data <- merge(fit$observed, fit$predicted, by = c("time", "name"))

   names(data) <- c("time", "variable", "observed", "predicted")

@@ -342,10 +253,11 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, ...) {
   ans$errmin <- mkinerrmin(object, alpha = 0.05)

 

   ans$parms.all <- object$parms.all

-  if (!is.null(object$ff))

-    ans$ff <- object$ff

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

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

+  ep <- endpoints(object)

+  if (!is.null(ep$ff))

+    ans$ff <- ep$ff

+  if(distimes) ans$distimes <- ep$distimes

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

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

   return(ans)  

 }