From 3a7751c62e369547db2048b7fc89c21b862010e9 Mon Sep 17 00:00:00 2001 From: jranke Date: Mon, 18 Jun 2012 18:57:07 +0000 Subject: 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 --- R/endpoints.R | 94 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ R/mkinerrmin.R | 2 +- R/mkinfit.R | 98 +++------------------------------------------------------- 3 files changed, 100 insertions(+), 94 deletions(-) create mode 100644 R/endpoints.R (limited to 'R') 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) } -- cgit v1.2.1