From a2f772990127891f9596b79771832bc23777a89a Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Tue, 22 Apr 2014 19:08:09 +0200 Subject: Possibility to fit without parameter transformation --- DESCRIPTION | 4 ++-- R/mkinfit.R | 34 ++++++++++++++++++++++++++-------- R/transform_odeparms.R | 41 +++++++++++++++++++++++++---------------- TODO | 1 + man/mkinfit.Rd | 18 ++++++++++++++++++ man/transform_odeparms.Rd | 35 ++++++++++++++++++++++++++--------- 6 files changed, 98 insertions(+), 35 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 563431e..f348248 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,8 +2,8 @@ Package: mkin Type: Package Title: Routines for fitting kinetic models with one or more state variables to chemical degradation data -Version: 0.9-25 -Date: 2014-03-24 +Version: 0.9-26 +Date: 2014-04-22 Authors@R: c(person("Johannes", "Ranke", role = c("aut", "cre", "cph"), email = "jranke@uni-bremen.de"), person("Katrin", "Lindenberger", role = "ctb"), diff --git a/R/mkinfit.R b/R/mkinfit.R index 326f404..a623146 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -30,6 +30,8 @@ mkinfit <- function(mkinmod, observed, method.ode = "lsoda", method.modFit = "Marq", control.modFit = list(), + transform_rates = TRUE, + transform_fractions = TRUE, plot = FALSE, quiet = FALSE, err = NULL, weight = "none", scaleVar = FALSE, atol = 1e-8, rtol = 1e-10, n.outtimes = 100, @@ -83,7 +85,9 @@ mkinfit <- function(mkinmod, observed, if(is.null(names(state.ini))) names(state.ini) <- mod_vars # Transform initial parameter values for fitting - transparms.ini <- transform_odeparms(parms.ini, mod_vars) + transparms.ini <- transform_odeparms(parms.ini, mod_vars, + transform_rates = transform_rates, + transform_fractions = transform_fractions) # Parameters to be optimised: # Kinetic parameters in parms.ini whose names are not in fixed_parms @@ -156,7 +160,9 @@ mkinfit <- function(mkinmod, observed, odeparms <- c(P[(length(state.ini.optim) + 1):length(P)], parms.fixed) - parms <- backtransform_odeparms(odeparms, mod_vars) + parms <- backtransform_odeparms(odeparms, mod_vars, + transform_rates = transform_rates, + transform_fractions = transform_fractions) # Solve the system with current transformed parameter values out <- mkinpredict(mkinmod, parms, odeini, outtimes, @@ -241,6 +247,8 @@ mkinfit <- function(mkinmod, observed, # We need to return some more data for summary and plotting fit$solution_type <- solution_type + fit$transform_rates <- transform_rates + fit$transform_fractions <- transform_fractions # We also need the model for summary and plotting fit$mkinmod <- mkinmod @@ -252,19 +260,27 @@ mkinfit <- function(mkinmod, observed, # Collect initial parameter values in two dataframes fit$start <- data.frame(value = c(state.ini.optim, - backtransform_odeparms(parms.optim, mod_vars))) + backtransform_odeparms(parms.optim, mod_vars, + transform_rates = transform_rates, + transform_fractions = transform_fractions))) fit$start$type = c(rep("state", length(state.ini.optim)), rep("deparm", length(parms.optim))) fit$start$transformed = c(state.ini.optim, parms.optim) fit$fixed <- data.frame(value = c(state.ini.fixed, - backtransform_odeparms(parms.fixed, mod_vars))) + backtransform_odeparms(parms.fixed, mod_vars, + transform_rates = transform_rates, + transform_fractions = transform_fractions))) fit$fixed$type = c(rep("state", length(state.ini.fixed)), rep("deparm", length(parms.fixed))) - bparms.optim = backtransform_odeparms(fit$par, mod_vars) + bparms.optim = backtransform_odeparms(fit$par, mod_vars, + transform_rates = transform_rates, + transform_fractions = transform_fractions) bparms.fixed = backtransform_odeparms(c(state.ini.fixed, parms.fixed), - mod_vars) + mod_vars, + transform_rates = transform_rates, + transform_fractions = transform_fractions) bparms.all = c(bparms.optim, bparms.fixed) # Collect observed, predicted, residuals and weighting @@ -328,8 +344,10 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, par.lower <- par.upper <- object$par par.lower[pname] <- param[pname, "Lower"] par.upper[pname] <- param[pname, "Upper"] - blci[pname] <- backtransform_odeparms(par.lower, mod_vars)[pname] - buci[pname] <- backtransform_odeparms(par.upper, mod_vars)[pname] + blci[pname] <- backtransform_odeparms(par.lower, mod_vars, + object$transform_rates, object$transform_fractions)[pname] + buci[pname] <- backtransform_odeparms(par.upper, mod_vars, + object$transform_rates, object$transform_fractions)[pname] } bparam <- cbind(object$bparms.optim, blci, buci) dimnames(bparam) <- list(pnames, c("Estimate", "Lower", "Upper")) diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R index f56478f..31200c7 100644 --- a/R/transform_odeparms.R +++ b/R/transform_odeparms.R @@ -1,6 +1,4 @@ -# $Id$ - -# Copyright (C) 2010-2013 Johannes Ranke +# Copyright (C) 2010-2014 Johannes Ranke # Contact: jranke@uni-bremen.de # This file is part of the R package mkin @@ -18,17 +16,21 @@ # You should have received a copy of the GNU General Public License along with # this program. If not, see -transform_odeparms <- function(parms, mod_vars) { +transform_odeparms <- function(parms, mod_vars, + transform_rates = TRUE, + transform_fractions = TRUE) +{ # Set up container for transformed parameters transparms <- parms - # Log transformation for rate constants + # Log transformation for rate constants if requested index_k <- grep("^k_", names(transparms)) if (length(index_k) > 0) { - transparms[index_k] <- log(parms[index_k]) + if(transform_rates) transparms[index_k] <- log(parms[index_k]) + else transparms[index_k] <- parms[index_k] } - # Go through state variables and apply isotropic logratio transformation + # Go through state variables and apply isotropic logratio transformation if requested for (box in mod_vars) { indices_f <- grep(paste("^f", box, sep = "_"), names(parms)) f_names <- grep(paste("^f", box, sep = "_"), names(parms), value = TRUE) @@ -37,32 +39,38 @@ transform_odeparms <- function(parms, mod_vars) { f <- parms[indices_f] trans_f <- ilr(c(f, 1 - sum(f))) names(trans_f) <- f_names - transparms[indices_f] <- trans_f + if(transform_fractions) transparms[indices_f] <- trans_f + else transparms[indices_f] <- f } } - # Transform parameters also for FOMC, DFOP and HS models + # Transform parameters also for FOMC, DFOP and HS models if requested for (pname in c("alpha", "beta", "k1", "k2", "tb")) { if (!is.na(parms[pname])) { - transparms[pname] <- log(parms[pname]) + transparms[pname] <- ifelse(transform_rates, log(parms[pname]), parms[pname]) + transparms[pname] <- ifelse(transform_rates, log(parms[pname]), parms[pname]) } } if (!is.na(parms["g"])) { g <- parms["g"] - transparms["g"] <- ilr(c(g, 1- g)) + transparms["g"] <- ifelse(transform_fractions, ilr(c(g, 1 - g)), g) } return(transparms) } -backtransform_odeparms <- function(transparms, mod_vars) { +backtransform_odeparms <- function(transparms, mod_vars, + transform_rates = TRUE, + transform_fractions = TRUE) +{ # Set up container for backtransformed parameters parms <- transparms # Exponential transformation for rate constants index_k <- grep("^k_", names(parms)) if (length(index_k) > 0) { - parms[index_k] <- exp(transparms[index_k]) + if(transform_rates) parms[index_k] <- exp(transparms[index_k]) + else parms[index_k] <- transparms[index_k] } # Go through state variables and apply inverse isotropic logratio transformation @@ -73,19 +81,20 @@ backtransform_odeparms <- function(transparms, mod_vars) { if (n_paths > 0) { f <- invilr(transparms[indices_f])[1:n_paths] # We do not need the last component names(f) <- f_names - parms[indices_f] <- f + if(transform_fractions) parms[indices_f] <- f + else parms[indices_f] <- transparms[indices_f] } } # Transform parameters also for DFOP and HS models for (pname in c("alpha", "beta", "k1", "k2", "tb")) { if (!is.na(transparms[pname])) { - parms[pname] <- exp(transparms[pname]) + parms[pname] <- ifelse(transform_rates, exp(transparms[pname]), transparms[pname]) } } if (!is.na(transparms["g"])) { g <- transparms["g"] - parms["g"] <- invilr(g)[1] + parms["g"] <- ifelse(transform_fractions, invilr(g)[1], g) } return(parms) diff --git a/TODO b/TODO index affc587..944ad57 100644 --- a/TODO +++ b/TODO @@ -3,6 +3,7 @@ TODO for version 1.0 - Support model definitions without pathway to sink in combination with formation fractions - Complete the main package vignette named mkin to include a method description - Calculate pseudoDT50 values as recommended by FOCUS +- Improve formatting of differential equations in the summary Nice to have: - Calculate confidence intervals for DT50 and DT90 values when only one parameter is involved diff --git a/man/mkinfit.Rd b/man/mkinfit.Rd index 0c5d76a..839580a 100644 --- a/man/mkinfit.Rd +++ b/man/mkinfit.Rd @@ -23,6 +23,8 @@ mkinfit(mkinmod, observed, method.ode = "lsoda", method.modFit = "Marq", control.modFit = list(), + transform_rates = TRUE, + transform_fractions = TRUE, plot = FALSE, quiet = FALSE, err = NULL, weight = "none", scaleVar = FALSE, atol = 1e-8, rtol = 1e-10, n.outtimes = 100, @@ -102,6 +104,22 @@ mkinfit(mkinmod, observed, Additional arguments passed to the optimisation method used by \code{\link{modFit}}. } + \item{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. + } + \item{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. + } \item{plot}{ Should the observed values and the numerical solutions be plotted at each stage of the optimisation? diff --git a/man/transform_odeparms.Rd b/man/transform_odeparms.Rd index 999beaa..e251cf6 100644 --- a/man/transform_odeparms.Rd +++ b/man/transform_odeparms.Rd @@ -13,20 +13,37 @@ negative, the \code{\link{ilr}} transformation is used. } \usage{ -transform_odeparms(parms, mod_vars) -backtransform_odeparms(transparms, mod_vars) +transform_odeparms(parms, mod_vars, transform_rates = TRUE, transform_fractions = TRUE) +backtransform_odeparms(transparms, mod_vars, + transform_rates = TRUE, transform_fractions = TRUE) } \arguments{ \item{parms}{ - Parameters of kinetic models as used in the differential equations. -} + Parameters of kinetic models as used in the differential equations. + } \item{transparms}{ - Transformed parameters of kinetic models as used in the fitting procedure. -} + Transformed parameters of kinetic models as used in the fitting procedure. + } \item{mod_vars}{ - Names of the state variables in the kinetic model. These are used for - grouping the formation fractions before \code{\link{ilr}} transformation. -} + Names of the state variables in the kinetic model. These are used for + grouping the formation fractions before \code{\link{ilr}} transformation. + } + \item{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. + } + \item{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. + } } \value{ A vector of transformed or backtransformed parameters with the same names -- cgit v1.2.1