aboutsummaryrefslogtreecommitdiff
path: root/R/transform_odeparms.R
diff options
context:
space:
mode:
Diffstat (limited to 'R/transform_odeparms.R')
-rw-r--r--R/transform_odeparms.R438
1 files changed, 258 insertions, 180 deletions
diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R
index f69f4ebd..28e58f87 100644
--- a/R/transform_odeparms.R
+++ b/R/transform_odeparms.R
@@ -1,180 +1,258 @@
-# Copyright (C) 2010-2014,2019 Johannes Ranke
-# Contact: jranke@uni-bremen.de
-
-# This file is part of the R package mkin
-
-# mkin is free software: you can redistribute it and/or modify it under the
-# terms of the GNU General Public License as published by the Free Software
-# Foundation, either version 3 of the License, or (at your option) any later
-# version.
-
-# This program is distributed in the hope that it will be useful, but WITHOUT
-# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
-# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-# details.
-
-# You should have received a copy of the GNU General Public License along with
-# this program. If not, see <http://www.gnu.org/licenses/>
-
-transform_odeparms <- function(parms, mkinmod,
- transform_rates = TRUE,
- transform_fractions = TRUE)
-{
- # We need the model specification for the names of the model
- # variables and the information on the sink
- spec = mkinmod$spec
-
- # Set up container for transformed parameters
- transparms <- numeric(0)
-
- # Do not transform initial values for state variables
- state.ini.optim <- parms[grep("_0$", names(parms))]
- transparms[names(state.ini.optim)] <- state.ini.optim
-
- # Log transformation for rate constants if requested
- k <- parms[grep("^k_", names(parms))]
- k__iore <- parms[grep("^k__iore_", names(parms))]
- k <- c(k, k__iore)
- if (length(k) > 0) {
- if(transform_rates) {
- transparms[paste0("log_", names(k))] <- log(k)
- } else transparms[names(k)] <- k
- }
-
- # Do not transform exponents in IORE models
- N <- parms[grep("^N", names(parms))]
- transparms[names(N)] <- N
-
- # Go through state variables and apply isometric logratio transformation to
- # formation fractions if requested
- mod_vars = names(spec)
- for (box in mod_vars) {
- f <- parms[grep(paste("^f", box, sep = "_"), names(parms))]
-
- 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
- } else {
- if (length(f) > 1) {
- trans_f <- ilr(f)
- trans_f_names <- paste("f", box, "ilr", 1:length(trans_f), sep = "_")
- transparms[trans_f_names] <- trans_f
- }
- }
- } else {
- transparms[names(f)] <- f
- }
- }
- }
-
- # Transform also FOMC parameters alpha and beta, DFOP and HS rates k1 and k2
- # and HS parameter tb as well as logistic model parameters kmax, k0 and r if
- # transformation of rates is requested
- for (pname in c("alpha", "beta", "k1", "k2", "tb", "kmax", "k0", "r")) {
- if (!is.na(parms[pname])) {
- if (transform_rates) {
- transparms[paste0("log_", pname)] <- log(parms[pname])
- } else {
- transparms[pname] <- parms[pname]
- }
- }
- }
-
- # DFOP parameter g is treated as a fraction
- if (!is.na(parms["g"])) {
- g <- parms["g"]
- if (transform_fractions) {
- transparms["g_ilr"] <- ilr(c(g, 1 - g))
- } else {
- transparms["g"] <- g
- }
- }
-
- return(transparms)
-}
-
-backtransform_odeparms <- function(transparms, mkinmod,
- transform_rates = TRUE,
- transform_fractions = TRUE)
-{
- # We need the model specification for the names of the model
- # variables and the information on the sink
- spec = mkinmod$spec
-
- # Set up container for backtransformed parameters
- parms <- numeric(0)
-
- # Do not transform initial values for state variables
- state.ini.optim <- transparms[grep("_0$", names(transparms))]
- parms[names(state.ini.optim)] <- state.ini.optim
-
- # Exponential transformation for rate constants
- if(transform_rates) {
- trans_k <- transparms[grep("^log_k_", names(transparms))]
- trans_k__iore <- transparms[grep("^log_k__iore_", names(transparms))]
- trans_k = c(trans_k, trans_k__iore)
- if (length(trans_k) > 0) {
- k_names <- gsub("^log_k", "k", names(trans_k))
- parms[k_names] <- exp(trans_k)
- }
- } else {
- trans_k <- transparms[grep("^k_", names(transparms))]
- parms[names(trans_k)] <- trans_k
- trans_k__iore <- transparms[grep("^k__iore_", names(transparms))]
- parms[names(trans_k__iore)] <- trans_k__iore
- }
-
- # Do not transform exponents in IORE models
- N <- transparms[grep("^N", names(transparms))]
- parms[names(N)] <- N
-
- # Go through state variables and apply inverse isometric logratio transformation
- 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]
- } else {
- parms[f_names] <- f
- }
- } else {
- parms[names(trans_f)] <- trans_f
- }
- }
- }
-
- # Transform parameters also for FOMC, DFOP, HS and logistic models
- for (pname in c("alpha", "beta", "k1", "k2", "tb", "kmax", "k0", "r")) {
- if (transform_rates) {
- pname_trans = paste0("log_", pname)
- if (!is.na(transparms[pname_trans])) {
- parms[pname] <- exp(transparms[pname_trans])
- }
- } else {
- if (!is.na(transparms[pname])) {
- parms[pname] <- transparms[pname]
- }
- }
- }
-
- # DFOP parameter g is treated as a fraction
- if (!is.na(transparms["g_ilr"])) {
- g_ilr <- transparms["g_ilr"]
- parms["g"] <- invilr(g_ilr)[1]
- }
- if (!is.na(transparms["g"])) {
- parms["g"] <- transparms["g"]
- }
-
- return(parms)
-}
-# vim: set ts=2 sw=2 expandtab:
+#' Functions to transform and backtransform kinetic parameters for fitting
+#'
+#' The transformations are intended to map parameters that should only take on
+#' restricted values to the full scale of real numbers. For kinetic rate
+#' constants and other paramters 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 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}}.
+#'
+#' @aliases transform_odeparms backtransform_odeparms
+#' @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
+#' the names of the model variables that are needed for grouping the
+#' formation fractions before \code{\link{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
+#' compliance with the assumption of normal distribution of the estimator. If
+#' TRUE, also alpha and beta parameters of the FOMC model are
+#' log-transformed, as well as k1 and k2 rate constants for the DFOP and HS
+#' models and the break point tb of the HS model.
+#' @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.
+#' @return A vector of transformed or backtransformed parameters with the same
+#' names as the original parameters.
+#' @author Johannes Ranke
+#' @examples
+#'
+#' SFO_SFO <- mkinmod(
+#' parent = list(type = "SFO", to = "m1", sink = TRUE),
+#' m1 = list(type = "SFO"))
+#' # Fit the model to the FOCUS example dataset D using defaults
+#' fit <- mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE)
+#' fit.s <- summary(fit)
+#' # Transformed and backtransformed parameters
+#' print(fit.s$par, 3)
+#' print(fit.s$bpar, 3)
+#'
+#' \dontrun{
+#' # Compare to the version without transforming rate parameters
+#' fit.2 <- mkinfit(SFO_SFO, FOCUS_2006_D, transform_rates = FALSE, quiet = TRUE)
+#' fit.2.s <- summary(fit.2)
+#' print(fit.2.s$par, 3)
+#' print(fit.2.s$bpar, 3)
+#' }
+#'
+#' initials <- fit$start$value
+#' names(initials) <- rownames(fit$start)
+#' transformed <- fit$start_transformed$value
+#' names(transformed) <- rownames(fit$start_transformed)
+#' transform_odeparms(initials, SFO_SFO)
+#' backtransform_odeparms(transformed, SFO_SFO)
+#'
+#' \dontrun{
+#' # The case of formation fractions
+#' SFO_SFO.ff <- mkinmod(
+#' parent = list(type = "SFO", to = "m1", sink = TRUE),
+#' m1 = list(type = "SFO"),
+#' use_of_ff = "max")
+#'
+#' fit.ff <- mkinfit(SFO_SFO.ff, FOCUS_2006_D, quiet = TRUE)
+#' fit.ff.s <- summary(fit.ff)
+#' print(fit.ff.s$par, 3)
+#' print(fit.ff.s$bpar, 3)
+#' initials <- c("f_parent_to_m1" = 0.5)
+#' transformed <- transform_odeparms(initials, SFO_SFO.ff)
+#' backtransform_odeparms(transformed, SFO_SFO.ff)
+#'
+#' # And without sink
+#' SFO_SFO.ff.2 <- mkinmod(
+#' parent = list(type = "SFO", to = "m1", sink = FALSE),
+#' m1 = list(type = "SFO"),
+#' use_of_ff = "max")
+#'
+#'
+#' fit.ff.2 <- mkinfit(SFO_SFO.ff.2, FOCUS_2006_D, quiet = TRUE)
+#' fit.ff.2.s <- summary(fit.ff.2)
+#' print(fit.ff.2.s$par, 3)
+#' print(fit.ff.2.s$bpar, 3)
+#' }
+#'
+#' @export transform_odeparms
+transform_odeparms <- function(parms, mkinmod,
+ transform_rates = TRUE,
+ transform_fractions = TRUE)
+{
+ # We need the model specification for the names of the model
+ # variables and the information on the sink
+ spec = mkinmod$spec
+
+ # Set up container for transformed parameters
+ transparms <- numeric(0)
+
+ # Do not transform initial values for state variables
+ state.ini.optim <- parms[grep("_0$", names(parms))]
+ transparms[names(state.ini.optim)] <- state.ini.optim
+
+ # Log transformation for rate constants if requested
+ k <- parms[grep("^k_", names(parms))]
+ k__iore <- parms[grep("^k__iore_", names(parms))]
+ k <- c(k, k__iore)
+ if (length(k) > 0) {
+ if(transform_rates) {
+ transparms[paste0("log_", names(k))] <- log(k)
+ } else transparms[names(k)] <- k
+ }
+
+ # Do not transform exponents in IORE models
+ N <- parms[grep("^N", names(parms))]
+ transparms[names(N)] <- N
+
+ # Go through state variables and apply isometric logratio transformation to
+ # formation fractions if requested
+ mod_vars = names(spec)
+ for (box in mod_vars) {
+ f <- parms[grep(paste("^f", box, sep = "_"), names(parms))]
+
+ 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
+ } else {
+ if (length(f) > 1) {
+ trans_f <- ilr(f)
+ trans_f_names <- paste("f", box, "ilr", 1:length(trans_f), sep = "_")
+ transparms[trans_f_names] <- trans_f
+ }
+ }
+ } else {
+ transparms[names(f)] <- f
+ }
+ }
+ }
+
+ # Transform also FOMC parameters alpha and beta, DFOP and HS rates k1 and k2
+ # and HS parameter tb as well as logistic model parameters kmax, k0 and r if
+ # transformation of rates is requested
+ for (pname in c("alpha", "beta", "k1", "k2", "tb", "kmax", "k0", "r")) {
+ if (!is.na(parms[pname])) {
+ if (transform_rates) {
+ transparms[paste0("log_", pname)] <- log(parms[pname])
+ } else {
+ transparms[pname] <- parms[pname]
+ }
+ }
+ }
+
+ # DFOP parameter g is treated as a fraction
+ if (!is.na(parms["g"])) {
+ g <- parms["g"]
+ if (transform_fractions) {
+ transparms["g_ilr"] <- ilr(c(g, 1 - g))
+ } else {
+ transparms["g"] <- g
+ }
+ }
+
+ return(transparms)
+}
+
+#' @describeIn transform_odeparms Backtransform the set of transformed parameters
+#' @export backtransform_odeparms
+backtransform_odeparms <- function(transparms, mkinmod,
+ transform_rates = TRUE,
+ transform_fractions = TRUE)
+{
+ # We need the model specification for the names of the model
+ # variables and the information on the sink
+ spec = mkinmod$spec
+
+ # Set up container for backtransformed parameters
+ parms <- numeric(0)
+
+ # Do not transform initial values for state variables
+ state.ini.optim <- transparms[grep("_0$", names(transparms))]
+ parms[names(state.ini.optim)] <- state.ini.optim
+
+ # Exponential transformation for rate constants
+ if(transform_rates) {
+ trans_k <- transparms[grep("^log_k_", names(transparms))]
+ trans_k__iore <- transparms[grep("^log_k__iore_", names(transparms))]
+ trans_k = c(trans_k, trans_k__iore)
+ if (length(trans_k) > 0) {
+ k_names <- gsub("^log_k", "k", names(trans_k))
+ parms[k_names] <- exp(trans_k)
+ }
+ } else {
+ trans_k <- transparms[grep("^k_", names(transparms))]
+ parms[names(trans_k)] <- trans_k
+ trans_k__iore <- transparms[grep("^k__iore_", names(transparms))]
+ parms[names(trans_k__iore)] <- trans_k__iore
+ }
+
+ # Do not transform exponents in IORE models
+ N <- transparms[grep("^N", names(transparms))]
+ parms[names(N)] <- N
+
+ # Go through state variables and apply inverse isometric logratio transformation
+ 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]
+ } else {
+ parms[f_names] <- f
+ }
+ } else {
+ parms[names(trans_f)] <- trans_f
+ }
+ }
+ }
+
+ # Transform parameters also for FOMC, DFOP, HS and logistic models
+ for (pname in c("alpha", "beta", "k1", "k2", "tb", "kmax", "k0", "r")) {
+ if (transform_rates) {
+ pname_trans = paste0("log_", pname)
+ if (!is.na(transparms[pname_trans])) {
+ parms[pname] <- exp(transparms[pname_trans])
+ }
+ } else {
+ if (!is.na(transparms[pname])) {
+ parms[pname] <- transparms[pname]
+ }
+ }
+ }
+
+ # DFOP parameter g is treated as a fraction
+ if (!is.na(transparms["g_ilr"])) {
+ g_ilr <- transparms["g_ilr"]
+ parms["g"] <- invilr(g_ilr)[1]
+ }
+ if (!is.na(transparms["g"])) {
+ parms["g"] <- transparms["g"]
+ }
+
+ return(parms)
+}
+# vim: set ts=2 sw=2 expandtab:

Contact - Imprint