First working fits with IORE model

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) {