From a5874ab7fce4616e80be69366ff0685332f47bf1 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Tue, 27 Oct 2020 15:34:14 +0100 Subject: Add summary method for nlme.mmkin objects Improve and update docs --- NAMESPACE | 3 + NEWS.md | 4 +- R/nlme.mmkin.R | 46 ++- R/plot.nlme.mmkin.R | 6 +- R/sigma_twocomp.R | 21 +- R/summary.nlme.mmkin.R | 233 +++++++++++++++ _pkgdown.yml | 1 + docs/dev/news/index.html | 3 +- docs/dev/pkgdown.yml | 2 +- docs/dev/reference/Rplot001.png | Bin 27839 -> 14324 bytes docs/dev/reference/index.html | 8 +- docs/dev/reference/nlme.mmkin.html | 96 ++++--- docs/dev/reference/plot.nlme.mmkin.html | 14 +- docs/dev/reference/sigma_twocomp.html | 11 +- docs/dev/reference/summary.nlme.mmkin.html | 446 +++++++++++++++++++++++++++++ docs/dev/sitemap.xml | 3 + man/nlme.mmkin.Rd | 5 + man/plot.nlme.mmkin.Rd | 8 +- man/sigma_twocomp.Rd | 11 +- man/summary.nlme.mmkin.Rd | 99 +++++++ test.log | 20 +- tests/testthat/FOCUS_2006_D.csf | 2 +- 22 files changed, 949 insertions(+), 93 deletions(-) create mode 100644 R/summary.nlme.mmkin.R create mode 100644 docs/dev/reference/summary.nlme.mmkin.html create mode 100644 man/summary.nlme.mmkin.Rd diff --git a/NAMESPACE b/NAMESPACE index 645c145e..9554a0d6 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -26,8 +26,10 @@ S3method(print,mmkin) S3method(print,nafta) S3method(print,nlme.mmkin) S3method(print,summary.mkinfit) +S3method(print,summary.nlme.mmkin) S3method(residuals,mkinfit) S3method(summary,mkinfit) +S3method(summary,nlme.mmkin) S3method(update,mkinfit) S3method(update,nlme.mmkin) export(CAKE_export) @@ -108,6 +110,7 @@ importFrom(stats,na.fail) importFrom(stats,nlminb) importFrom(stats,nobs) importFrom(stats,optimize) +importFrom(stats,predict) importFrom(stats,pt) importFrom(stats,qchisq) importFrom(stats,qf) diff --git a/NEWS.md b/NEWS.md index ea4b5318..706e0ae8 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,8 +1,6 @@ # mkin 0.9.50.4 (unreleased) -- 'plot' method for 'nlme.mmkin' objects - -- 'print' method for 'mmkin' objects +- 'plot', 'summary' and 'print' methods for 'nlme.mmkin' objects - 'saemix_model', 'saemix_data': Helper functions to fit nonlinear mixed-effects models for mmkin row objects using the saemix package diff --git a/R/nlme.mmkin.R b/R/nlme.mmkin.R index d3369cf5..6d24a044 100644 --- a/R/nlme.mmkin.R +++ b/R/nlme.mmkin.R @@ -42,6 +42,9 @@ get_deg_func <- function() { #' @importFrom stats na.fail as.formula #' @return Upon success, a fitted nlme.mmkin object, which is an nlme object #' with additional elements +#' @note As the object inherits from [nlme::nlme], there is a wealth of +#' methods that will automatically work on 'nlme.mmkin' objects, such as +#' [nlme::intervals()], [nlme::anova.lme()] and [nlme::coef.lme()]. #' @export #' @seealso \code{\link{nlme_function}} #' @examples @@ -141,8 +144,8 @@ nlme.mmkin <- function(model, data = sys.frame(sys.parent()), thisCall[["model"]] <- this_model - mean_dp <- mean_degparms(model) - dp_names <- names(mean_dp) + mean_dp_start <- mean_degparms(model) + dp_names <- names(mean_dp_start) thisCall[["data"]] <- nlme_data(model) @@ -175,10 +178,21 @@ nlme.mmkin <- function(model, data = sys.frame(sys.parent()), thisCall[["control"]] <- control } - val <- do.call("nlme.formula", thisCall) + fit_time <- system.time(val <- do.call("nlme.formula", thisCall)) + val$time <- fit_time + + val$mean_dp_start <- mean_dp_start val$mmkin_orig <- model val$data <- thisCall[["data"]] val$mkinmod <- model[[1]]$mkinmod + val$err_mode <- error_model + val$transform_rates <- model[[1]]$transform_rates + val$transform_fractions <- model[[1]]$transform_fractions + val$solution_type <- model[[1]]$solution_type + val$date.fit <- date() + val$nlmeversion <- as.character(utils::packageVersion("nlme")) + val$mkinversion <- as.character(utils::packageVersion("mkin")) + val$Rversion <- paste(R.version$major, R.version$minor, sep=".") class(val) <- c("nlme.mmkin", "nlme", "lme") return(val) } @@ -186,10 +200,30 @@ nlme.mmkin <- function(model, data = sys.frame(sys.parent()), #' @export #' @rdname nlme.mmkin #' @param x An nlme.mmkin object to print -#' @param ... Further arguments as in the generic print.nlme.mmkin <- function(x, ...) { - x$call$data <- "Not shown" - NextMethod("print", x) + cat( "Kinetic nonlinear mixed-effects model fit by " ) + cat( if(x$method == "REML") "REML\n" else "maximum likelihood\n") + cat("\nStructural model:\n") + diffs <- x$mmkin_orig[[1]]$mkinmod$diffs + nice_diffs <- gsub("^(d.*) =", "\\1/dt =", diffs) + writeLines(strwrap(nice_diffs, exdent = 11)) + cat("\nData:\n") + cat(nrow(x$data), "observations of", + length(unique(x$data$name)), "variable(s) grouped in", + length(unique(x$data$ds)), "datasets\n") + cat("\nLog-", if(x$method == "REML") "restricted-" else "", + "likelihood: ", format(x$logLik), "\n", sep = "") + fixF <- x$call$fixed + cat("\nFixed effects:\n", + deparse( + if(inherits(fixF, "formula") || is.call(fixF) || is.name(fixF)) + x$call$fixed + else + lapply(fixF, function(el) as.name(deparse(el)))), "\n") + print(fixef(x), ...) + cat("\n") + print(summary(x$modelStruct), sigma = x$sigma, ...) + invisible(x) } #' @export diff --git a/R/plot.nlme.mmkin.R b/R/plot.nlme.mmkin.R index afb682a7..05a17a22 100644 --- a/R/plot.nlme.mmkin.R +++ b/R/plot.nlme.mmkin.R @@ -11,6 +11,8 @@ if(getRversion() >= '2.15.1') utils::globalVariables("ds") #' @param rel.height.legend The relative height of the legend shown on top #' @param rel.height.bottom The relative height of the bottom plot row #' @param ymax Vector of maximum y axis values +#' @param ncol.legend Number of columns to use in the legend +#' @param nrow.legend Number of rows to use in the legend #' @param \dots Further arguments passed to \code{\link{plot.mkinfit}} and #' \code{\link{mkinresplot}}. #' @param resplot Should the residuals plotted against time or against @@ -28,7 +30,8 @@ if(getRversion() >= '2.15.1') utils::globalVariables("ds") #' names(ds) <- paste0("ds ", 6:10) #' dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"), #' A1 = mkinsub("SFO"), quiet = TRUE) -#' f <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE, cores = 1) +#' \dontrun{ +#' f <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE) #' plot(f[, 3:4], standardized = TRUE) #' #' library(nlme) @@ -36,6 +39,7 @@ if(getRversion() >= '2.15.1') utils::globalVariables("ds") #' # tolerance in order to speed up the fit for this example evaluation #' f_nlme <- nlme(f, control = list(pnlsMaxIter = 120, tolerance = 1e-3)) #' plot(f_nlme) +#' } #' @export plot.nlme.mmkin <- function(x, i = 1:ncol(x$mmkin_orig), obs_vars = names(x$mkinmod$map), diff --git a/R/sigma_twocomp.R b/R/sigma_twocomp.R index e8a92ced..e7f4368b 100644 --- a/R/sigma_twocomp.R +++ b/R/sigma_twocomp.R @@ -1,16 +1,16 @@ #' Two-component error model -#' +#' #' Function describing the standard deviation of the measurement error in #' dependence of the measured value \eqn{y}: -#' +#' #' \deqn{\sigma = \sqrt{ \sigma_{low}^2 + y^2 * {rsd}_{high}^2}} sigma = #' sqrt(sigma_low^2 + y^2 * rsd_high^2) -#' +#' #' This is the error model used for example by Werner et al. (1978). The model #' proposed by Rocke and Lorenzato (1995) can be written in this form as well, #' but assumes approximate lognormal distribution of errors for high values of #' y. -#' +#' #' @param y The magnitude of the observed value #' @param sigma_low The asymptotic minimum of the standard deviation for low #' observed values @@ -20,7 +20,7 @@ #' @references Werner, Mario, Brooks, Samuel H., and Knott, Lancaster B. (1978) #' Additive, Multiplicative, and Mixed Analytical Errors. Clinical Chemistry #' 24(11), 1895-1898. -#' +#' #' Rocke, David M. and Lorenzato, Stefan (1995) A two-component model for #' measurement error in analytical chemistry. Technometrics 37(2), 176-184. #' @examples @@ -36,15 +36,8 @@ #' data = d_syn, na.action = na.omit, #' start = list(parent_0 = 100, lrc = -3)) #' if (length(findFunction("varConstProp")) > 0) { -#' f_gnls_tc <- gnls(value ~ SSasymp(time, 0, parent_0, lrc), -#' data = d_syn, na.action = na.omit, -#' start = list(parent_0 = 100, lrc = -3), -#' weights = varConstProp()) -#' f_gnls_tc_sf <- gnls(value ~ SSasymp(time, 0, parent_0, lrc), -#' data = d_syn, na.action = na.omit, -#' start = list(parent_0 = 100, lrc = -3), -#' control = list(sigma = 1), -#' weights = varConstProp()) +#' f_gnls_tc <- update(f_gnls, weights = varConstProp()) +#' f_gnls_tc_sf <- update(f_gnls_tc, control = list(sigma = 1)) #' } #' f_mkin <- mkinfit("SFO", d_syn, error_model = "const", quiet = TRUE) #' f_mkin_tc <- mkinfit("SFO", d_syn, error_model = "tc", quiet = TRUE) diff --git a/R/summary.nlme.mmkin.R b/R/summary.nlme.mmkin.R new file mode 100644 index 00000000..9fdd3f73 --- /dev/null +++ b/R/summary.nlme.mmkin.R @@ -0,0 +1,233 @@ +#' Summary method for class "nlme.mmkin" +#' +#' Lists model equations, initial parameter values, optimised parameters +#' for fixed effects (population), random effects (deviations from the +#' population mean) and residual error model, as well as the resulting +#' endpoints such as formation fractions and DT50 values. Optionally +#' (default is FALSE), the data are listed in full. +#' +#' @param object an object of class [nlme.mmkin] +#' @param x an object of class [summary.nlme.mmkin] +#' @param data logical, indicating whether the full data should be included in +#' the summary. +#' @param verbose Should the summary be verbose? +#' @param distimes logical, indicating whether DT50 and DT90 values should be +#' included. +#' @param alpha error level for confidence interval estimation from the t +#' distribution +#' @param digits Number of digits to use for printing +#' @param \dots optional arguments passed to methods like \code{print}. +#' @return The summary function returns a list based on the [nlme] object +#' obtained in the fit, with at least the following additional components +#' \item{nlmeversion, mkinversion, Rversion}{The nlme, mkin and R versions used} +#' \item{date.fit, date.summary}{The dates where the fit and the summary were +#' produced} +#' \item{diffs}{The differential equations used in the degradation model} +#' \item{use_of_ff}{Was maximum or minimum use made of formation fractions} +#' \item{data}{The data} +#' \item{confint_trans}{Transformed parameters as used in the optimisation, with confidence intervals} +#' \item{confint_back}{Backtransformed parameters, with confidence intervals if available} +#' \item{ff}{The estimated formation fractions derived from the fitted +#' model.} +#' \item{distimes}{The DT50 and DT90 values for each observed variable.} +#' \item{SFORB}{If applicable, eigenvalues of SFORB components of the model.} +#' The print method is called for its side effect, i.e. printing the summary. +#' @importFrom stats predict +#' @author Johannes Ranke for the mkin specific parts +#' José Pinheiro and Douglas Bates for the components inherited from nlme +#' @examples +#' +#' # Generate five datasets following SFO kinetics +#' sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120) +#' dt50_sfo_in_pop <- 50 +#' k_in_pop <- log(2) / dt50_sfo_in_pop +#' set.seed(1234) +#' k_in <- rlnorm(5, log(k_in_pop), 0.5) +#' 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) +#' names(ds_sfo_mean) <- paste("ds", 1:5) +#' +#' 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]] +#' }) +#' +#' # Evaluate using mmkin and nlme +#' library(nlme) +#' f_mmkin <- mmkin("SFO", ds_sfo_syn, quiet = TRUE, error_model = "tc", cores = 1) +#' f_nlme <- nlme(f_mmkin) +#' summary(f_nlme, data = TRUE) +#' +#' @export +summary.nlme.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes = TRUE, alpha = 0.05, ...) { + + mod_vars <- names(object$mkinmod$diffs) + + confint_trans <- intervals(object, which = "fixed", level = 1 - alpha)$fixed + attr(confint_trans, "label") <- NULL + pnames <- rownames(confint_trans) + confint_trans[, "est."] + bp <- backtransform_odeparms(confint_trans[, "est."], object$mkinmod, + object$transform_rates, object$transform_fractions) + bpnames <- names(bp) + + # variance-covariance estimates for fixed effects (from summary.lme) + fixed <- fixef(object) + stdFixed <- sqrt(diag(as.matrix(object$varFix))) + object$corFixed <- array( + t(object$varFix/stdFixed)/stdFixed, + dim(object$varFix), + list(names(fixed), names(fixed))) + + # 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) + for (box in mod_vars) { # Figure out sets of fractions to skip + f_names <- grep(paste("^f", box, sep = "_"), pnames, value = TRUE) + n_paths <- length(f_names) + if (n_paths > 1) f_names_skip <- c(f_names_skip, f_names) + } + + confint_back <- matrix(NA, nrow = length(bp), ncol = 3, + dimnames = list(bpnames, colnames(confint_trans))) + confint_back[, "est."] <- bp + + for (pname in pnames) { + if (!pname %in% f_names_skip) { + par.lower <- confint_trans[pname, "lower"] + par.upper <- confint_trans[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) + confint_back[names(bpl), "lower"] <- bpl + confint_back[names(bpu), "upper"] <- bpu + } + } + + object$confint_trans <- confint_trans + object$confint_back <- confint_back + + object$date.summary = date() + object$use_of_ff = object$mkinmod$use_of_ff + object$error_model_algorithm = object$mmkin_orig[[1]]$error_model_algorithm + err_mod = object$mmkin_orig[[1]]$err_mod + + object$diffs <- object$mkinmod$diffs + object$print_data <- data + if (data) { + object$data[["observed"]] <- object$data[["value"]] + object$data[["value"]] <- NULL + object$data[["predicted"]] <- predict(object) + object$data[["residual"]] <- residuals(object, type = "response") + object$data[["std"]] <- object$sigma <- 1/attr(object$modelStruct$varStruct, "weights") + object$data[["standardized"]] <- residuals(object, type = "pearson") + } + object$verbose <- verbose + + object$fixed <- object$mmkin_orig[[1]]$fixed + object$AIC = AIC(object) + object$BIC = BIC(object) + object$logLik = logLik(object) + + ep <- endpoints(object) + if (length(ep$ff) != 0) + object$ff <- ep$ff + if (distimes) object$distimes <- ep$distimes + if (length(ep$SFORB) != 0) object$SFORB <- ep$SFORB + class(object) <- c("summary.nlme.mmkin", "nlme.mmkin", "nlme", "lme") + return(object) +} + +#' @rdname summary.nlme.mmkin +#' @export +print.summary.nlme.mmkin <- function(x, digits = max(3, getOption("digits") - 3), verbose = x$verbose, ...) { + cat("nlme version used for fitting: ", x$nlmeversion, "\n") + cat("mkin version used for pre-fitting: ", x$mkinversion, "\n") + cat("R version used for fitting: ", x$Rversion, "\n") + + cat("Date of fit: ", x$date.fit, "\n") + cat("Date of summary:", x$date.summary, "\n") + + cat("\nEquations:\n") + nice_diffs <- gsub("^(d.*) =", "\\1/dt =", x[["diffs"]]) + writeLines(strwrap(nice_diffs, exdent = 11)) + + cat("\nData:\n") + cat(nrow(x$data), "observations of", + length(unique(x$data$name)), "variable(s) grouped in", + length(unique(x$data$ds)), "datasets\n") + + cat("\nModel predictions using solution type", x$solution_type, "\n") + + cat("\nFitted in", x$time[["elapsed"]], "s using", x$numIter, "iterations\n") + + cat("\nVariance model: ") + cat(switch(x$err_mod, + const = "Constant variance", + obs = "Variance unique to each observed variable", + tc = "Two-component variance function"), "\n") + + cat("\nMean of starting values for individual parameters:\n") + print(x$mean_dp_start) + + cat("\nFixed degradation parameter values:\n") + if(length(x$fixed$value) == 0) cat("None\n") + else print(x$fixed) + + cat("\nResults:\n\n") + print(data.frame(AIC = x$AIC, BIC = x$BIC, logLik = x$logLik, + row.names = " ")) + + cat("\nOptimised, transformed parameters with symmetric confidence intervals:\n") + print(x$confint_trans) + + if (nrow(x$confint_trans) > 1) { + corr <- x$corFixed + class(corr) <- "correlation" + print(corr, title = "\nCorrelation:", ...) + } + + cat("\nBacktransformed parameters with asymmetric confidence intervals:\n") + print(x$confint_back) + + print(summary(x$modelStruct), sigma = x$sigma, + reEstimates = x$coef$random, verbose = verbose, ...) + + printSFORB <- !is.null(x$SFORB) + if(printSFORB){ + cat("\nEstimated Eigenvalues of SFORB model(s):\n") + print(x$SFORB, digits=digits,...) + } + + printff <- !is.null(x$ff) + if(printff){ + cat("\nResulting formation fractions:\n") + print(data.frame(ff = x$ff), digits=digits,...) + } + + printdistimes <- !is.null(x$distimes) + if(printdistimes){ + cat("\nEstimated disappearance times:\n") + print(x$distimes, digits=digits,...) + } + + if (x$print_data){ + cat("\nData:\n") + print(format(x$data, digits = digits, ...), row.names = FALSE) + } + + invisible(x) +} diff --git a/_pkgdown.yml b/_pkgdown.yml index 4631dc97..ac100ac1 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -43,6 +43,7 @@ reference: contents: - nlme.mmkin - plot.nlme.mmkin + - summary.nlme.mmkin - nlme_function - get_deg_func - saemix_model diff --git a/docs/dev/news/index.html b/docs/dev/news/index.html index c29e413d..4e31f4b7 100644 --- a/docs/dev/news/index.html +++ b/docs/dev/news/index.html @@ -146,8 +146,7 @@ mkin 0.9.50.4 (unreleased) Unreleased diff --git a/docs/dev/pkgdown.yml b/docs/dev/pkgdown.yml index 8f493a24..657aa128 100644 --- a/docs/dev/pkgdown.yml +++ b/docs/dev/pkgdown.yml @@ -10,7 +10,7 @@ articles: web_only/NAFTA_examples: NAFTA_examples.html web_only/benchmarks: benchmarks.html web_only/compiled_models: compiled_models.html -last_built: 2020-10-26T13:18Z +last_built: 2020-10-27T14:34Z urls: reference: https://pkgdown.jrwb.de/mkin/reference article: https://pkgdown.jrwb.de/mkin/articles diff --git a/docs/dev/reference/Rplot001.png b/docs/dev/reference/Rplot001.png index cfc5bc2b..f001da49 100644 Binary files a/docs/dev/reference/Rplot001.png and b/docs/dev/reference/Rplot001.png differ diff --git a/docs/dev/reference/index.html b/docs/dev/reference/index.html index 24056025..fa3ec868 100644 --- a/docs/dev/reference/index.html +++ b/docs/dev/reference/index.html @@ -308,7 +308,7 @@ of an mmkin object

