Support formation fractions without sink pathway, updates

4 files changed, 215 insertions, 102 deletions

diff --git a/R/endpoints.R b/R/endpoints.R
index a5bebd71..676831a0 100644
--- a/R/endpoints.R
+++ b/R/endpoints.R
@@ -1,8 +1,8 @@
 endpoints <- function(fit) {

-  # Calculate dissipation times DT50 and DT90 and, if necessary, formation

-  # fractions and SFORB eigenvalues from optimised parameters

+  # Calculate dissipation times DT50 and DT90 and formation

+  # fractions as well as SFORB eigenvalues from optimised parameters

   # Additional DT50 values are calculated from the FOMC DT90 and k1 and k2 from HS and DFOP,

-  # as well as from Eigenvalues b1 and b2 of the SFORB model

+  # as well as from Eigenvalues b1 and b2 of any SFORB models

   ep <- list()

   obs_vars <- fit$obs_vars

   parms.all <- fit$bparms.ode

@@ -16,15 +16,17 @@ endpoints <- function(fit) {
 

     # Get formation fractions if directly fitted, and calculate remaining fraction to sink

     f_names = grep(paste("f", obs_var, sep = "_"), names(parms.all), value=TRUE)

-    f_values = parms.all[f_names]

-    f_to_sink = 1 - sum(f_values)

-    names(f_to_sink) = ifelse(type == "SFORB", 

-                            paste(obs_var, "free", "sink", sep = "_"), 

-                            paste(obs_var, "sink", sep = "_"))

-    for (f_name in f_names) {

-      ep$ff[[sub("f_", "", sub("_to_", "_", f_name))]] = f_values[[f_name]]

+    if (length(f_names) > 0) {

+      f_values = parms.all[f_names]

+      f_to_sink = 1 - sum(f_values)

+      names(f_to_sink) = ifelse(type == "SFORB", 

+                              paste(obs_var, "free", "sink", sep = "_"), 

+                              paste(obs_var, "sink", sep = "_"))

+      for (f_name in f_names) {

+        ep$ff[[sub("f_", "", sub("_to_", "_", f_name))]] = f_values[[f_name]]

+      }

+      ep$ff = append(ep$ff, f_to_sink)

     }

-    ep$ff = append(ep$ff, f_to_sink)

 

     # Get the rest

     if (type == "SFO") {

diff --git a/R/mkinfit.R b/R/mkinfit.R
index f1a3e6bb..2480c135 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -59,6 +59,31 @@ mkinfit <- function(mkinmod, observed,
          " not used in the model")

   }

 

+  # Warn that the sum of formation fractions may exceed one they are not

+  # fitted in the transformed way

+  if (mkinmod$use_of_ff == "max" & transform_fractions == FALSE) {

+    warning("The sum of formation fractions may exceed one if you do not use ",

+            "transform_fractions = TRUE." )

+    for (box in mod_vars) {

+      # Stop if formation fractions are not transformed and we have no sink

+      if (mkinmod$spec[[box]]$sink == FALSE) {

+        stop("If formation fractions are not transformed during the fitting, ",

+             "it is not supported to turn off pathways to sink.\n ",

+             "Consider turning on the transformation of formation fractions or ",

+             "setting up a model with use_of_ff = 'min'.\n")

+      }

+    }

+  }

+

+  # Do not allow fixing formation fractions if we are using the ilr transformation,

+  # this is not supported

+  if (transform_fractions == TRUE && length(fixed_parms) > 0) {

+    if (grepl("^f_", fixed_parms)) {

+      stop("Fixing formation fractions is not supported when using the ilr ",

+           "transformation.")

+    }

+ }

+

   # Set initial parameter values, including a small increment (salt)

   # to avoid linear dependencies (singular matrix) in Eigenvalue based solutions

   k_salt = 0

@@ -72,8 +97,6 @@ mkinfit <- function(mkinmod, observed,
     # 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 formation fractions

-    if (substr(parmname, 1, 2) == "f_") parms.ini[parmname] = 0.2

     # Default values for the FOMC, DFOP and HS models

     if (parmname == "alpha") parms.ini[parmname] = 1

     if (parmname == "beta") parms.ini[parmname] = 10

@@ -82,19 +105,47 @@ mkinfit <- function(mkinmod, observed,
     if (parmname == "tb") parms.ini[parmname] = 5

     if (parmname == "g") parms.ini[parmname] = 0.5

   }

+  # Default values for formation fractions in case they are used

+  if (mkinmod$use_of_ff == "max") {

+    for (box in mod_vars) {

+      f_names <- mkinmod$parms[grep(paste0("^f_", box), mkinmod$parms)]

+      f_default_names <- intersect(f_names, defaultpar.names)

+      f_specified_names <- setdiff(f_names, defaultpar.names)

+      sum_f_specified = sum(parms.ini[f_specified_names])

+      if (sum_f_specified > 1) {

+        stop("Starting values for the formation fractions originating from ",

+             box, " sum up to more than 1.")

+      }

+      if (mkinmod$spec[[box]]$sink) n_unspecified = length(f_default_names) + 1

+      else {

+        n_unspecified = length(f_default_names)

+        # When we have no sink and only one pathway, we get an implicitly

+        # fixed parameter

+        if (length(f_names) == 1) fixed_parms = c(fixed_parms, f_names)

+      }

+      parms.ini[f_default_names] <- (1 - sum_f_specified) / n_unspecified

+    }

+  }

 

   # Name the inital state variable values if they are not named yet

   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, mkinmod,

                                        transform_rates = transform_rates,

                                        transform_fractions = transform_fractions)

 

   # Parameters to be optimised:

   # Kinetic parameters in parms.ini whose names are not in fixed_parms

-  parms.fixed <- transparms.ini[fixed_parms]

-  parms.optim <- transparms.ini[setdiff(names(transparms.ini), fixed_parms)]

+  parms.fixed <- parms.ini[fixed_parms]

+  parms.optim <- parms.ini[setdiff(names(parms.ini), fixed_parms)]

+

+  transparms.fixed <- transform_odeparms(parms.fixed, mkinmod,

+                                       transform_rates = transform_rates,

+                                       transform_fractions = transform_fractions)

+  transparms.optim <- transform_odeparms(parms.optim, mkinmod,

+                                       transform_rates = transform_rates,

+                                       transform_fractions = transform_fractions)

 

   # Inital state variables in state.ini whose names are not in fixed_initials

   state.ini.fixed <- state.ini[fixed_initials]

@@ -160,9 +211,9 @@ mkinfit <- function(mkinmod, observed,
       names(odeini) <- state.ini.fixed.boxnames

     }

 

-    odeparms <- c(P[(length(state.ini.optim) + 1):length(P)], parms.fixed)

+    odeparms <- c(P[(length(state.ini.optim) + 1):length(P)], transparms.fixed)

 

-    parms <- backtransform_odeparms(odeparms, mod_vars,

+    parms <- backtransform_odeparms(odeparms, mkinmod,

                                     transform_rates = transform_rates,

                                     transform_fractions = transform_fractions)

 

@@ -210,9 +261,9 @@ mkinfit <- function(mkinmod, observed,
     return(mC)

   }

 

-  lower <- rep(-Inf, length(c(state.ini.optim, parms.optim)))

-  upper <- rep(Inf, length(c(state.ini.optim, parms.optim)))

-  names(lower) <- names(upper) <- c(names(state.ini.optim), names(parms.optim))

+  lower <- rep(-Inf, length(c(state.ini.optim, transparms.optim)))

+  upper <- rep(Inf, length(c(state.ini.optim, transparms.optim)))

+  names(lower) <- names(upper) <- c(names(state.ini.optim), names(transparms.optim))

   if (!transform_rates) {

     index_k <- grep("^k_", names(lower))

     lower[index_k] <- 0

@@ -229,7 +280,7 @@ mkinfit <- function(mkinmod, observed,
     upper[other_fraction_parms] <- 1

   }

 

-  fit <- modFit(cost, c(state.ini.optim, parms.optim), 

+  fit <- modFit(cost, c(state.ini.optim, transparms.optim), 

                 method = method.modFit, control = control.modFit, 

                 lower = lower, upper = upper, ...)

 

@@ -280,33 +331,32 @@ mkinfit <- function(mkinmod, observed,
   fit$obs_vars <- obs_vars

   fit$predicted <- mkin_wide_to_long(out_predicted, time = "time")

 

-  # Collect initial parameter values in two dataframes

+  # Backtransform parameters

+  bparms.optim = backtransform_odeparms(fit$par, fit$mkinmod,

+                                        transform_rates = transform_rates,

+                                        transform_fractions = transform_fractions)

+  # As backtransform_odeparms does not know about fixed values, it

+  # generates a formation fraction even it is an implicitly fixed one.

+  # This needs to be removed from bparms.optim

+  bparms.optim = bparms.optim[setdiff(names(bparms.optim), names(parms.fixed))]

+  bparms.fixed = c(state.ini.fixed, parms.fixed)

+  bparms.all = c(bparms.optim, parms.fixed)

+

+  # Collect initial parameter values in three dataframes

   fit$start <- data.frame(value = c(state.ini.optim, 

-		  backtransform_odeparms(parms.optim, mod_vars,

-                             transform_rates = transform_rates,

-                             transform_fractions = transform_fractions)))

+                                    parms.optim))

   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$start$lower_bound = lower

-  fit$start$upper_bound = upper

-

-  fit$fixed <- data.frame(value = c(state.ini.fixed, 

-      backtransform_odeparms(parms.fixed, mod_vars,

-                             transform_rates = transform_rates,

-                             transform_fractions = transform_fractions)))

+

+  fit$start_transformed = data.frame(

+      value = c(state.ini.optim, transparms.optim),

+      lower = lower,

+      upper = upper)

+

+  fit$fixed <- data.frame(value = c(state.ini.fixed, parms.fixed))

   fit$fixed$type = c(rep("state", length(state.ini.fixed)), 

                      rep("deparm", length(parms.fixed)))

 

-  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,

-                                        transform_rates = transform_rates,

-                                        transform_fractions = transform_fractions)

-  bparms.all = c(bparms.optim, bparms.fixed)

-

   # Collect observed, predicted, residuals and weighting

   data <- merge(fit$observed, fit$predicted, by = c("time", "name"))

   data$name <- ordered(data$name, levels = obs_vars)

@@ -327,10 +377,13 @@ mkinfit <- function(mkinmod, observed,
   fit$reweight.tol <- reweight.tol

   fit$reweight.max.iter <- reweight.max.iter

 

-  # Return all backtransformed parameters for summary

+  # Return different sets of backtransformed parameters for summary and plotting

   fit$bparms.optim <- bparms.optim 

   fit$bparms.fixed <- bparms.fixed

-  fit$bparms.ode <- bparms.all[mkinmod$parms] # Return ode parameters for further fitting

+

+  # Return ode parameters for further fitting

+  fit$bparms.ode <- bparms.all[mkinmod$parms] 

+

   fit$date <- date()

 

   class(fit) <- c("mkinfit", "modFit") 

@@ -340,6 +393,7 @@ mkinfit <- function(mkinmod, observed,
 summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, ...) {

   param  <- object$par

   pnames <- names(param)

+  bpnames <- names(object$bparms.optim)

   p      <- length(param)

   mod_vars <- names(object$mkinmod$diffs)

   covar  <- try(solve(0.5*object$hessian), silent = TRUE)   # unscaled covariance

@@ -365,7 +419,8 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05,
   dimnames(param) <- list(pnames, c("Estimate", "Std. Error", "Lower", "Upper",

                                     "t value", "Pr(>|t|)", "Pr(>t)"))

 

-  blci <- buci <- numeric()

+  bparam <- cbind(Estimate = object$bparms.optim, Lower = NA, Upper = NA)

+

   # Transform boundaries of CI for one parameter at a time,

   # with the exception of sets of formation fractions (single fractions are OK).

   f_names_skip <- character(0)

@@ -376,20 +431,24 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05,
   }

 

   for (pname in pnames) {

-    par.lower <- par.upper <- object$par

-    par.lower[pname] <- param[pname, "Lower"]

-    par.upper[pname] <- param[pname, "Upper"]

     if (!pname %in% f_names_skip) {

-      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]

-    } else {

-      blci[pname] <- buci[pname] <- NA

-    }

+      par.lower <- param[pname, "Lower"]

+      par.upper <- param[pname, "Upper"]

+      names(par.lower) <- names(par.upper) <- pname

+      bpl <- backtransform_odeparms(par.lower, object$mkinmod, 

+                                            object$transform_rates, 

+                                            object$transform_fractions)

+      bpu <- backtransform_odeparms(par.upper, object$mkinmod,

+                                            object$transform_rates, 

+                                            object$transform_fractions)

+      # Again, we have to remove formation fractions that were implicitly generated

+      bpl <- bpl[intersect(bpnames, names(bpl))]

+      bpu <- bpu[intersect(bpnames, names(bpu))]

+

+      bparam[names(bpl), "Lower"] <- bpl

+      bparam[names(bpu), "Upper"] <- bpu

+    } 

   }

-  bparam <- cbind(object$bparms.optim, blci, buci)

-  dimnames(bparam) <- list(pnames, c("Estimate", "Lower", "Upper"))

 

   ans <- list(

     version = as.character(packageVersion("mkin")),

@@ -416,6 +475,7 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05,
   ans$diffs <- object$mkinmod$diffs

   if(data) ans$data <- object$data

   ans$start <- object$start

+  ans$start_transformed <- object$start_transformed

 

   ans$fixed <- object$fixed

 

@@ -451,9 +511,12 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), .
                                       x$reweight.method)

   cat("\n")

 

-  cat("\nStarting values for optimised parameters:\n")

+  cat("\nStarting values for parameters to be optimised:\n")

   print(x$start)

 

+  cat("\nStarting values for the transformed parameters actually optimised:\n")

+  print(x$start_transformed)

+

   cat("\nFixed parameter values:\n")

   if(length(x$fixed$value) == 0) cat("None\n")

   else print(x$fixed)

@@ -487,8 +550,8 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), .
   }

 

   printff <- !is.null(x$ff)

