From d518c4cfa22994db5ba81a6b01c6cb4c4186872e Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 15 Apr 2020 20:42:29 +0200 Subject: Adapt endpoint() to also work for nlme.mmkin objects --- NEWS.md | 2 +- R/endpoints.R | 39 +++++++++++++++++++++++++++++++-------- R/nlme.mmkin.R | 5 +++++ docs/news/index.html | 2 +- docs/reference/endpoints.html | 29 ++++++++++++++++++++++------- docs/reference/nlme.mmkin.html | 29 ++++++++++++++++++++++++++--- man/endpoints.Rd | 12 ++++++++++-- man/nlme.mmkin.Rd | 5 +++++ 8 files changed, 101 insertions(+), 22 deletions(-) diff --git a/NEWS.md b/NEWS.md index a4976488..37800e6d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,6 +1,6 @@ # mkin 0.9.49.10 (unreleased) -- 'nlme.mmkin': An nlme method for mmkin row objects and an associated class with plot method +- 'nlme.mmkin': An nlme method for mmkin row objects and an associated S3 class with print, plot, anova and endpoint methods - 'mean_degparms, nlme_data, nlme_function': Three new functions to facilitate building nlme models from mmkin row objects diff --git a/R/endpoints.R b/R/endpoints.R index f7ee483a..586ef9ff 100644 --- a/R/endpoints.R +++ b/R/endpoints.R @@ -7,9 +7,13 @@ #' 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}}. +#' @param fit An object of class \code{\link{mkinfit}} or +#' \code{\link{nlme.mmkin}} #' @importFrom stats optimize -#' @return A list with the components mentioned above. +#' @return A list with a matrix of dissipation times named distimes, +#' and, if applicable, a vector of formation fractions named ff +#' and, if the SFORB model was in use, a vector of eigenvalues +#' of these SFORB models, equivalent to DFOP rate constants #' @note The function is used internally by \code{\link{summary.mkinfit}}. #' @author Johannes Ranke #' @keywords manip @@ -17,6 +21,10 @@ #' #' fit <- mkinfit("FOMC", FOCUS_2006_C, quiet = TRUE) #' endpoints(fit) +#' \dontrun{ +#' fit_2 <- mkinfit("SFORB", FOCUS_2006_C, quiet = TRUE) +#' endpoints(fit_2) +#' } #' #' @export endpoints <- function(fit) { @@ -25,8 +33,22 @@ endpoints <- function(fit) { # 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) + if (inherits(fit, "mkinfit")) { + mkinmod <- fit$mkinmod + parms.all <- c(fit$bparms.optim, fit$bparms.fixed) + } else { + if (inherits(fit, "nlme.mmkin")) { + mkinmod <- fit$mmkin_orig[[1]]$mkinmod + bparms.optim <- backtransform_odeparms(fit$coefficients$fixed, + mkinmod, + transform_rates = fit$mmkin_orig[[1]]$transform_rates, + transform_fractions = fit$mmkin_orig[[1]]$transform_fractions) + parms.all <- c(bparms.optim, fit$bparms.fixed) + } else { + stop("Only implemented for mkinfit and nlme.mmkin objects") + } + } + obs_vars <- names(mkinmod$spec) ep$ff <- vector() ep$SFORB <- vector() ep$distimes <- data.frame( @@ -34,7 +56,7 @@ endpoints <- function(fit) { DT90 = rep(NA, length(obs_vars)), row.names = obs_vars) for (obs_var in obs_vars) { - type = names(fit$mkinmod$map[[obs_var]])[1] + type = names(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) @@ -56,7 +78,7 @@ endpoints <- function(fit) { k_tot = sum(parms.all[k_names]) DT50 = log(2)/k_tot DT90 = log(10)/k_tot - if (fit$mkinmod$use_of_ff == "min") { + if (mkinmod$use_of_ff == "min" && length(obs_vars) > 1) { for (k_name in k_names) { ep$ff[[sub("k_", "", k_name)]] = parms.all[[k_name]] / k_tot @@ -78,7 +100,7 @@ endpoints <- function(fit) { 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]] + source_name = mkinmod$map[[1]][[1]] c0 = parms.all[paste(source_name, "0", sep = "_")] alpha = 1 / (n - 1) beta = (c0^(1 - n))/(k * (n - 1)) @@ -86,7 +108,7 @@ endpoints <- function(fit) { 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") { + if (mkinmod$use_of_ff == "min") { for (k_name in k_names) { ep$ff[[sub("k_", "", k_name)]] = parms.all[[k_name]] / k_tot @@ -196,6 +218,7 @@ endpoints <- function(fit) { } ep$distimes[obs_var, c("DT50", "DT90")] = c(DT50, DT90) } + if (length(ep$ff) == 0) ep$ff <- NULL if (length(ep$SFORB) == 0) ep$SFORB <- NULL return(ep) } diff --git a/R/nlme.mmkin.R b/R/nlme.mmkin.R index 6e3467ed..1d6c2e75 100644 --- a/R/nlme.mmkin.R +++ b/R/nlme.mmkin.R @@ -31,8 +31,10 @@ #' function(x) subset(x$data[c("name", "time", "value")], name == "parent")) #' f <- mmkin("SFO", ds, quiet = TRUE, cores = 1) #' library(nlme) +#' endpoints(f[[1]]) #' f_nlme <- nlme(f) #' print(f_nlme) +#' endpoints(f_nlme) #' f_nlme_2 <- nlme(f, start = c(parent_0 = 100, log_k_parent_sink = 0.1)) #' update(f_nlme_2, random = parent_0 ~ 1) #' \dontrun{ @@ -77,6 +79,9 @@ #' #' anova(f_nlme_dfop_sfo, f_nlme_fomc_sfo, f_nlme_sfo_sfo) #' anova(f_nlme_dfop_sfo, f_nlme_sfo_sfo) # if we ignore FOMC +#' +#' endpoints(f_nlme_sfo_sfo) +#' endpoints(f_nlme_dfop_sfo) #' } # Code inspired by nlme.nlsList nlme.mmkin <- function(model, data = sys.frame(sys.parent()), diff --git a/docs/news/index.html b/docs/news/index.html index 5fd97344..cb18664e 100644 --- a/docs/news/index.html +++ b/docs/news/index.html @@ -134,7 +134,7 @@ mkin 0.9.49.10 (unreleased) Unreleased