diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/IORE.solution.R | 4 | ||||
| -rw-r--r-- | R/mkinfit.R | 4 | ||||
| -rw-r--r-- | R/mkinmod.R | 28 | ||||
| -rw-r--r-- | R/mkinpredict.R | 6 | ||||
| -rw-r--r-- | R/transform_odeparms.R | 14 |
5 files changed, 48 insertions, 8 deletions
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) {
|