From 316f1360d1e12eefe491f86b0bbcd6dcf091c736 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Sat, 12 Jul 2014 08:32:08 +0200 Subject: First working fits with IORE model --- R/IORE.solution.R | 4 ++++ R/mkinfit.R | 4 +++- R/mkinmod.R | 28 ++++++++++++++++++++++------ R/mkinpredict.R | 6 ++++++ R/transform_odeparms.R | 14 +++++++++++++- 5 files changed, 48 insertions(+), 8 deletions(-) create mode 100644 R/IORE.solution.R (limited to 'R') diff --git a/R/IORE.solution.R b/R/IORE.solution.R new file mode 100644 index 00000000..9546ce56 --- /dev/null +++ b/R/IORE.solution.R @@ -0,0 +1,4 @@ +IORE.solution <- function(t, parent.0, k.iore, N) +{ + parent = (parent.0^(1 - N) - (1 - N) * k.iore * t)^(1/(1 - N)) +} diff --git a/R/mkinfit.R b/R/mkinfit.R index b7ca1d74..26fc0e6b 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -90,13 +90,15 @@ mkinfit <- function(mkinmod, observed, defaultpar.names <- setdiff(mkinmod$parms, names(parms.ini)) for (parmname in defaultpar.names) { # Default values for rate constants, depending on the parameterisation - if (substr(parmname, 1, 2) == "k_") { + if (grepl("^k", parmname)) { parms.ini[parmname] = 0.1 + k_salt k_salt = k_salt + 1e-4 } # Default values for rate constants for reversible binding if (grepl("free_bound$", parmname)) parms.ini[parmname] = 0.1 if (grepl("bound_free$", parmname)) parms.ini[parmname] = 0.02 + # Default values for IORE exponents + if (grepl("^N", parmname)) parms.ini[parmname] = 1 # Default values for the FOMC, DFOP and HS models if (parmname == "alpha") parms.ini[parmname] = 1 if (parmname == "beta") parms.ini[parmname] = 10 diff --git a/R/mkinmod.R b/R/mkinmod.R index eb290719..321887fc 100644 --- a/R/mkinmod.R +++ b/R/mkinmod.R @@ -41,9 +41,8 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL) map <- list() # }}} - # Do not return a coefficient matrix mat when FOMC, DFOP or HS are used for - # the parent compound {{{ - if(spec[[1]]$type %in% c("FOMC", "DFOP", "HS")) { + # Do not return a coefficient matrix mat when FOMC, IORE, DFOP or HS is used for the parent {{{ + if(spec[[1]]$type %in% c("FOMC", "IORE", "DFOP", "HS")) { mat = FALSE } else mat = TRUE #}}} @@ -55,8 +54,8 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL) # Check the type component of the compartment specification {{{ if(is.null(spec[[varname]]$type)) stop( "Every part of the model specification must be a list containing a type component") - if(!spec[[varname]]$type %in% c("SFO", "FOMC", "DFOP", "HS", "SFORB")) stop( - "Available types are SFO, FOMC, DFOP, HS and SFORB only") + if(!spec[[varname]]$type %in% c("SFO", "FOMC", "IORE", "DFOP", "HS", "SFORB")) stop( + "Available types are SFO, FOMC, IORE, DFOP, HS and SFORB only") if(spec[[varname]]$type %in% c("FOMC", "DFOP", "HS") & match(varname, obs_vars) != 1) { stop(paste("Types FOMC, DFOP and HS are only implemented for the first compartment,", "which is assumed to be the source compartment")) @@ -65,6 +64,7 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL) new_boxes <- switch(spec[[varname]]$type, SFO = varname, FOMC = varname, + IORE = varname, DFOP = varname, HS = varname, SFORB = paste(varname, c("free", "bound"), sep = "_") @@ -85,20 +85,36 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL) # Turn on sink if this is not explicitly excluded by the user by # specifying sink=FALSE if(is.null(spec[[varname]]$sink)) spec[[varname]]$sink <- TRUE - if(spec[[varname]]$type %in% c("SFO", "SFORB")) { # {{{ Add SFO or SFORB decline term + if(spec[[varname]]$type %in% c("SFO", "IORE", "SFORB")) { # {{{ Add decline term if (use_of_ff == "min") { # Minimum use of formation fractions if(spec[[varname]]$sink) { # If sink is required, add first-order sink term k_compound_sink <- paste("k", box_1, "sink", sep = "_") + if(spec[[varname]]$type == "IORE") { + k_compound_sink <- paste("k.iore", box_1, "sink", sep = "_") + } parms <- c(parms, k_compound_sink) decline_term <- paste(k_compound_sink, "*", box_1) + if(spec[[varname]]$type == "IORE") { + N <- paste("N", box_1, sep = "_") + parms <- c(parms, N) + decline_term <- paste0(decline_term, "^", N) + } } else { # otherwise no decline term needed here decline_term = "0" } } else { k_compound <- paste("k", box_1, sep = "_") + if(spec[[varname]]$type == "IORE") { + k_compound <- paste("k.iore", box_1, sep = "_") + } parms <- c(parms, k_compound) decline_term <- paste(k_compound, "*", box_1) + if(spec[[varname]]$type == "IORE") { + N <- paste("N", box_1, sep = "_") + parms <- c(parms, N) + decline_term <- paste0(decline_term, "^", N) + } } } #}}} if(spec[[varname]]$type == "FOMC") { # {{{ Add FOMC decline term diff --git a/R/mkinpredict.R b/R/mkinpredict.R index 14efc568..43b5d2f3 100644 --- a/R/mkinpredict.R +++ b/R/mkinpredict.R @@ -50,6 +50,12 @@ mkinpredict <- function(mkinmod, odeparms, odeini, FOMC = FOMC.solution(outtimes, evalparse(parent.name), evalparse("alpha"), evalparse("beta")), + IORE = IORE.solution(outtimes, + evalparse(parent.name), + ifelse(mkinmod$use_of_ff == "min", + evalparse(paste("k.iore", parent.name, "sink", sep="_")), + evalparse(paste("k.iore", parent.name, sep="_"))), + evalparse("N_parent")), DFOP = DFOP.solution(outtimes, evalparse(parent.name), evalparse("k1"), evalparse("k2"), diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R index 912a5c0a..a36b7eae 100644 --- a/R/transform_odeparms.R +++ b/R/transform_odeparms.R @@ -33,6 +33,8 @@ transform_odeparms <- function(parms, mkinmod, # 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) @@ -40,6 +42,10 @@ transform_odeparms <- function(parms, mkinmod, 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) @@ -98,8 +104,10 @@ backtransform_odeparms <- function(transparms, mkinmod, # 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)) + k_names <- gsub("^log_k", "k", names(trans_k)) parms[k_names] <- exp(trans_k) } } else { @@ -107,6 +115,10 @@ backtransform_odeparms <- function(transparms, mkinmod, parms[names(trans_k)] <- trans_k } + # 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) { -- cgit v1.2.1