Mixed models

-

Create and work with nonlinear mixed models

+

Create and work with nonlinear mixed effects models

@@ -330,6 +330,12 @@ of an mmkin object

Plot a fitted nonlinear mixed model obtained via an mmkin row object

+ +

summary(<nlme.mmkin>) print(<summary.nlme.mmkin>)

+ +

Summary method for class "nlme.mmkin"

+ +

nlme_function() mean_degparms() nlme_data()

diff --git a/docs/dev/reference/nlme.mmkin.html b/docs/dev/reference/nlme.mmkin.html index 90aec9be..3fa4d97b 100644 --- a/docs/dev/reference/nlme.mmkin.html +++ b/docs/dev/reference/nlme.mmkin.html @@ -152,7 +152,7 @@ have been obtained by fitting the same model to a list of datasets.

# S3 method for mmkin
-nlme(
+nlme(
   model,
   data = sys.frame(sys.parent()),
   fixed,
@@ -255,6 +255,11 @@ parameters taken from the mmkin object are used

 
     
Upon success, a fitted nlme.mmkin object, which is an nlme object
 with additional elements

+    
Note

+
+    
As the object inherits from nlme::nlme, there is a wealth of
+methods that will automatically work on 'nlme.mmkin' objects, such as
+nlme::intervals(), nlme::anova.lme() and nlme::coef.lme().

     
See also

 
     
nlme_function

@@ -270,11 +275,20 @@ with additional elements

 
#>             df      AIC
 #> f_nlme_sfo   5 625.0539
 #> f_nlme_dfop  9 495.1270
print(f_nlme_dfop)
-
#> Nonlinear mixed-effects model fit by maximum likelihood
-#>   Model: value ~ (mkin::get_deg_func())(name, time, parent_0, log_k1,      log_k2, g_ilr) 
-#>   Data: "Not shown" 
-#>   Log-likelihood: -238.5635
-#>   Fixed: list(parent_0 ~ 1, log_k1 ~ 1, log_k2 ~ 1, g_ilr ~ 1) 
+
#> Kinetic nonlinear mixed-effects model fit by maximum likelihood
+#> 
+#> Structural model:
+#> d_parent/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
+#>            time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time)))
+#>            * parent
+#> 
+#> Data:
+#> 90 observations of 1 variable(s) grouped in 5 datasets
+#> 
+#> Log-likelihood: -238.5635
+#> 
+#> Fixed effects:
+#>  list(parent_0 ~ 1, log_k1 ~ 1, log_k2 ~ 1, g_ilr ~ 1) 
 #>    parent_0      log_k1      log_k2       g_ilr 
 #> 94.17015133 -1.80015306 -4.14738870  0.02290935 
 #> 
@@ -284,9 +298,7 @@ with additional elements

 #>  Structure: Diagonal
 #>         parent_0    log_k1  log_k2     g_ilr Residual
 #> StdDev: 2.488249 0.8447273 1.32965 0.3289311 2.321364
-#> 
-#> Number of Observations: 90
-#> Number of Groups: 5 
plot(f_nlme_dfop)
+#> 
plot(f_nlme_dfop)
 
endpoints(f_nlme_dfop)
 
#> $distimes
 #>            DT50     DT90 DT50back  DT50_k1  DT50_k2
@@ -295,11 +307,18 @@ with additional elements

 # \dontrun{
   f_nlme_2 <- nlme(f["SFO", ], start = c(parent_0 = 100, log_k_parent = 0.1))
   update(f_nlme_2, random = parent_0 ~ 1)
-
#> Nonlinear mixed-effects model fit by maximum likelihood
-#>   Model: value ~ (mkin::get_deg_func())(name, time, parent_0, log_k_parent) 
-#>   Data: "Not shown" 
-#>   Log-likelihood: -404.3729
-#>   Fixed: list(parent_0 ~ 1, log_k_parent ~ 1) 
+
#> Kinetic nonlinear mixed-effects model fit by maximum likelihood
+#> 
+#> Structural model:
+#> d_parent/dt = - k_parent * parent
+#> 
+#> Data:
+#>  observations of 0 variable(s) grouped in 0 datasets
+#> 
+#> Log-likelihood: -404.3729
+#> 
+#> Fixed effects:
+#>  list(parent_0 ~ 1, log_k_parent ~ 1) 
 #>     parent_0 log_k_parent 
 #>    75.933480    -3.555983 
 #> 
@@ -307,9 +326,7 @@ with additional elements

 #>  Formula: parent_0 ~ 1 | ds
 #>            parent_0 Residual
 #> StdDev: 0.002416792 21.63027
-#> 
-#> Number of Observations: 90
-#> Number of Groups: 5 
  ds_2 <- lapply(experimental_data_for_UBA_2019[6:10],
+#> 
  ds_2 <- lapply(experimental_data_for_UBA_2019[6:10],
    function(x) x$data[c("name", "time", "value")])
   m_sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"),
     A1 = mkinsub("SFO"), use_of_ff = "min", quiet = TRUE)
@@ -395,11 +412,20 @@ with additional elements

     AIC(f_nlme_sfo, f_nlme_sfo_tc, f_nlme_dfop, f_nlme_dfop_tc)
     print(f_nlme_dfop_tc)
   }
-
#> Nonlinear mixed-effects model fit by maximum likelihood
-#>   Model: value ~ (mkin::get_deg_func())(name, time, parent_0, log_k1,      log_k2, g_ilr) 
-#>   Data: "Not shown" 
-#>   Log-likelihood: -238.4298
-#>   Fixed: list(parent_0 ~ 1, log_k1 ~ 1, log_k2 ~ 1, g_ilr ~ 1) 
+
#> Kinetic nonlinear mixed-effects model fit by maximum likelihood
+#> 
+#> Structural model:
+#> d_parent/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
+#>            time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time)))
+#>            * parent
+#> 
+#> Data:
+#> 90 observations of 1 variable(s) grouped in 5 datasets
+#> 
+#> Log-likelihood: -238.4298
+#> 
+#> Fixed effects:
+#>  list(parent_0 ~ 1, log_k1 ~ 1, log_k2 ~ 1, g_ilr ~ 1) 
 #>    parent_0      log_k1      log_k2       g_ilr 
 #> 94.04774463 -1.82339924 -4.16715509  0.04020161 
 #> 
@@ -415,19 +441,25 @@ with additional elements

 #>  Formula: ~fitted(.) 
 #>  Parameter estimates:
 #>      const       prop 
-#> 2.23222625 0.01262414 
-#> Number of Observations: 90
-#> Number of Groups: 5 

+#> 2.23222625 0.01262414 

   f_2_obs <- mmkin(list("SFO-SFO" = m_sfo_sfo,
    "DFOP-SFO" = m_dfop_sfo),
     ds_2, quiet = TRUE, error_model = "obs")
   f_nlme_sfo_sfo_obs <- nlme(f_2_obs["SFO-SFO", ])
   print(f_nlme_sfo_sfo_obs)
-
#> Nonlinear mixed-effects model fit by maximum likelihood
-#>   Model: value ~ (mkin::get_deg_func())(name, time, parent_0, log_k_parent_sink,      log_k_parent_A1, log_k_A1_sink) 
-#>   Data: "Not shown" 
-#>   Log-likelihood: -472.976
-#>   Fixed: list(parent_0 ~ 1, log_k_parent_sink ~ 1, log_k_parent_A1 ~ 1,      log_k_A1_sink ~ 1) 
+
#> Kinetic nonlinear mixed-effects model fit by maximum likelihood
+#> 
+#> Structural model:
+#> d_parent/dt = - k_parent_sink * parent - k_parent_A1 * parent
+#> d_A1/dt = + k_parent_A1 * parent - k_A1_sink * A1
+#> 
+#> Data:
+#> 170 observations of 2 variable(s) grouped in 5 datasets
+#> 
+#> Log-likelihood: -472.976
+#> 
+#> Fixed effects:
+#>  list(parent_0 ~ 1, log_k_parent_sink ~ 1, log_k_parent_A1 ~ 1,      log_k_A1_sink ~ 1) 
 #>          parent_0 log_k_parent_sink   log_k_parent_A1     log_k_A1_sink 
 #>         87.975536         -3.669816         -4.164127         -4.645073 
 #> 
@@ -443,9 +475,7 @@ with additional elements

 #>  Formula: ~1 | name 
 #>  Parameter estimates:
 #>    parent        A1 
-#> 1.0000000 0.2050003 
-#> Number of Observations: 170
-#> Number of Groups: 5 
  # The same with DFOP-SFO does not converge, apparently the variances of
+#> 1.0000000 0.2050003 
  # The same with DFOP-SFO does not converge, apparently the variances of
   # parent and A1 are too similar in this case, so that the model is
   # overparameterised
   #f_nlme_dfop_sfo_obs <- nlme(f_2_obs["DFOP-SFO", ], control = list(maxIter = 100))
diff --git a/docs/dev/reference/plot.nlme.mmkin.html b/docs/dev/reference/plot.nlme.mmkin.html
index afd9d8d0..267bef05 100644
--- a/docs/dev/reference/plot.nlme.mmkin.html
+++ b/docs/dev/reference/plot.nlme.mmkin.html
@@ -214,6 +214,14 @@ predicted values?

       
Maximum absolute value of the residuals. This is used for the
 scaling of the y axis and defaults to "auto".

     
+    
+      ncol.legend
+      
Number of columns to use in the legend

+    
+    
+      nrow.legend
+      
Number of rows to use in the legend

+    
     
       rel.height.legend
       
The relative height of the legend shown on top

@@ -259,7 +267,8 @@ corresponding model prediction lines for the different datasets.

 names(ds) <- paste0("ds ", 6:10)
 dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"),
   A1 = mkinsub("SFO"), quiet = TRUE)
-f <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE, cores = 1)
+# \dontrun{
+f <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE)
 plot(f[, 3:4], standardized = TRUE)
 

 library(nlme)
@@ -267,7 +276,8 @@ corresponding model prediction lines for the different datasets.

 # tolerance in order to speed up the fit for this example evaluation
 f_nlme <- nlme(f, control = list(pnlsMaxIter = 120, tolerance = 1e-3))
 plot(f_nlme)
-
+
# } +