From 0a3eb0893cb4bd1b12f07a79069d1c7f5e991495 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Fri, 25 Oct 2019 00:37:42 +0200 Subject: Use roxygen for functions and methods --- R/endpoints.R | 377 +++++++++++++++++++++++++++++++--------------------------- 1 file changed, 199 insertions(+), 178 deletions(-) (limited to 'R/endpoints.R') diff --git a/R/endpoints.R b/R/endpoints.R index f8a44c4d..14beadea 100644 --- a/R/endpoints.R +++ b/R/endpoints.R @@ -1,178 +1,199 @@ -endpoints <- function(fit) { - # Calculate dissipation times DT50 and DT90 and formation - # fractions as well as SFORB eigenvalues from optimised parameters - # Additional DT50 values are calculated from the FOMC DT90 and k1 and k2 from - # HS and DFOP, as well as from Eigenvalues b1 and b2 of any SFORB models - ep <- list() - obs_vars <- fit$obs_vars - parms.all <- c(fit$bparms.optim, fit$bparms.fixed) - ep$ff <- vector() - ep$SFORB <- vector() - ep$distimes <- data.frame(DT50 = rep(NA, length(obs_vars)), - DT90 = rep(NA, length(obs_vars)), - row.names = obs_vars) - for (obs_var in obs_vars) { - type = names(fit$mkinmod$map[[obs_var]])[1] - - # Get formation fractions if directly fitted, and calculate remaining fraction to sink - f_names = grep(paste("^f", obs_var, sep = "_"), names(parms.all), value=TRUE) - if (length(f_names) > 0) { - f_values = parms.all[f_names] - f_to_sink = 1 - sum(f_values) - names(f_to_sink) = ifelse(type == "SFORB", - paste(obs_var, "free", "sink", sep = "_"), - paste(obs_var, "sink", sep = "_")) - for (f_name in f_names) { - ep$ff[[sub("f_", "", sub("_to_", "_", f_name))]] = f_values[[f_name]] - } - ep$ff = append(ep$ff, f_to_sink) - } - - # Get the rest - 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 - if (fit$mkinmod$use_of_ff == "min") { - 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) - DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011 - ep$distimes[obs_var, c("DT50back")] = DT50_back - } - if (type == "IORE") { - k_names = grep(paste("^k__iore", obs_var, sep="_"), names(parms.all), value=TRUE) - k_tot = sum(parms.all[k_names]) - # From the NAFTA kinetics guidance, p. 5 - n = parms.all[paste("N", obs_var, sep = "_")] - k = k_tot - # Use the initial concentration of the parent compound - source_name = fit$mkinmod$map[[1]][[1]] - c0 = parms.all[paste(source_name, "0", sep = "_")] - alpha = 1 / (n - 1) - beta = (c0^(1 - n))/(k * (n - 1)) - DT50 = beta * (2^(1/alpha) - 1) - DT90 = beta * (10^(1/alpha) - 1) - DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011 - ep$distimes[obs_var, c("DT50back")] = DT50_back - if (fit$mkinmod$use_of_ff == "min") { - for (k_name in k_names) - { - ep$ff[[sub("k_", "", k_name)]] = parms.all[[k_name]] / k_tot - } - } - } - if (type == "DFOP") { - k1 = parms.all["k1"] - k2 = parms.all["k2"] - g = parms.all["g"] - f <- function(log_t, x) { - t <- exp(log_t) - fraction <- g * exp( - k1 * t) + (1 - g) * exp( - k2 * t) - (fraction - (1 - x/100))^2 - } - DT50_k1 = log(2)/k1 - DT50_k2 = log(2)/k2 - DT90_k1 = log(10)/k1 - DT90_k2 = log(10)/k2 - - DT50 <- try(exp(optimize(f, c(log(DT50_k1), log(DT50_k2)), x=50)$minimum), - silent = TRUE) - DT90 <- try(exp(optimize(f, c(log(DT90_k1), log(DT90_k2)), x=90)$minimum), - silent = TRUE) - if (inherits(DT50, "try-error")) DT50 = NA - if (inherits(DT90, "try-error")) DT90 = NA - - ep$distimes[obs_var, c("DT50_k1")] = DT50_k1 - ep$distimes[obs_var, c("DT50_k2")] = DT50_k2 - } - 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) - DT50_k1 = log(2)/k1 - DT50_k2 = log(2)/k2 - ep$distimes[obs_var, c("DT50_k1")] = DT50_k1 - ep$distimes[obs_var, c("DT50_k2")] = DT50_k2 - } - 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 - - DT50_b1 = log(2)/b1 - DT50_b2 = log(2)/b2 - DT90_b1 = log(10)/b1 - DT90_b2 = log(10)/b2 - - 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(log_t) (SFORB_fraction(exp(log_t)) - 0.5)^2 - log_DT50 <- try(optimize(f_50, c(log(DT50_b1), log(DT50_b2)))$minimum, - silent = TRUE) - f_90 <- function(log_t) (SFORB_fraction(exp(log_t)) - 0.1)^2 - log_DT90 <- try(optimize(f_90, c(log(DT90_b1), log(DT90_b2)))$minimum, - silent = TRUE) - - DT50 = if (inherits(log_DT50, "try-error")) NA - else exp(log_DT50) - DT90 = if (inherits(log_DT90, "try-error")) NA - else exp(log_DT90) - - 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(paste("DT50", obs_var, "b1", sep = "_"))] = DT50_b1 - ep$distimes[obs_var, c(paste("DT50", obs_var, "b2", sep = "_"))] = DT50_b2 - } - if (type == "logistic") { - # FOCUS kinetics (2014) p. 67 - kmax = parms.all["kmax"] - k0 = parms.all["k0"] - r = parms.all["r"] - DT50 = (1/r) * log(1 - ((kmax/k0) * (1 - 2^(r/kmax)))) - DT90 = (1/r) * log(1 - ((kmax/k0) * (1 - 10^(r/kmax)))) - - DT50_k0 = log(2)/k0 - DT50_kmax = log(2)/kmax - ep$distimes[obs_var, c("DT50_k0")] = DT50_k0 - ep$distimes[obs_var, c("DT50_kmax")] = DT50_kmax - } - ep$distimes[obs_var, c("DT50", "DT90")] = c(DT50, DT90) - } - return(ep) -} +#' Function to calculate endpoints for further use from kinetic models fitted +#' with mkinfit +#' +#' This function calculates DT50 and DT90 values as well as formation fractions +#' from kinetic models fitted with mkinfit. If the SFORB model was specified +#' for one of the parents or metabolites, the Eigenvalues are returned. These +#' are equivalent to the rate constantes of the DFOP model, but with the +#' advantage that the SFORB model can also be used for metabolites. +#' +#' @param fit An object of class \code{\link{mkinfit}}. +#' @importFrom stats optimize +#' @return A list with the components mentioned above. +#' @note The function is used internally by \code{\link{summary.mkinfit}}. +#' @author Johannes Ranke +#' @keywords manip +#' @examples +#' +#' fit <- mkinfit("FOMC", FOCUS_2006_C, quiet = TRUE) +#' endpoints(fit) +#' +#' @export +endpoints <- function(fit) { + # Calculate dissipation times DT50 and DT90 and formation + # fractions as well as SFORB eigenvalues from optimised parameters + # Additional DT50 values are calculated from the FOMC DT90 and k1 and k2 from + # HS and DFOP, as well as from Eigenvalues b1 and b2 of any SFORB models + ep <- list() + obs_vars <- fit$obs_vars + parms.all <- c(fit$bparms.optim, fit$bparms.fixed) + ep$ff <- vector() + ep$SFORB <- vector() + ep$distimes <- data.frame(DT50 = rep(NA, length(obs_vars)), + DT90 = rep(NA, length(obs_vars)), + row.names = obs_vars) + for (obs_var in obs_vars) { + type = names(fit$mkinmod$map[[obs_var]])[1] + + # Get formation fractions if directly fitted, and calculate remaining fraction to sink + f_names = grep(paste("^f", obs_var, sep = "_"), names(parms.all), value=TRUE) + if (length(f_names) > 0) { + f_values = parms.all[f_names] + f_to_sink = 1 - sum(f_values) + names(f_to_sink) = ifelse(type == "SFORB", + paste(obs_var, "free", "sink", sep = "_"), + paste(obs_var, "sink", sep = "_")) + for (f_name in f_names) { + ep$ff[[sub("f_", "", sub("_to_", "_", f_name))]] = f_values[[f_name]] + } + ep$ff = append(ep$ff, f_to_sink) + } + + # Get the rest + 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 + if (fit$mkinmod$use_of_ff == "min") { + 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) + DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011 + ep$distimes[obs_var, c("DT50back")] = DT50_back + } + if (type == "IORE") { + k_names = grep(paste("^k__iore", obs_var, sep="_"), names(parms.all), value=TRUE) + k_tot = sum(parms.all[k_names]) + # From the NAFTA kinetics guidance, p. 5 + n = parms.all[paste("N", obs_var, sep = "_")] + k = k_tot + # Use the initial concentration of the parent compound + source_name = fit$mkinmod$map[[1]][[1]] + c0 = parms.all[paste(source_name, "0", sep = "_")] + alpha = 1 / (n - 1) + beta = (c0^(1 - n))/(k * (n - 1)) + DT50 = beta * (2^(1/alpha) - 1) + DT90 = beta * (10^(1/alpha) - 1) + DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011 + ep$distimes[obs_var, c("DT50back")] = DT50_back + if (fit$mkinmod$use_of_ff == "min") { + for (k_name in k_names) + { + ep$ff[[sub("k_", "", k_name)]] = parms.all[[k_name]] / k_tot + } + } + } + if (type == "DFOP") { + k1 = parms.all["k1"] + k2 = parms.all["k2"] + g = parms.all["g"] + f <- function(log_t, x) { + t <- exp(log_t) + fraction <- g * exp( - k1 * t) + (1 - g) * exp( - k2 * t) + (fraction - (1 - x/100))^2 + } + DT50_k1 = log(2)/k1 + DT50_k2 = log(2)/k2 + DT90_k1 = log(10)/k1 + DT90_k2 = log(10)/k2 + + DT50 <- try(exp(optimize(f, c(log(DT50_k1), log(DT50_k2)), x=50)$minimum), + silent = TRUE) + DT90 <- try(exp(optimize(f, c(log(DT90_k1), log(DT90_k2)), x=90)$minimum), + silent = TRUE) + if (inherits(DT50, "try-error")) DT50 = NA + if (inherits(DT90, "try-error")) DT90 = NA + + ep$distimes[obs_var, c("DT50_k1")] = DT50_k1 + ep$distimes[obs_var, c("DT50_k2")] = DT50_k2 + } + 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) + DT50_k1 = log(2)/k1 + DT50_k2 = log(2)/k2 + ep$distimes[obs_var, c("DT50_k1")] = DT50_k1 + ep$distimes[obs_var, c("DT50_k2")] = DT50_k2 + } + 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 + + DT50_b1 = log(2)/b1 + DT50_b2 = log(2)/b2 + DT90_b1 = log(10)/b1 + DT90_b2 = log(10)/b2 + + 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(log_t) (SFORB_fraction(exp(log_t)) - 0.5)^2 + log_DT50 <- try(optimize(f_50, c(log(DT50_b1), log(DT50_b2)))$minimum, + silent = TRUE) + f_90 <- function(log_t) (SFORB_fraction(exp(log_t)) - 0.1)^2 + log_DT90 <- try(optimize(f_90, c(log(DT90_b1), log(DT90_b2)))$minimum, + silent = TRUE) + + DT50 = if (inherits(log_DT50, "try-error")) NA + else exp(log_DT50) + DT90 = if (inherits(log_DT90, "try-error")) NA + else exp(log_DT90) + + 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(paste("DT50", obs_var, "b1", sep = "_"))] = DT50_b1 + ep$distimes[obs_var, c(paste("DT50", obs_var, "b2", sep = "_"))] = DT50_b2 + } + if (type == "logistic") { + # FOCUS kinetics (2014) p. 67 + kmax = parms.all["kmax"] + k0 = parms.all["k0"] + r = parms.all["r"] + DT50 = (1/r) * log(1 - ((kmax/k0) * (1 - 2^(r/kmax)))) + DT90 = (1/r) * log(1 - ((kmax/k0) * (1 - 10^(r/kmax)))) + + DT50_k0 = log(2)/k0 + DT50_kmax = log(2)/kmax + ep$distimes[obs_var, c("DT50_k0")] = DT50_k0 + ep$distimes[obs_var, c("DT50_kmax")] = DT50_kmax + } + ep$distimes[obs_var, c("DT50", "DT90")] = c(DT50, DT90) + } + return(ep) +} -- cgit v1.2.1