-  if(printff & x$use_of_ff == "min"){

-    cat("\nEstimated formation fractions:\n")

+  if(printff){

+    cat("\nResulting formation fractions:\n")

     print(data.frame(ff = x$ff), digits=digits,...)

   }

 

diff --git a/R/mkinmod.R b/R/mkinmod.R
index 49ec4a44..fe5e8142 100644
--- a/R/mkinmod.R
+++ b/R/mkinmod.R
@@ -164,9 +164,6 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL)
           diffs[[target_box]] <- paste(diffs[[target_box]], "+", 

             k_from_to, "*", origin_box)

         } else {

-          if (!spec[[varname]]$sink) {

-            stop("Turning off the sink when using formation fractions is not supported")

-          }

           fraction_to_target = paste("f", origin_box, "to", target, sep="_")

           parms <- c(parms, fraction_to_target)

           diffs[[target_box]] <- paste(diffs[[target_box]], "+", 

diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R
index 6932fdef..4774fcf6 100644
--- a/R/transform_odeparms.R
+++ b/R/transform_odeparms.R
@@ -16,31 +16,52 @@
 # You should have received a copy of the GNU General Public License along with

 # this program. If not, see <http://www.gnu.org/licenses/>

 

-transform_odeparms <- function(parms, mod_vars, 

+transform_odeparms <- function(parms, mkinmod,

                                transform_rates = TRUE, 

                                transform_fractions = TRUE) 

 {

+  # We need the model specification for the names of the model 

+  # variables and the information on the sink

+  spec = mkinmod$spec

+

   # Set up container for transformed parameters

-  transparms <- parms

+  transparms <- numeric(0)

+

+  # Do not transform initial values for state variables

+  state.ini.optim <- parms[grep("_0$", names(parms))]

+  transparms[names(state.ini.optim)] <- state.ini.optim

 

   # Log transformation for rate constants if requested

-  index_k <- grep("^k_", names(transparms))

-  if (length(index_k) > 0) {

-    if(transform_rates) transparms[index_k] <- log(parms[index_k])

-    else transparms[index_k] <- parms[index_k]

+  k <- parms[grep("^k_", names(parms))]

+  if (length(k) > 0) {

+    if(transform_rates) {

+      transparms[paste0("log_", names(k))] <- log(k)

+    }

+    else transparms[names(k)] <- k

   }

 

-  # Go through state variables and apply isometric logratio transformation if requested

+  # Go through state variables and apply isometric logratio transformation to

+  # formation fractions if requested

+  mod_vars = names(spec)

   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)

-    n_paths <- length(indices_f)

-    if (n_paths > 0) {

-      f <- parms[indices_f]

-      trans_f <- ilr(c(f, 1 - sum(f)))

-      names(trans_f) <- f_names

-      if(transform_fractions) transparms[indices_f] <- trans_f

-      else transparms[indices_f] <- f

+    f <- parms[grep(paste("^f", box, sep = "_"), names(parms))]

+

+    if (length(f) > 0) {

+      if(transform_fractions) {

+        if (spec[[box]]$sink) {

+          trans_f <- ilr(c(f, 1 - sum(f)))

+          trans_f_names <- paste("f", box, "ilr", 1:length(trans_f), sep = "_")

+          transparms[trans_f_names] <- trans_f

+        } else {

+          if (length(f) > 1) {

+            trans_f <- ilr(f)

+            trans_f_names <- paste("f", box, "ilr", 1:length(trans_f), sep = "_")

+            transparms[trans_f_names] <- trans_f

+          }

+        }

+      } else {

+        transparms[names(f)] <- f

+      }

     }

   }

 

@@ -48,55 +69,85 @@ transform_odeparms <- function(parms, mod_vars,
   # and HS parameter tb if transformation of rates is requested

   for (pname in c("alpha", "beta", "k1", "k2", "tb")) {

     if (!is.na(parms[pname])) {

-      transparms[pname] <- ifelse(transform_rates, log(parms[pname]), parms[pname])

-      transparms[pname] <- ifelse(transform_rates, log(parms[pname]), parms[pname])

+      transparms[paste0("log_", pname)] <- ifelse(transform_rates, log(parms[pname]), parms[pname])

     } 

   }

   if (!is.na(parms["g"])) {

     g <- parms["g"]

-    transparms["g"] <- ifelse(transform_fractions, ilr(c(g, 1 - g)), g)

+    transparms["g_ilr"] <- ifelse(transform_fractions, ilr(c(g, 1 - g)), g)

   }

 

   return(transparms)

 }

 

-backtransform_odeparms <- function(transparms, mod_vars,

+backtransform_odeparms <- function(transparms, mkinmod,

                                    transform_rates = TRUE,

                                    transform_fractions = TRUE)

 {

+  # We need the model specification for the names of the model 

+  # variables and the information on the sink

+  spec = mkinmod$spec

+

   # Set up container for backtransformed parameters

-  parms <- transparms

+  parms <- numeric(0)

+

+  # Do not transform initial values for state variables

+  state.ini.optim <- transparms[grep("_0$", names(transparms))]

+  parms[names(state.ini.optim)] <- state.ini.optim

 

   # Exponential transformation for rate constants

-  index_k <- grep("^k_", names(parms))

-  if (length(index_k) > 0) {

-    if(transform_rates) parms[index_k] <- exp(transparms[index_k])

-    else parms[index_k] <- transparms[index_k]

+  if(transform_rates) {

+    trans_k <- transparms[grep("^log_k_", names(transparms))]

+    if (length(trans_k) > 0) {

+      k_names <- gsub("^log_k_", "k_", names(trans_k))

+      parms[k_names] <- exp(trans_k)

+    }

+  } else {

+    trans_k <- transparms[grep("^k_", names(transparms))]

+    parms[names(trans_k)] <- trans_k

   }

 

   # Go through state variables and apply inverse isometric logratio transformation

+  mod_vars = names(spec)

   for (box in mod_vars) {

-    indices_f <- grep(paste("^f", box, sep = "_"), names(transparms))

-    f_names <- grep(paste("^f", box, sep = "_"), names(transparms), value = TRUE)

-    n_paths <- length(indices_f)

-    if (n_paths > 0) {

-      f <- invilr(transparms[indices_f])[1:n_paths] # We do not need the last component

-      names(f) <- f_names

-      if(transform_fractions) parms[indices_f] <- f

-      else parms[indices_f] <- transparms[indices_f]

+    # Get the names as used in the model

+    f_names = grep(paste("^f", box, sep = "_"), mkinmod$parms, value = TRUE)

+    # Get the formation fraction parameters

+    trans_f = transparms[grep(paste("^f", box, sep = "_"), names(transparms))]

+

+    # If we have one formation fraction parameter, but no optimised parameter,

+    # the one must be unity

+    if (length(trans_f) == 0 & length(f_names == 1)) {

+      parms[f_names] = 1

+    }

+    if (length(trans_f) > 0) {

+      if(transform_fractions) {

+        f <- invilr(trans_f)

+        if (spec[[box]]$sink) {

+          parms[f_names] <- f[1:length(f)-1]

+        } else {

+          parms[f_names] <- f

+        }

+      } else {

+        parms[names(trans_f)] <- trans_f

+      }

     }

   }

 

   # Transform parameters also for FOMC, DFOP and HS models

   for (pname in c("alpha", "beta", "k1", "k2", "tb")) {

-    if (!is.na(transparms[pname])) {

-      parms[pname] <- ifelse(transform_rates, exp(transparms[pname]), transparms[pname])

+    pname_trans = paste0("log_", pname)

+    if (!is.na(transparms[pname_trans])) {

+      parms[pname] <- ifelse(transform_rates, 

+                             exp(transparms[pname_trans]), 

+                             transparms[pname])

     } 

   }

-  if (!is.na(transparms["g"])) {

-    g <- transparms["g"]

-    parms["g"] <- ifelse(transform_fractions, invilr(g)[1], g)

+  if (!is.na(transparms["g_ilr"])) {

+    g_ilr <- transparms["g_ilr"]

+    parms["g"] <- ifelse(transform_fractions, invilr(g_ilr)[1], g_ilr)

   }

 

   return(parms)

 }

+# vim: set ts=2 sw=2 expandtab: