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 --- R/mkinfit.R | 34 ++++++++++++++++++++++++++-------- R/transform_odeparms.R | 41 +++++++++++++++++++++++++---------------- 2 files changed, 51 insertions(+), 24 deletions(-) (limited to 'R') 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) -- cgit v1.2.1