Fast analytical solutions for saemix, update.mmkin

Also, use logit transformation for g and for solitary formation
fractions, addressing #10.

4 files changed, 47 insertions, 53 deletions

diff --git a/man/mkinfit.Rd b/man/mkinfit.Rd
index 8f10ea0a..768b85d3 100644
--- a/man/mkinfit.Rd
+++ b/man/mkinfit.Rd
@@ -119,12 +119,11 @@ models and the break point tb of the HS model.  If FALSE, zero is used as
 a lower bound for the rates in the optimisation.}
 
 \item{transform_fractions}{Boolean specifying if formation fractions
-constants should be transformed in the model specification used in the
-fitting for better compliance with the assumption of normal distribution
-of the estimator. The default (TRUE) is to do transformations. If TRUE,
-the g parameter of the DFOP and HS models are also transformed, as they
-can also be seen as compositional data. The transformation used for these
-transformations is the \code{\link[=ilr]{ilr()}} transformation.}
+should be transformed in the model specification used in the fitting for
+better compliance with the assumption of normal distribution of the
+estimator. The default (TRUE) is to do transformations. If TRUE,
+the g parameter of the DFOP model is also transformed. Transformations
+are described in \link{transform_odeparms}.}
 
 \item{quiet}{Suppress printing out the current value of the negative
 log-likelihood after each improvement?}
@@ -233,15 +232,14 @@ SFO_SFO <- mkinmod(
 
 # Fit the model quietly to the FOCUS example dataset D using defaults
 fit <- mkinfit(SFO_SFO, FOCUS_D, quiet = TRUE)
-# Since mkin 0.9.50.3, we get a warning about non-normality of residuals,
-# so we try an alternative error model
+plot_sep(fit)
+# As lower parent values appear to have lower variance, we try an alternative error model
 fit.tc <- mkinfit(SFO_SFO, FOCUS_D, quiet = TRUE, error_model = "tc")
 # This avoids the warning, and the likelihood ratio test confirms it is preferable
 lrtest(fit.tc, fit)
 # We can also allow for different variances of parent and metabolite as error model
 fit.obs <- mkinfit(SFO_SFO, FOCUS_D, quiet = TRUE, error_model = "obs")
-# This also avoids the warning about non-normality, but the two-component error model
-# has significantly higher likelihood
+# The two-component error model has significantly higher likelihood
 lrtest(fit.obs, fit.tc)
 parms(fit.tc)
 endpoints(fit.tc)
diff --git a/man/mmkin.Rd b/man/mmkin.Rd
index e06f6cd1..80360cd7 100644
--- a/man/mmkin.Rd
+++ b/man/mmkin.Rd
@@ -8,7 +8,7 @@ more datasets}
 mmkin(
   models = c("SFO", "FOMC", "DFOP"),
   datasets,
-  cores = detectCores(),
+  cores = parallel::detectCores(),
   cluster = NULL,
   ...
 )
diff --git a/man/saemix.Rd b/man/saemix.Rd
index 962526ba..d4a8d0a4 100644
--- a/man/saemix.Rd
+++ b/man/saemix.Rd
@@ -38,46 +38,39 @@ mmkin. Starting variances of the random effects (argument omega.init) are the
 variances of the deviations of the parameters from these mean values.
 }
 \examples{
+\dontrun{
+library(saemix)
 ds <- lapply(experimental_data_for_UBA_2019[6:10],
  function(x) subset(x$data[c("name", "time", "value")]))
 names(ds) <- paste("Dataset", 6:10)
-sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"),
+f_mmkin_parent_p0_fixed <- mmkin("FOMC", ds, cores = 1,
+  state.ini = c(parent = 100), fixed_initials = "parent", quiet = TRUE)
+m_saemix_p0_fixed <- saemix_model(f_mmkin_parent_p0_fixed["FOMC", ])
+d_saemix_parent <- saemix_data(f_mmkin_parent_p0_fixed)
+saemix_options <- list(seed = 123456, displayProgress = FALSE,
+  save = FALSE, save.graphs = FALSE, nbiter.saemix = c(200, 80))
+f_saemix_p0_fixed <- saemix(m_saemix_p0_fixed, d_saemix_parent, saemix_options)
+
+f_mmkin_parent <- mmkin(c("SFO", "FOMC", "DFOP"), ds, quiet = TRUE)
+m_saemix_sfo <- saemix_model(f_mmkin_parent["SFO", ])
+m_saemix_fomc <- saemix_model(f_mmkin_parent["FOMC", ])
+m_saemix_dfop <- saemix_model(f_mmkin_parent["DFOP", ])
+d_saemix_parent <- saemix_data(f_mmkin_parent["SFO", ])
+f_saemix_sfo <- saemix(m_saemix_sfo, d_saemix_parent, saemix_options)
+f_saemix_fomc <- saemix(m_saemix_fomc, d_saemix_parent, saemix_options)
+f_saemix_dfop <- saemix(m_saemix_dfop, d_saemix_parent, saemix_options)
+compare.saemix(list(f_saemix_sfo, f_saemix_fomc, f_saemix_dfop))
+f_mmkin_parent_tc <- update(f_mmkin_parent, error_model = "tc")
+m_saemix_fomc_tc <- saemix_model(f_mmkin_parent_tc["FOMC", ])
+f_saemix_fomc_tc <- saemix(m_saemix_fomc_tc, d_saemix_parent, saemix_options)
+compare.saemix(list(f_saemix_fomc, f_saemix_fomc_tc))
+
+dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"),
   A1 = mkinsub("SFO"))
-\dontrun{
-f_mmkin <- mmkin(list("SFO-SFO" = sfo_sfo), ds, quiet = TRUE)
-library(saemix)
-m_saemix <- saemix_model(f_mmkin, cores = 1)
+f_mmkin <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE)
+m_saemix <- saemix_model(f_mmkin)
 d_saemix <- saemix_data(f_mmkin)
-saemix_options <- list(seed = 123456,
-  save = FALSE, save.graphs = FALSE, displayProgress = FALSE,
-  nbiter.saemix = c(200, 80))
 f_saemix <- saemix(m_saemix, d_saemix, saemix_options)
-plot(f_saemix, plot.type = "convergence")
-}
-# Synthetic data with two-component error
-sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120)
-dt50_sfo_in <- c(80, 90, 100, 111.111, 125)
-k_in <- log(2) / dt50_sfo_in
-
-SFO <- mkinmod(parent = mkinsub("SFO"))
 
-pred_sfo <- function(k) {
-  mkinpredict(SFO, c(k_parent = k),
-    c(parent = 100), sampling_times)
-}
-
-ds_sfo_mean <- lapply(k_in, pred_sfo)
-set.seed(123456L)
-ds_sfo_syn <- lapply(ds_sfo_mean, function(ds) {
-  add_err(ds, sdfunc = function(value) sqrt(1^2 + value^2 * 0.07^2),
-     n = 1)[[1]]
- })
-\dontrun{
-f_mmkin_syn <- mmkin("SFO", ds_sfo_syn, error_model = "tc", quiet = TRUE)
-# plot(f_mmkin_syn)
-m_saemix_tc <- saemix_model(f_mmkin_syn, cores = 1)
-d_saemix_tc <- saemix_data(f_mmkin_syn)
-f_saemix_tc <- saemix(m_saemix_tc, d_saemix_tc, saemix_options)
-plot(f_saemix_tc, plot.type = "convergence")
 }
 }
diff --git a/man/transform_odeparms.Rd b/man/transform_odeparms.Rd
index 710b0875..f38bb051 100644
--- a/man/transform_odeparms.Rd
+++ b/man/transform_odeparms.Rd
@@ -23,9 +23,9 @@ backtransform_odeparms(
 \item{parms}{Parameters of kinetic models as used in the differential
 equations.}
 
-\item{mkinmod}{The kinetic model of class \code{\link{mkinmod}}, containing
+\item{mkinmod}{The kinetic model of class \link{mkinmod}, containing
 the names of the model variables that are needed for grouping the
-formation fractions before \code{\link{ilr}} transformation, the parameter
+formation fractions before \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
@@ -38,10 +38,13 @@ models and the break point tb of the HS model.}
 \item{transform_fractions}{Boolean specifying if formation fractions
 constants should be transformed in the model specification used in the
 fitting for better compliance with the assumption of normal distribution
-of the estimator. The default (TRUE) is to do transformations. The g
-parameter of the DFOP and HS models are also transformed, as they can also
-be seen as compositional data. The transformation used for these
-transformations is the \code{\link{ilr}} transformation.}
+of the estimator. The default (TRUE) is to do transformations.
+The g parameter of the DFOP model is also seen as a fraction.
+If a single fraction is transformed (g parameter of DFOP or only a single
+target variable e.g. a single metabolite plus a pathway to sink), a
+logistic transformation is used \code{\link[stats:Logistic]{stats::qlogis()}}. In other cases, i.e. if
+two or more formation fractions need to be transformed whose sum cannot
+exceed one, the \link{ilr} transformation is used.}
 
 \item{transparms}{Transformed parameters of kinetic models as used in the
 fitting procedure.}
@@ -55,12 +58,12 @@ restricted values to the full scale of real numbers. For kinetic rate
 constants and other parameters that can only take on positive values, a
 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 \link{ilr} transformation is used.
 }
 \details{
 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}}.
+This is no problem for the internal use in \link{mkinfit}.
 }
 \examples{