From b5b446b718b15ccaae5b197e147fc1358f0f564e Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Fri, 6 Nov 2020 00:03:29 +0100 Subject: Fast analytical solutions for saemix, update.mmkin Also, use logit transformation for g and for solitary formation fractions, addressing #10. --- man/mkinfit.Rd | 18 +++++++------- man/mmkin.Rd | 2 +- man/saemix.Rd | 61 +++++++++++++++++++++-------------------------- man/transform_odeparms.Rd | 19 ++++++++------- 4 files changed, 47 insertions(+), 53 deletions(-) (limited to 'man') 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{ -- cgit v1.2.1