Support formation fractions without sink pathway, updates

12 files changed, 317 insertions, 197 deletions

diff --git a/GNUmakefile b/GNUmakefile
index 969b4ff0..693ec892 100644
--- a/GNUmakefile
+++ b/GNUmakefile
@@ -26,7 +26,7 @@ pkgfiles = NEWS \
 	   TODO \
 	   vignettes/*
 
-all: NEWS check
+all: NEWS check clean
 
 # convert markdown to R's NEWS format (from knitr package)
 NEWS: NEWS.md
@@ -57,7 +57,12 @@ check: build
 	"$(RBIN)/R" CMD check --as-cran --no-tests --no-build-vignettes $(TGZ)
 
 check-no-vignettes: build-no-vignettes
-	"$(RBIN)/R" CMD check --as-cran --no-tests $(TGZVNR)
+	mv $(TGZVNR) $(TGZ)
+	"$(RBIN)/R" CMD check --as-cran --no-tests $(TGZ)
+	mv $(TGZ) $(TGZVNR)
+
+clean: 
+	$(RM) -r $(PKGNAME).Rcheck/
 
 test: install-no-vignettes
 	cd tests;\
@@ -74,9 +79,6 @@ move-sd:
 	rm -rf $(SDDIR)/*;\
 		cp -r inst/web/* $(SDDIR)
 
-#------------------------------------------------------------------------------
-# Packaging Tasks
-#------------------------------------------------------------------------------
 winbuilder: build
 	date
 	@echo "Uploading to R-release on win-builder"
diff --git a/NEWS.md b/NEWS.md
index 8d65b41f..6c23f316 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,13 +1,29 @@
 # CHANGES in mkin VERSION 0.9-30 
 
-- The ChangeLog was renamed to NEWS.md and the entries were converted to markdown syntax compatible with the `tools::news()` function built into R.
+## NEW FEATURES
+
+- It is now possible to use formation fractions in combination with turning off the sink in `mkinmod()`.
+
+## MAJOR CHANGES
+
+- The original and the transformed parameters now have different names (e.g. `k_parent` and `log_k_parent`. They also differ in how many they are when we have formation fractions but no pathway to sink.
 
 - The order of some of the information blocks in `print.summary.mkinfit.R()` has been ordered in a more logical way
 
-- The vignette FOCUS_Z has been slightly amended to use the new versions of DT50 values calculated since mkin 0.9-29.
+## MINOR CHANGES
+
+- The vignette FOCUS_Z has been simplified to use formation fractions with turning off the sink, and slightly amended to use the new versions of DT50 values calculated since mkin 0.9-29.
 
 - All vignettes have been rebuilt so they reflect all changes
 
+- The ChangeLog was renamed to NEWS.md and the entries were converted to markdown syntax compatible with the `tools::news()` function built into R.
+
+- The test suite was overhauled. Tests of the DFOP and SFORB models with dataset FOCUS_2006_A were removed, as they were too much dependent on the optimisation algorithm and/or starting parameters, because the dataset is SFO (compare kinfit vignette).
+
+- Also, the Schaefer complex case can now be fitted using formation fractions, and with the 'Port' optimisation method we also fit A2 in the same way as published in the Piacenza paper.
+
+- Some more checks were introduced to `mkinfit()`, leading to warnings or stopping execution if unsupported combinations of methods and parameters are requested.
+
 # CHANGES in mkin VERSION 0.9-29
 
 - R/mkinresplot.R: Make it possible to specify `xlim`
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:

diff --git a/TODO b/TODO
index 27d84f69..af4ff82d 100644
--- a/TODO
+++ b/TODO
@@ -1,11 +1,7 @@
 TODO for version 1.0
 - Think about what a user would expect from version 1.0
-- Support model definitions without pathway to sink in combination with
-  formation fractions
 - Complete the main package vignette named mkin to include a method description
 - Improve formatting of differential equations in the summary
-- Rename formation fractions during transformation and backtransformation as there is no 
-  one to one relationship
 
 Nice to have:
 - Calculate confidence intervals for DT50 and DT90 values when only one
diff --git a/inst/unitTests/runit.mkinfit.R b/inst/unitTests/runit.mkinfit.R
index 61a0b303..fdbc86e0 100644
--- a/inst/unitTests/runit.mkinfit.R
+++ b/inst/unitTests/runit.mkinfit.R
@@ -100,15 +100,9 @@ test.FOCUS_2006_DFOP <- function()
   DFOP <- mkinmod(parent = list(type = "DFOP"))

 

   # FOCUS_2006_A

-  fit.A.DFOP <- mkinfit(DFOP, FOCUS_2006_A, quiet=TRUE)

-

-  median.A.DFOP <- as.numeric(lapply(subset(FOCUS_2006_DFOP_ref_A_to_B, dataset == "A", 

-                   c(M0, k1, k2, f, DT50, DT90)), "median"))

-

-  fit.A.DFOP.r <- as.numeric(c(fit.A.DFOP$bparms.optim, endpoints(fit.A.DFOP)$distimes[c("DT50", "DT90")]))

-  dev.A.DFOP <- abs(round(100 * ((median.A.DFOP - fit.A.DFOP.r)/median.A.DFOP), digits=1))

-  # about 6.7% deviation for parameter f, the others are < 0.1%

-  checkIdentical(dev.A.DFOP < c(1, 1, 1, 10, 1, 1), rep(TRUE, length(dev.A.DFOP)))

+  # Results were too much dependent on algorithm, as this dataset

+  # is pretty much SFO. "Port" gave a lower deviance, but deviated from the

+  # median of FOCUS_2006 solutions

 

   # FOCUS_2006_B

   fit.B.DFOP <- mkinfit(DFOP, FOCUS_2006_B, quiet=TRUE)

@@ -157,10 +151,10 @@ test.FOCUS_2006_HS <- function()
   median.C.HS <- as.numeric(lapply(subset(FOCUS_2006_HS_ref_A_to_F, dataset == "C", 

                    c(M0, k1, k2, tb, DT50, DT90)), "median"))

 

-  fit.A.HS.r <- as.numeric(c(fit.A.HS$bparms.optim, endpoints(fit.A.HS)$distimes[c("DT50", "DT90")]))

-  dev.A.HS <- abs(round(100 * ((median.A.HS - fit.A.HS.r)/median.A.HS), digits=1))

+  fit.C.HS.r <- as.numeric(c(fit.C.HS$bparms.optim, endpoints(fit.C.HS)$distimes[c("DT50", "DT90")]))

+  dev.C.HS <- abs(round(100 * ((median.C.HS - fit.C.HS.r)/median.C.HS), digits=1))

   # deviation <= 0.1%

-  checkIdentical(dev.A.HS < 1, rep(TRUE, length(dev.A.HS)))

+  checkIdentical(dev.C.HS < 1, rep(TRUE, length(dev.C.HS)))

 } # }}}

 

 # Test SFORB model against DFOP solutions to a relative tolerance of 1% # {{{

@@ -169,33 +163,7 @@ test.FOCUS_2006_SFORB <- function()
   SFORB <- mkinmod(parent = list(type = "SFORB"))

 

   # FOCUS_2006_A

-  fit.A.SFORB.1 <- mkinfit(SFORB, FOCUS_2006_A, quiet=TRUE)

-  fit.A.SFORB.2 <- mkinfit(SFORB, FOCUS_2006_A, solution_type = "deSolve", quiet=TRUE)

-

-  median.A.SFORB <- as.numeric(lapply(subset(FOCUS_2006_DFOP_ref_A_to_B, dataset == "A", 

-                   c(M0, k1, k2, DT50, DT90)), "median"))

-

-  fit.A.SFORB.1.r <- as.numeric(c(

-                      parent_0 = fit.A.SFORB.1$bparms.optim[[1]], 

-                      k1 = endpoints(fit.A.SFORB.1)$SFORB[[1]],

-                      k2 = endpoints(fit.A.SFORB.1)$SFORB[[2]],

-                      endpoints(fit.A.SFORB.1)$distimes[c("DT50", "DT90")]))

-  dev.A.SFORB.1 <- abs(round(100 * ((median.A.SFORB - fit.A.SFORB.1.r)/median.A.SFORB), digits=1))

-  # The first Eigenvalue is a lot different from k1 in the DFOP fit

-  # The explanation is that the dataset is simply SFO

-  dev.A.SFORB.1 <- dev.A.SFORB.1[c(1, 3, 4, 5)]

-  checkIdentical(dev.A.SFORB.1 < 1, rep(TRUE, length(dev.A.SFORB.1)))

-

-  fit.A.SFORB.2.r <- as.numeric(c(

-                      parent_0 = fit.A.SFORB.2$bparms.optim[[1]], 

-                      k1 = endpoints(fit.A.SFORB.2)$SFORB[[1]],

-                      k2 = endpoints(fit.A.SFORB.2)$SFORB[[2]],

-                      endpoints(fit.A.SFORB.2)$distimes[c("DT50", "DT90")]))

-  dev.A.SFORB.2 <- abs(round(100 * ((median.A.SFORB - fit.A.SFORB.2.r)/median.A.SFORB), digits=1))

-  # The first Eigenvalue is a lot different from k1 in the DFOP fit

-  # The explanation is that the dataset is simply SFO

-  dev.A.SFORB.2 <- dev.A.SFORB.2[c(1, 3, 4, 5)]

-  checkIdentical(dev.A.SFORB.2 < 1, rep(TRUE, length(dev.A.SFORB.2)))

+  # Again it does not make a lot of sense to use a SFO dataset for this

 

   # FOCUS_2006_B

   fit.B.SFORB.1 <- mkinfit(SFORB, FOCUS_2006_B, quiet=TRUE)

@@ -263,32 +231,33 @@ test.mkinfit.schaefer07_complex_example <- function()
     A1 = list(type = "SFO", to = "A2"),

     B1 = list(type = "SFO"),

     C1 = list(type = "SFO"),

-    A2 = list(type = "SFO"))

+    A2 = list(type = "SFO"), use_of_ff = "max")

   

-  fit <- mkinfit(schaefer07_complex_model, 

+  # If we use the default algorithm 'Marq' we need to give a good starting 

+  # estimate for k_A2 in order to find the solution published by Schaefer et al.

+  fit <- mkinfit(schaefer07_complex_model, method.modFit = "Port",

     mkin_wide_to_long(schaefer07_complex_case, time = "time"))

   s <- summary(fit)

   r <- schaefer07_complex_results

   attach(as.list(fit$bparms.optim))

-  k_parent <- sum(k_parent_A1, k_parent_B1, k_parent_C1)

   r$mkin <- c(

     k_parent,

     s$distimes["parent", "DT50"],

     s$ff["parent_A1"],

-    sum(k_A1_sink, k_A1_A2),

+    k_A1,

     s$distimes["A1", "DT50"],

     s$ff["parent_B1"],

-    k_B1_sink,

+    k_B1,

     s$distimes["B1", "DT50"],

     s$ff["parent_C1"],

-    k_C1_sink,

+    k_C1,

     s$distimes["C1", "DT50"],

     s$ff["A1_A2"],

-    k_A2_sink,

+    k_A2,

     s$distimes["A2", "DT50"])

   r$means <- (r$KinGUI + r$ModelMaker)/2

   r$mkin.deviation <- abs(round(100 * ((r$mkin - r$means)/r$means), digits=1))

-  checkIdentical(r$mkin.deviation[1:11] < 10, rep(TRUE, 11))

+  checkIdentical(r$mkin.deviation < 10, rep(TRUE, 14))

 } # }}}

 

 # vim: set foldmethod=marker ts=2 sw=2 expandtab:

diff --git a/man/mkinfit.Rd b/man/mkinfit.Rd
index 823ceeec..bd7f73b7 100644
--- a/man/mkinfit.Rd
+++ b/man/mkinfit.Rd
@@ -94,11 +94,16 @@ mkinfit(mkinmod, observed,
     "lsoda" is performant, but sometimes fails to converge.
   }
   \item{method.modFit}{
-    The optimisation method passed to \code{\link{modFit}}. The default "Marq"
+    The optimisation method passed to \code{\link{modFit}}.  The default "Marq"
     is the Levenberg Marquardt algorithm \code{\link{nls.lm}} from the package
-    \code{minpack.lm}. Often other methods need more iterations to find the
-    same result. When using "Pseudo", "upper" and "lower" need to be 
-    specified for the transformed parameters.
+    \code{minpack.lm} and usually needs the least number of iterations.
+
+    For more complex problems where local minima occur, the "Port" algorithm is
+    recommended as it is less prone to get trapped in local minima and depends
+    less on starting values for parameters.  However, it needs more iterations.
+
+    When using "Pseudo", "upper" and "lower" need to be specified for the
+    transformed parameters.
   }
   \item{control.modFit}{
     Additional arguments passed to the optimisation method used by
diff --git a/man/mkinparplot.Rd b/man/mkinparplot.Rd
index 60103b89..f40a96d1 100644
--- a/man/mkinparplot.Rd
+++ b/man/mkinparplot.Rd
@@ -25,11 +25,10 @@
 

 \examples{

 model <- mkinmod(

-  T245 = list(type = "SFO", to = c("phenol")),

+  T245 = list(type = "SFO", to = c("phenol"), sink = FALSE),

   phenol = list(type = "SFO", to = c("anisole")),

   anisole = list(type = "SFO"), use_of_ff = "max")

-fit <- mkinfit(model, subset(mccall81_245T, soil == "Commerce"), 

-  parms.ini = c(f_phenol_to_anisole = 1), fixed_parms = "f_phenol_to_anisole")

+fit <- mkinfit(model, subset(mccall81_245T, soil == "Commerce"))

 \dontrun{mkinparplot(fit)}

 }

 \keyword{ hplot }

diff --git a/man/transform_odeparms.Rd b/man/transform_odeparms.Rd
index c52fb4f1..ea0b5024 100644
--- a/man/transform_odeparms.Rd
+++ b/man/transform_odeparms.Rd
@@ -11,10 +11,15 @@
   simple log transformation is used. For compositional parameters, such as
   the formations fractions that should always sum up to 1 and can not be
   negative, the \code{\link{ilr}} transformation is used.
+
+  The transformation of sets of formation fractions is fragile, as it supposes
+  the same ordering of the components in forward and backward transformation.
+  This is no problem for the internal use in \code{\link{mkinfit}}.
 }
 \usage{
-transform_odeparms(parms, mod_vars, transform_rates = TRUE, transform_fractions = TRUE)
-backtransform_odeparms(transparms, mod_vars, 
+transform_odeparms(parms, mkinmod, 
+                   transform_rates = TRUE, transform_fractions = TRUE)
+backtransform_odeparms(transparms, mkinmod, 
                        transform_rates = TRUE, transform_fractions = TRUE)
 }
 \arguments{
@@ -24,9 +29,11 @@ backtransform_odeparms(transparms, mod_vars,
   \item{transparms}{
     Transformed parameters of kinetic models as used in the fitting procedure.
   }
-  \item{mod_vars}{
-    Names of the state variables in the kinetic model. These are used for 
-    grouping the formation fractions before \code{\link{ilr}} transformation.  
+  \item{mkinmod}{
+    The kinetic model of class \code{\link{mkinmod}}, containing the names of
+    the model variables that are needed for grouping the formation fractions
+    before \code{\link{ilr}} transformation, the parameter names and 
+    the information if the pathway to sink is included in the model. 
   }
   \item{transform_rates}{
     Boolean specifying if kinetic rate constants should be transformed in the
@@ -59,11 +66,36 @@ SFO_SFO <- mkinmod(
   m1 = list(type = "SFO"))
 # Fit the model to the FOCUS example dataset D using defaults
 fit <- mkinfit(SFO_SFO, FOCUS_2006_D)
+fit.2 <- mkinfit(SFO_SFO, FOCUS_2006_D, transform_rates = FALSE)
 summary(fit, data=FALSE) # See transformed and backtransformed parameters
+summary(fit.2, data=FALSE) 
 initials <- fit$start$value
-transformed <- fit$start$transformed
-names(initials) <- names(transformed) <- rownames(fit$start)
-transform_odeparms(initials, c("parent", "m1"))
-backtransform_odeparms(transformed, c("parent", "m1"))
+names(initials) <- rownames(fit$start)
+transformed <- fit$start_transformed$value
+names(transformed) <- rownames(fit$start_transformed)
+transform_odeparms(initials, SFO_SFO)
+backtransform_odeparms(transformed, SFO_SFO)
+
+# The case of formation fractions
+SFO_SFO.ff <- mkinmod(
+  parent = list(type = "SFO", to = "m1", sink = TRUE),
+  m1 = list(type = "SFO"),
+  use_of_ff = "max")
+
+fit.ff <- mkinfit(SFO_SFO.ff, FOCUS_2006_D)
+summary(fit.ff, data = FALSE)
+initials <- c("f_parent_to_m1" = 0.5)
+transformed <- transform_odeparms(initials, SFO_SFO.ff)
+backtransform_odeparms(transformed, SFO_SFO.ff)
+
+# And without sink
+SFO_SFO.ff.2 <- mkinmod(
+  parent = list(type = "SFO", to = "m1", sink = FALSE),
+  m1 = list(type = "SFO"),
+  use_of_ff = "max")
+
+
+fit.ff.2 <- mkinfit(SFO_SFO.ff.2, FOCUS_2006_D)
+summary(fit.ff.2, data = FALSE)
 }
 \keyword{ manip }
diff --git a/vignettes/FOCUS_Z.Rnw b/vignettes/FOCUS_Z.Rnw
index 40ddfc44..23b229ab 100644
--- a/vignettes/FOCUS_Z.Rnw
+++ b/vignettes/FOCUS_Z.Rnw
@@ -108,24 +108,14 @@ The simplified model is obtained by setting the list component \texttt{sink} to
 
 <<FOCUS_2006_Z_fits_3, echo=TRUE, fig.height=4>>=
 Z.3 <- mkinmod(Z0 = list(type = "SFO", to = "Z1", sink = FALSE),
-               Z1 = list(type = "SFO"))
+               Z1 = list(type = "SFO"), use_of_ff = "max")
 m.Z.3 <- mkinfit(Z.3, FOCUS_2006_Z_mkin, quiet = TRUE)
 plot(m.Z.3)
 summary(m.Z.3, data = FALSE)
 @
 
-This model definition is not supported when formation fractions 
-are used, but the formation fraction can be fixed to unity.
-
-<<FOCUS_2006_Z_fits_3.ff, echo=TRUE, fig.height=4>>=
-Z.3.ff <- mkinmod(Z0 = list(type = "SFO", to = "Z1"),
-               Z1 = list(type = "SFO"), use_of_ff = "max")
-m.Z.3.ff <- mkinfit(Z.3.ff, FOCUS_2006_Z_mkin, 
-                    parms.ini = c(f_Z0_to_Z1 = 1),
-                    fixed_parms = "f_Z0_to_Z1",
-                    quiet = TRUE)
-summary(m.Z.3.ff, data = FALSE)
-@
+As there is only one transformation product for Z0 and no pathway
+to sink, the formation fraction is internally fixed to unity.
 
 \section{Including metabolites Z2 and Z3}
 
@@ -172,17 +162,15 @@ mkinresplot(m.Z.FOCUS, "Z3", lpos = "bottomright")
 
 We can also investigate the confidence interval for the formation
 fraction from Z1 to Z2 by specifying the model using formation
-fractions, and fixing only the formation fraction from Z0 to Z1
-to unity.
+fractions.
 
 <<FOCUS_2006_Z_fits_6_ff, echo=TRUE, fig.height=4>>=
-Z.FOCUS.ff <- mkinmod(Z0 = list(type = "SFO", to = "Z1"),
-                   Z1 = list(type = "SFO", to = "Z2"),
+Z.FOCUS.ff <- mkinmod(Z0 = list(type = "SFO", to = "Z1", sink = FALSE),
+                   Z1 = list(type = "SFO", to = "Z2", sink = FALSE),
                    Z2 = list(type = "SFO", to = "Z3"),
-                   Z3 = list(type = "SFO"), use_of_ff = "max")
-m.Z.FOCUS.ff <- mkinfit(Z.FOCUS.ff, FOCUS_2006_Z_mkin,
-                        parms.ini = c(f_Z0_to_Z1 = 1),
-                        fixed_parms = c("f_Z0_to_Z1"), quiet = TRUE) 
+                   Z3 = list(type = "SFO"), 
+                   use_of_ff = "max")
+m.Z.FOCUS.ff <- mkinfit(Z.FOCUS.ff, FOCUS_2006_Z_mkin, quiet = TRUE)
 plot(m.Z.FOCUS.ff)
 summary(m.Z.FOCUS.ff, data = FALSE)
 @