From b5b446b718b15ccaae5b197e147fc1358f0f564e Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Fri, 6 Nov 2020 00:03:29 +0100 Subject: Fast analytical solutions for saemix, update.mmkin Also, use logit transformation for g and for solitary formation fractions, addressing #10. --- R/transform_odeparms.R | 61 ++++++++++++++++++++++++++++++++------------------ 1 file changed, 39 insertions(+), 22 deletions(-) (limited to 'R/transform_odeparms.R') diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R index 0a25ee8c..f21d31fc 100644 --- a/R/transform_odeparms.R +++ b/R/transform_odeparms.R @@ -5,19 +5,19 @@ #' constants and other parameters that can only take on positive values, a #' simple log transformation is used. For compositional parameters, such as the #' formations fractions that should always sum up to 1 and can not be negative, -#' the \code{\link{ilr}} transformation is used. +#' the [ilr] transformation is used. #' #' The transformation of sets of formation fractions is fragile, as it supposes #' the same ordering of the components in forward and backward transformation. -#' This is no problem for the internal use in \code{\link{mkinfit}}. +#' This is no problem for the internal use in [mkinfit]. #' #' @param parms Parameters of kinetic models as used in the differential #' equations. #' @param transparms Transformed parameters of kinetic models as used in the #' fitting procedure. -#' @param mkinmod The kinetic model of class \code{\link{mkinmod}}, containing +#' @param mkinmod The kinetic model of class [mkinmod], containing #' the names of the model variables that are needed for grouping the -#' formation fractions before \code{\link{ilr}} transformation, the parameter +#' formation fractions before [ilr] transformation, the parameter #' names and the information if the pathway to sink is included in the model. #' @param transform_rates Boolean specifying if kinetic rate constants should #' be transformed in the model specification used in the fitting for better @@ -28,11 +28,15 @@ #' @param transform_fractions Boolean specifying if formation fractions #' constants should be transformed in the model specification used in the #' fitting for better compliance with the assumption of normal distribution -#' of the estimator. The default (TRUE) is to do transformations. The g -#' parameter of the DFOP and HS models are also transformed, as they can also -#' be seen as compositional data. The transformation used for these -#' transformations is the \code{\link{ilr}} transformation. +#' of the estimator. The default (TRUE) is to do transformations. +#' The g parameter of the DFOP model is also seen as a fraction. +#' If a single fraction is transformed (g parameter of DFOP or only a single +#' target variable e.g. a single metabolite plus a pathway to sink), a +#' logistic transformation is used [stats::qlogis()]. In other cases, i.e. if +#' two or more formation fractions need to be transformed whose sum cannot +#' exceed one, the [ilr] transformation is used. #' @return A vector of transformed or backtransformed parameters +#' @importFrom stats plogis qlogis #' @author Johannes Ranke #' @examples #' @@ -91,8 +95,7 @@ #' #' @export transform_odeparms transform_odeparms <- function(parms, mkinmod, - transform_rates = TRUE, - transform_fractions = TRUE) + transform_rates = TRUE, transform_fractions = TRUE) { # We need the model specification for the names of the model # variables and the information on the sink @@ -119,8 +122,7 @@ transform_odeparms <- function(parms, mkinmod, N <- parms[grep("^N", names(parms))] transparms[names(N)] <- N - # Go through state variables and apply isometric logratio transformation to - # formation fractions if requested + # Go through state variables and transform formation fractions if requested mod_vars = names(spec) for (box in mod_vars) { f <- parms[grep(paste("^f", box, sep = "_"), names(parms))] @@ -128,9 +130,14 @@ transform_odeparms <- function(parms, mkinmod, if (length(f) > 0) { if(transform_fractions) { if (spec[[box]]$sink) { - trans_f <- ilr(c(f, 1 - sum(f))) - trans_f_names <- paste("f", box, "ilr", 1:length(trans_f), sep = "_") - transparms[trans_f_names] <- trans_f + if (length(f) == 1) { + trans_f_name <- paste("f", box, "qlogis", sep = "_") + transparms[trans_f_name] <- qlogis(f) + } else { + trans_f <- ilr(c(f, 1 - sum(f))) + trans_f_names <- paste("f", box, "ilr", 1:length(trans_f), sep = "_") + transparms[trans_f_names] <- trans_f + } } else { if (length(f) > 1) { trans_f <- ilr(f) @@ -161,7 +168,7 @@ transform_odeparms <- function(parms, mkinmod, if (!is.na(parms["g"])) { g <- parms["g"] if (transform_fractions) { - transparms["g_ilr"] <- ilr(c(g, 1 - g)) + transparms["g_qlogis"] <- qlogis(g) } else { transparms["g"] <- g } @@ -207,20 +214,25 @@ backtransform_odeparms <- function(transparms, mkinmod, N <- transparms[grep("^N", names(transparms))] parms[names(N)] <- N - # Go through state variables and apply inverse isometric logratio transformation + # Go through state variables and apply inverse transformations to formation fractions mod_vars = names(spec) for (box in mod_vars) { # Get the names as used in the model f_names = grep(paste("^f", box, sep = "_"), mkinmod$parms, value = TRUE) + # Get the formation fraction parameters trans_f = transparms[grep(paste("^f", box, sep = "_"), names(transparms))] if (length(trans_f) > 0) { if(transform_fractions) { - f <- invilr(trans_f) - if (spec[[box]]$sink) { - parms[f_names] <- f[1:length(f)-1] + if (any(grepl("qlogis", names(trans_f)))) { + parms[f_names] <- plogis(trans_f) } else { - parms[f_names] <- f + f <- invilr(trans_f) + if (spec[[box]]$sink) { + parms[f_names] <- f[1:length(f)-1] + } else { + parms[f_names] <- f + } } } else { parms[names(trans_f)] <- trans_f @@ -242,7 +254,12 @@ backtransform_odeparms <- function(transparms, mkinmod, } } - # DFOP parameter g is treated as a fraction + # DFOP parameter g is now transformed using qlogis + if (!is.na(transparms["g_qlogis"])) { + g_qlogis <- transparms["g_qlogis"] + parms["g"] <- plogis(g_qlogis) + } + # In earlier times we used ilr for g, so we keep this around if (!is.na(transparms["g_ilr"])) { g_ilr <- transparms["g_ilr"] parms["g"] <- invilr(g_ilr)[1] -- cgit v1.2.1