From 675a733fa2acc08daabb9b8b571c7d658f281f73 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Tue, 26 May 2020 18:38:51 +0200 Subject: Use all cores per default, confint tolerance Also, use more intelligent starting values for the variance of the random effects for saemix. While this does not appear to speed up the convergence, it shows where this variance is greatly reduced by using mixed-effects models as opposed to the separate independent fits. --- NAMESPACE | 1 + NEWS.md | 4 + R/confint.mkinfit.R | 86 +++---- R/mmkin.R | 5 +- R/nlme.R | 2 +- R/parms.mkinfit.R | 4 +- R/saemix.R | 31 ++- build.log | 2 +- check.log | 2 +- docs/news/index.html | 2 + docs/pkgdown.yml | 2 +- docs/reference/confint.mkinfit.html | 20 +- docs/reference/mmkin.html | 16 +- docs/reference/nlme.html | 4 +- docs/reference/parms.html | 6 +- docs/reference/saemix-1.png | Bin 37443 -> 31551 bytes docs/reference/saemix-2.png | Bin 38557 -> 58815 bytes docs/reference/saemix.html | 135 ++++------- man/confint.mkinfit.Rd | 7 +- man/mmkin.Rd | 5 +- man/nlme.Rd | 2 +- man/parms.Rd | 4 +- man/saemix.Rd | 23 +- test.log | 20 +- tests/testthat/FOCUS_2006_D.csf | 2 +- vignettes/FOCUS_D.html | 10 +- vignettes/FOCUS_L.html | 453 +++++++++++++++++++----------------- vignettes/mkin.html | 4 +- vignettes/twa.html | 19 +- 29 files changed, 466 insertions(+), 405 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index d2298dec..f8be4fae 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -116,5 +116,6 @@ importFrom(stats,qt) importFrom(stats,residuals) importFrom(stats,rnorm) importFrom(stats,update) +importFrom(stats,var) importFrom(utils,getFromNamespace) importFrom(utils,write.table) diff --git a/NEWS.md b/NEWS.md index a3bd2bc2..8f88b64d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,6 +4,10 @@ - 'saemix_model', 'saemix_data': Helper functions to fit nonlinear mixed-effects models for mmkin row objects using the saemix package +- 'mmkin' and 'confint(method = 'profile'): Use all cores detected by parallel::detectCores() per default + +- 'confint(method = 'profile'): Choose accuracy based on 'rel_tol' argument, relative to the bounds obtained by the quadratic approximation + # mkin 0.9.50.2 (2020-05-12) - Increase tolerance for a platform specific test results on the Solaris test machine on CRAN diff --git a/R/confint.mkinfit.R b/R/confint.mkinfit.R index 78dda95d..53eb45ee 100644 --- a/R/confint.mkinfit.R +++ b/R/confint.mkinfit.R @@ -27,7 +27,10 @@ #' @param backtransform If we approximate the likelihood in terms of the #' transformed parameters, should we backtransform the parameters with #' their confidence intervals? -#' @param cores The number of cores to be used for multicore processing. +#' @param rel_tol If the method is 'profile', what should be the accuracy +#' of the lower and upper bounds, relative to the estimate obtained from +#' the quadratic method? +#' @param cores The number of cores to be used for multicore processing. #' On Windows machines, cores > 1 is currently not supported. #' @param quiet Should we suppress the message "Profiling the likelihood" #' @return A matrix with columns giving lower and upper confidence limits for @@ -121,7 +124,7 @@ confint.mkinfit <- function(object, parm, level = 0.95, alpha = 1 - level, cutoff, method = c("quadratic", "profile"), transformed = TRUE, backtransform = TRUE, - cores = round(detectCores()/2), quiet = FALSE, ...) + cores = parallel::detectCores(), rel_tol = 0.01, quiet = FALSE, ...) { tparms <- parms(object, transformed = TRUE) bparms <- parms(object, transformed = FALSE) @@ -140,50 +143,50 @@ confint.mkinfit <- function(object, parm, a <- c(alpha / 2, 1 - (alpha / 2)) - if (method == "quadratic") { + quantiles <- qt(a, object$df.residual) - quantiles <- qt(a, object$df.residual) - - covar_pnames <- if (missing(parm)) { - if (transformed) tpnames else bpnames - } else { - parm - } + covar_pnames <- if (missing(parm)) { + if (transformed) tpnames else bpnames + } else { + parm + } - return_parms <- if (backtransform) bparms[return_pnames] - else tparms[return_pnames] + return_parms <- if (backtransform) bparms[return_pnames] + else tparms[return_pnames] - covar_parms <- if (transformed) tparms[covar_pnames] - else bparms[covar_pnames] + covar_parms <- if (transformed) tparms[covar_pnames] + else bparms[covar_pnames] - if (transformed) { - covar <- try(solve(object$hessian), silent = TRUE) - } else { - covar <- try(solve(object$hessian_notrans), silent = TRUE) - } + if (transformed) { + covar <- try(solve(object$hessian), silent = TRUE) + } else { + covar <- try(solve(object$hessian_notrans), silent = TRUE) + } - # If inverting the covariance matrix failed or produced NA values - if (!is.numeric(covar) | is.na(covar[1])) { - ses <- lci <- uci <- rep(NA, p) - } else { - ses <- sqrt(diag(covar))[covar_pnames] - lci <- covar_parms + quantiles[1] * ses - uci <- covar_parms + quantiles[2] * ses - if (transformed & backtransform) { - lci_back <- backtransform_odeparms(lci, - object$mkinmod, object$transform_rates, object$transform_fractions) - uci_back <- backtransform_odeparms(uci, - object$mkinmod, object$transform_rates, object$transform_fractions) - - return_errparm_names <- intersect(names(object$errparms), return_pnames) - lci <- c(lci_back, lci[return_errparm_names]) - uci <- c(uci_back, uci[return_errparm_names]) - } + # If inverting the covariance matrix failed or produced NA values + if (!is.numeric(covar) | is.na(covar[1])) { + ses <- lci <- uci <- rep(NA, p) + } else { + ses <- sqrt(diag(covar))[covar_pnames] + lci <- covar_parms + quantiles[1] * ses + uci <- covar_parms + quantiles[2] * ses + if (transformed & backtransform) { + lci_back <- backtransform_odeparms(lci, + object$mkinmod, object$transform_rates, object$transform_fractions) + uci_back <- backtransform_odeparms(uci, + object$mkinmod, object$transform_rates, object$transform_fractions) + + return_errparm_names <- intersect(names(object$errparms), return_pnames) + lci <- c(lci_back, lci[return_errparm_names]) + uci <- c(uci_back, uci[return_errparm_names]) } - ci <- cbind(lower = lci, upper = uci) } + ci <- cbind(lower = lci, upper = uci) if (method == "profile") { + + ci_quadratic <- ci + if (!quiet) message("Profiling the likelihood") lci <- uci <- rep(NA, p) @@ -215,9 +218,14 @@ confint.mkinfit <- function(object, parm, (cutoff - (object$logLik - profile_ll(x)))^2 } - lci_pname <- optimize(cost, lower = 0, upper = all_parms[pname])$minimum + lower_quadratic <- ci_quadratic["lower"][pname] + upper_quadratic <- ci_quadratic["upper"][pname] + ltol <- if (!is.na(lower_quadratic)) rel_tol * lower_quadratic else .Machine$double.eps^0.25 + utol <- if (!is.na(upper_quadratic)) rel_tol * upper_quadratic else .Machine$double.eps^0.25 + lci_pname <- optimize(cost, lower = 0, upper = all_parms[pname], tol = ltol)$minimum uci_pname <- optimize(cost, lower = all_parms[pname], - upper = ifelse(grepl("^f_|^g$", pname), 1, 15 * all_parms[pname]))$minimum + upper = ifelse(grepl("^f_|^g$", pname), 1, 15 * all_parms[pname]), + tol = utol)$minimum return(c(lci_pname, uci_pname)) } ci <- t(parallel::mcmapply(get_ci, profile_pnames, mc.cores = cores)) diff --git a/R/mmkin.R b/R/mmkin.R index dbb61b78..37c4e87d 100644 --- a/R/mmkin.R +++ b/R/mmkin.R @@ -12,7 +12,8 @@ #' @param cores The number of cores to be used for multicore processing. This #' is only used when the \code{cluster} argument is \code{NULL}. On Windows #' machines, cores > 1 is not supported, you need to use the \code{cluster} -#' argument to use multiple logical processors. +#' argument to use multiple logical processors. Per default, all cores +#' detected by [parallel::detectCores()] are used. #' @param cluster A cluster as returned by \code{\link{makeCluster}} to be used #' for parallel execution. #' @param \dots Further arguments that will be passed to \code{\link{mkinfit}}. @@ -62,7 +63,7 @@ #' #' @export mmkin mmkin <- function(models = c("SFO", "FOMC", "DFOP"), datasets, - cores = round(detectCores()/2), cluster = NULL, ...) + cores = detectCores(), cluster = NULL, ...) { parent_models_available = c("SFO", "FOMC", "DFOP", "HS", "SFORB", "IORE", "logistic") n.m <- length(models) diff --git a/R/nlme.R b/R/nlme.R index 3ee7b9fd..20987064 100644 --- a/R/nlme.R +++ b/R/nlme.R @@ -125,7 +125,7 @@ nlme_function <- function(object) { #' @return If random is FALSE (default), a named vector containing mean values #' of the fitted degradation model parameters. If random is TRUE, a list with #' fixed and random effects, in the format required by the start argument of -#' nlme for the case of a single grouping variable ds? +#' nlme for the case of a single grouping variable ds. #' @param random Should a list with fixed and random effects be returned? #' @export mean_degparms <- function(object, random = FALSE) { diff --git a/R/parms.mkinfit.R b/R/parms.mkinfit.R index aae6fa52..a1f2d209 100644 --- a/R/parms.mkinfit.R +++ b/R/parms.mkinfit.R @@ -21,11 +21,13 @@ #' ds <- lapply(experimental_data_for_UBA_2019[6:10], #' function(x) subset(x$data[c("name", "time", "value")])) #' names(ds) <- paste("Dataset", 6:10) -#' fits <- mmkin(c("SFO", "FOMC", "DFOP"), ds, quiet = TRUE) +#' \dontrun{ +#' fits <- mmkin(c("SFO", "FOMC", "DFOP"), ds, quiet = TRUE, cores = 1) #' parms(fits["SFO", ]) #' parms(fits[, 2]) #' parms(fits) #' parms(fits, transformed = TRUE) +#' } #' @export parms <- function(object, ...) { diff --git a/R/saemix.R b/R/saemix.R index 69e5fc50..24c0f0d0 100644 --- a/R/saemix.R +++ b/R/saemix.R @@ -5,25 +5,37 @@ #' list of mkinfit objects that have been obtained by fitting the same model to #' a list of datasets. #' +#' Starting values for the fixed effects (population mean parameters, argument psi0 of +#' [saemix::saemixModel()] are the mean values of the parameters found using +#' mmkin. Starting variances of the random effects (argument omega.init) are the +#' variances of the deviations of the parameters from these mean values. +#' #' @param object An mmkin row object containing several fits of the same model to different datasets +#' @param cores The number of cores to be used for multicore processing. +#' On Windows machines, cores > 1 is currently not supported. #' @rdname saemix #' @importFrom saemix saemixData saemixModel +#' @importFrom stats var #' @examples #' 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"), #' A1 = mkinsub("SFO")) -#' f_mmkin <- mmkin(list("SFO-SFO" = sfo_sfo), ds, quiet = TRUE, cores = 5) +#' \dontrun{ +#' f_mmkin <- mmkin(list("SFO-SFO" = sfo_sfo), ds, quiet = TRUE) +#' library(saemix) #' m_saemix <- saemix_model(f_mmkin) #' d_saemix <- saemix_data(f_mmkin) -#' saemix_options <- list(seed = 123456, save = FALSE, save.graphs = FALSE) -#' \dontrun{ -#' saemix(m_saemix, d_saemix, saemix_options) +#' 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") #' } #' @return An [saemix::SaemixModel] object. #' @export -saemix_model <- function(object) { +saemix_model <- function(object, cores = parallel::detectCores()) { if (nrow(object) > 1) stop("Only row objects allowed") mkin_model <- object[[1]]$mkinmod @@ -81,14 +93,19 @@ saemix_model <- function(object) { out_values <- out_wide[out_index] } return(out_values) - }, mc.cores = 15) + }, mc.cores = cores) res <- unlist(res_list) return(res) } + raneff_0 <- mean_degparms(object, random = TRUE)$random$ds + var_raneff_0 <- apply(raneff_0, 2, var) + res <- saemixModel(model_function, psi0, "Mixed model generated from mmkin object", - transform.par = rep(0, length(degparms_optim))) + transform.par = rep(0, length(degparms_optim)), + omega.init = diag(var_raneff_0) + ) return(res) } diff --git a/build.log b/build.log index bd53dcee..b94e6450 100644 --- a/build.log +++ b/build.log @@ -5,5 +5,5 @@ * creating vignettes ... OK * checking for LF line-endings in source and make files and shell scripts * checking for empty or unneeded directories -* building ‘mkin_0.9.50.2.tar.gz’ +* building ‘mkin_0.9.50.3.tar.gz’ diff --git a/check.log b/check.log index 17413d04..cae31a24 100644 --- a/check.log +++ b/check.log @@ -5,7 +5,7 @@ * using options ‘--no-tests --as-cran’ * checking for file ‘mkin/DESCRIPTION’ ... OK * checking extension type ... Package -* this is package ‘mkin’ version ‘0.9.50.2’ +* this is package ‘mkin’ version ‘0.9.50.3’ * package encoding: UTF-8 * checking CRAN incoming feasibility ... Note_to_CRAN_maintainers Maintainer: ‘Johannes Ranke ’ diff --git a/docs/news/index.html b/docs/news/index.html index 149fc98e..c26652e9 100644 --- a/docs/news/index.html +++ b/docs/news/index.html @@ -148,6 +148,8 @@
diff --git a/docs/pkgdown.yml b/docs/pkgdown.yml index 0bb01ef4..a8b168ce 100644 --- a/docs/pkgdown.yml +++ b/docs/pkgdown.yml @@ -10,7 +10,7 @@ articles: NAFTA_examples: web_only/NAFTA_examples.html benchmarks: web_only/benchmarks.html compiled_models: web_only/compiled_models.html -last_built: 2020-05-25T10:48Z +last_built: 2020-05-26T16:38Z urls: reference: https://pkgdown.jrwb.de/mkin/reference article: https://pkgdown.jrwb.de/mkin/articles diff --git a/docs/reference/confint.mkinfit.html b/docs/reference/confint.mkinfit.html index 190494bc..0686c7bb 100644 --- a/docs/reference/confint.mkinfit.html +++ b/docs/reference/confint.mkinfit.html @@ -79,7 +79,7 @@ method of Venzon and Moolgavkar (1988)." /> mkin - 0.9.50.2 + 0.9.50.3
@@ -116,6 +116,9 @@ method of Venzon and Moolgavkar (1988)." />
  • Example evaluation of NAFTA SOP Attachment examples
    • +
  • + Some benchmark timings +
  • @@ -168,7 +171,8 @@ method of Venzon and Moolgavkar (1988).

    method = c("quadratic", "profile"), transformed = TRUE, backtransform = TRUE, - cores = round(detectCores()/2), + cores = parallel::detectCores(), + rel_tol = 0.01, quiet = FALSE, ... ) @@ -222,6 +226,12 @@ their confidence intervals?

    cores

    The number of cores to be used for multicore processing. On Windows machines, cores > 1 is currently not supported.

    + + + rel_tol +

    If the method is 'profile', what should be the accuracy +of the lower and upper bounds, relative to the estimate obtained from +the quadratic method?

    quiet @@ -270,13 +280,13 @@ Profile-Likelihood Based Confidence Intervals, Applied Statistics, 37, SFO_SFO.ff <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO"), use_of_ff = "max", quiet = TRUE) f_d_1 <- mkinfit(SFO_SFO, subset(FOCUS_2006_D, value != 0), quiet = TRUE) -system.time(ci_profile <- confint(f_d_1, method = "profile", cores = 1, quiet = TRUE))
    #> User System verstrichen -#> 3.410 0.000 3.412
    # Using more cores does not save much time here, as parent_0 takes up most of the time +system.time(ci_profile <- confint(f_d_1, method = "profile", cores = 1, quiet = TRUE))
    #> user system elapsed +#> 3.689 0.991 3.361
    # Using more cores does not save much time here, as parent_0 takes up most of the time # If we additionally exclude parent_0 (the confidence of which is often of # minor interest), we get a nice performance improvement from about 50 # seconds to about 12 seconds if we use at least four cores system.time(ci_profile_no_parent_0 <- confint(f_d_1, method = "profile", - c("k_parent_sink", "k_parent_m1", "k_m1_sink", "sigma"), cores = n_cores))
    #> Profiling the likelihood
    #> Warning: scheduled cores 1, 2, 3 encountered errors in user code, all values of the jobs will be affected
    #> Error in dimnames(x) <- dn: Länge von 'dimnames' [2] ungleich der Arrayausdehnung
    #> Timing stopped at: 0.008 0.044 0.201
    ci_profile
    #> 2.5% 97.5% + c("k_parent_sink", "k_parent_m1", "k_m1_sink", "sigma"), cores = n_cores))
    #> Profiling the likelihood
    #> Warning: scheduled cores 2, 1, 3 encountered errors in user code, all values of the jobs will be affected
    #> Error in dimnames(x) <- dn: length of 'dimnames' [2] not equal to array extent
    #> Timing stopped at: 0.007 0.042 0.193
    ci_profile
    #> 2.5% 97.5% #> parent_0 96.456003640 1.027703e+02 #> k_parent 0.090911032 1.071578e-01 #> k_m1 0.003892605 6.702778e-03 diff --git a/docs/reference/mmkin.html b/docs/reference/mmkin.html index 3be3b4b9..9628c017 100644 --- a/docs/reference/mmkin.html +++ b/docs/reference/mmkin.html @@ -75,7 +75,7 @@ datasets specified in its first two arguments." /> mkin - 0.9.50.2 + 0.9.50.3
    @@ -112,6 +112,9 @@ datasets specified in its first two arguments." />
  • Example evaluation of NAFTA SOP Attachment examples
    • +
  • + Some benchmark timings +
  • @@ -152,7 +155,7 @@ datasets specified in its first two arguments.

    mmkin(
   models = c("SFO", "FOMC", "DFOP"),
   datasets,
-  cores = round(detectCores()/2),
+  cores = detectCores(),
   cluster = NULL,
   ...
 )
    @@ -176,7 +179,8 @@ data for mkinfit.

    The number of cores to be used for multicore processing. This is only used when the cluster argument is NULL. On Windows machines, cores > 1 is not supported, you need to use the cluster -argument to use multiple logical processors.

    +argument to use multiple logical processors. Per default, all cores +detected by parallel::detectCores() are used.

    cluster @@ -215,9 +219,9 @@ plotting.

    time_default <- system.time(fits.0 <- mmkin(models, datasets, quiet = TRUE)) time_1 <- system.time(fits.4 <- mmkin(models, datasets, cores = 1, quiet = TRUE))
    #> Warning: Optimisation did not converge: #> false convergence (8)
    -time_default
    #> User System verstrichen -#> 4.520 0.374 1.284
    time_1
    #> User System verstrichen -#> 5.076 0.004 5.083
    +time_default
    #> user system elapsed +#> 4.457 0.561 1.328
    time_1
    #> user system elapsed +#> 5.031 0.004 5.038
    endpoints(fits.0[["SFO_lin", 2]])
    #> $ff #> parent_M1 parent_sink M1_M2 M1_sink #> 0.7340478 0.2659522 0.7505691 0.2494309 diff --git a/docs/reference/nlme.html b/docs/reference/nlme.html index b2d415dc..3462e52e 100644 --- a/docs/reference/nlme.html +++ b/docs/reference/nlme.html @@ -75,7 +75,7 @@ datasets. They are used internally by the nlme.mmkin() method." /> mkin - 0.9.50.2 + 0.9.50.3
    @@ -178,7 +178,7 @@ datasets. They are used internally by the nlme.m

    If random is FALSE (default), a named vector containing mean values of the fitted degradation model parameters. If random is TRUE, a list with fixed and random effects, in the format required by the start argument of -nlme for the case of a single grouping variable ds?

    +nlme for the case of a single grouping variable ds.

    A groupedData object

    See also

    diff --git a/docs/reference/parms.html b/docs/reference/parms.html index 432bbc88..2fe91c26 100644 --- a/docs/reference/parms.html +++ b/docs/reference/parms.html @@ -195,7 +195,8 @@ such matrices is returned.

    ds <- lapply(experimental_data_for_UBA_2019[6:10], function(x) subset(x$data[c("name", "time", "value")])) names(ds) <- paste("Dataset", 6:10) -fits <- mmkin(c("SFO", "FOMC", "DFOP"), ds, quiet = TRUE) +# \dontrun{ +fits <- mmkin(c("SFO", "FOMC", "DFOP"), ds, quiet = TRUE, cores = 1) parms(fits["SFO", ])
    #> Dataset 6 Dataset 7 Dataset 8 Dataset 9 Dataset 10 #> parent_0 88.52275400 82.666781678 86.8547308 91.7779306 82.14809450 #> k_parent_sink 0.05794659 0.009647805 0.2102974 0.1232258 0.00720421 @@ -259,7 +260,8 @@ such matrices is returned.

    #> log_k2 -3.5206791 -5.85402317 -2.5794240 -3.4233253 -5.676532 #> g_ilr -0.1463234 0.07627854 0.4719196 0.4477805 -0.460676 #> sigma 1.3569047 2.22130220 1.3416908 2.8715985 1.942068 -#>
    +#>
    # } +
    - 
    saemix_model(object)
+    
    saemix_model(object, cores = parallel::detectCores())
 
 saemix_data(object, ...)
    
 
@@ -164,6 +164,11 @@ a list of datasets.
    
       object
       
    An mmkin row object containing several fits of the same model to different datasets
    
     
+    
+      cores
+      
    The number of cores to be used for multicore processing.
+On Windows machines, cores > 1 is currently not supported.
    
+    
     
       ...
       
    Further parameters passed to saemix::saemixData
    
@@ -174,21 +179,22 @@ a list of datasets.
    
 
     
    An saemix::SaemixModel object.
    
 
    An saemix::SaemixData object.
    
+    
    Details
    
+
+    
    Starting values for the fixed effects (population mean parameters, argument psi0 of
+saemix::saemixModel() are the mean values of the parameters found using
+mmkin. Starting variances of the random effects (argument omega.init) are the
+variances of the deviations of the parameters from these mean values.
    
 
     
    Examples
    
     
    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"),
-  A1 = mkinsub("SFO"))
    #> Successfully compiled differential equation model from auto-generated C code.
    f_mmkin <- mmkin(list("SFO-SFO" = sfo_sfo), ds, quiet = TRUE, cores = 5)
-# \dontrun{
-if (require(saemix)) {
-  m_saemix <- saemix_model(f_mmkin)
-  d_saemix <- saemix_data(f_mmkin)
-  saemix_options <- list(seed = 123456, save = FALSE, save.graphs = FALSE)
-    saemix(m_saemix, d_saemix, saemix_options)
-}
    #> Loading required package: saemix
    #> Package saemix, version 3.1.9000
-#>   please direct bugs, questions and feedback to emmanuelle.comets@inserm.fr
    #> 
+  A1 = mkinsub("SFO"))
    #> Successfully compiled differential equation model from auto-generated C code.
    # \dontrun{
+f_mmkin <- mmkin(list("SFO-SFO" = sfo_sfo), ds, quiet = TRUE)
+library(saemix)
    #> Package saemix, version 3.1.9000
+#>   please direct bugs, questions and feedback to emmanuelle.comets@inserm.fr
    m_saemix <- saemix_model(f_mmkin)
    #> 
 #> 
 #> The following SaemixModel object was successfully created:
 #> 
@@ -224,12 +230,12 @@ a list of datasets.
    
 #>             out_values <- out_wide[out_index]
 #>         }
 #>         return(out_values)
-#>     }, mc.cores = 15)
+#>     }, mc.cores = cores)
 #>     res <- unlist(res_list)
 #>     return(res)
 #> }
-#> <bytecode: 0x555559875398>
-#> <environment: 0x55555973a248>
+#> <bytecode: 0x555559668108>
+#> <environment: 0x555559677c08>
 #>   Nb of parameters: 4 
 #>       parameter names:  parent_0 log_k_parent log_k_A1 f_parent_ilr_1 
 #>       distribution:
@@ -248,8 +254,7 @@ a list of datasets.
    
 #>     No covariate in the model.
 #>     Initial values
 #>              parent_0 log_k_parent  log_k_A1 f_parent_ilr_1
-#> Pop.CondInit 86.53449    -3.207005 -3.060308      -1.920449
-#> 
+#> Pop.CondInit 86.53449    -3.207005 -3.060308      -1.920449
    d_saemix <- saemix_data(f_mmkin)
    #> 
 #> 
 #> The following SaemixData object was successfully created:
 #> 
@@ -257,12 +262,14 @@ a list of datasets.
    
 #>     longitudinal data for use with the SAEM algorithm
 #> Dataset ds_saemix 
 #>     Structured data: value ~ time + name | ds 
-#>     X variable for graphs: time () 
-#> Running main SAEM algorithm
-#> [1] "Mon May 25 12:48:51 2020"
-#> .
    #> .
    #> .
    #> .
    #> 
+#>     X variable for graphs: time () 
    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)
    #> Running main SAEM algorithm
+#> [1] "Tue May 26 18:25:16 2020"
+#> ..
 #>     Minimisation finished
-#> [1] "Mon May 25 12:56:39 2020"
    #> Nonlinear mixed-effects model fit by the SAEM algorithm
+#> [1] "Tue May 26 18:31:09 2020"
    #> Nonlinear mixed-effects model fit by the SAEM algorithm
 #> -----------------------------------
 #> ----          Data             ----
 #> -----------------------------------
@@ -322,12 +329,12 @@ a list of datasets.
    
 #>             out_values <- out_wide[out_index]
 #>         }
 #>         return(out_values)
-#>     }, mc.cores = 15)
+#>     }, mc.cores = cores)
 #>     res <- unlist(res_list)
 #>     return(res)
 #> }
-#> <bytecode: 0x555559875398>
-#> <environment: 0x55555973a248>
+#> <bytecode: 0x555559668108>
+#> <environment: 0x555559677c08>
 #>   Nb of parameters: 4 
 #>       parameter names:  parent_0 log_k_parent log_k_A1 f_parent_ilr_1 
 #>       distribution:
@@ -353,7 +360,7 @@ a list of datasets.
    
 #>     Estimation of individual parameters (MAP)
 #>     Estimation of standard errors and linearised log-likelihood
 #>     Estimation of log-likelihood by importance sampling
-#>     Number of iterations:  K1=300, K2=100 
+#>     Number of iterations:  K1=200, K2=80 
 #>     Number of chains:  10 
 #>     Seed:  123456 
 #>     Number of MCMC iterations for IS:  5000 
@@ -369,19 +376,19 @@ a list of datasets.
    
 #> -----------------  Fixed effects  ------------------
 #> ----------------------------------------------------
 #>      Parameter      Estimate SE   CV(%)
-#> [1,] parent_0       86.21    1.51  1.7 
+#> [1,] parent_0       86.14    1.61  1.9 
 #> [2,] log_k_parent   -3.21    0.59 18.5 
-#> [3,] log_k_A1       -4.64    0.29  6.3 
-#> [4,] f_parent_ilr_1 -0.32    0.30 93.2 
-#> [5,] a.1             4.69    0.27  5.8 
+#> [3,] log_k_A1       -4.66    0.30  6.4 
+#> [4,] f_parent_ilr_1 -0.33    0.30 91.7 
+#> [5,] a.1             4.68    0.27  5.8 
 #> ----------------------------------------------------
 #> -----------  Variance of random effects  -----------
 #> ----------------------------------------------------
 #>                Parameter             Estimate SE   CV(%)
-#> parent_0       omega2.parent_0       6.07     7.08 117  
-#> log_k_parent   omega2.log_k_parent   1.75     1.11  63  
+#> parent_0       omega2.parent_0       7.71     8.14 106  
+#> log_k_parent   omega2.log_k_parent   1.76     1.12  63  
 #> log_k_A1       omega2.log_k_A1       0.26     0.26 101  
-#> f_parent_ilr_1 omega2.f_parent_ilr_1 0.38     0.27  71  
+#> f_parent_ilr_1 omega2.f_parent_ilr_1 0.39     0.28  71  
 #> ----------------------------------------------------
 #> ------  Correlation matrix of random effects  ------
 #> ----------------------------------------------------
@@ -399,66 +406,16 @@ a list of datasets.
    
 #> ---------------  Statistical criteria  -------------
 #> ----------------------------------------------------
 #> Likelihood computed by linearisation
-#>       -2LL= 1064.397 
-#>       AIC = 1082.397 
-#>       BIC = 1078.882 
+#>       -2LL= 1064.364 
+#>       AIC = 1082.364 
+#>       BIC = 1078.848 
 #> 
 #> Likelihood computed by importance sampling
-#>       -2LL= 1063.161 
-#>       AIC = 1081.161 
-#>       BIC = 1077.646 
-#> ----------------------------------------------------
    #> Nonlinear mixed-effects model fit by the SAEM algorithm
-#> -----------------------------------------
-#> ----         Data and Model          ----
-#> -----------------------------------------
-#> Data
-#>     Dataset ds_saemix 
-#>     Longitudinal data: value ~ time + name | ds 
-#> 
-#> Model:
-#>     Mixed model generated from mmkin object 
-#>      4 parameters: parent_0 log_k_parent log_k_A1 f_parent_ilr_1 
-#>      error model: constant 
-#>      No covariate
-#> 
-#> Key options
-#>     Estimation of individual parameters (MAP)
-#>     Estimation of standard errors and linearised log-likelihood
-#>     Estimation of log-likelihood by importance sampling
-#>     Number of iterations:  K1=300, K2=100 
-#>     Number of chains:  10 
-#>     Seed:  123456 
-#>     Number of MCMC iterations for IS:  5000 
-#>     Input/output
-#>         results not saved
-#>         no graphs
-#> ----------------------------------------------------
-#> ----                  Results                   ----
-#> Fixed effects
-#>  Parameter      Estimate   SE   CV(%)
-#>  parent_0       86.214   1.506  1.75 
-#>  log_k_parent   -3.210   0.593 18.47 
-#>  log_k_A1       -4.643   0.294  6.34 
-#>  f_parent_ilr_1 -0.322   0.300 93.24 
-#>  a.1             4.689   0.270  5.76 
-#> 
-#> Variance of random effects
-#>  Parameter             Estimate   SE   CV(%)
-#>  omega2.parent_0       6.068    7.078 116.7 
-#>  omega2.log_k_parent   1.752    1.111  63.4 
-#>  omega2.log_k_A1       0.256    0.257 100.5 
-#>  omega2.f_parent_ilr_1 0.385    0.273  70.8 
-#> 
-#> Statistical criteria
-#> Likelihood computed by linearisation
-#>       -2LL= 1064.397 
-#>        AIC= 1082.397 
-#>        BIC= 1078.882 
-#> Likelihood computed by importance sampling
-#>       -2LL= 1063.161 
-#>        AIC= 1081.161 
-#>        BIC= 1077.646 
    # }
-
    
+#>       -2LL= 1063.462 
+#>       AIC = 1081.462 
+#>       BIC = 1077.947 
+#> ----------------------------------------------------
    plot(f_saemix, plot.type = "convergence")
    #> Plotting convergence plots
    # }
+
    @@ -439,10 +439,10 @@ print(FOCUS_2006_D)

    A comprehensive report of the results is obtained using the summary method for mkinfit objects.

    summary(fit)
    -
    ## mkin version used for fitting:    0.9.50 
+
    ## mkin version used for fitting:    0.9.50.3 
 ## R version used for fitting:       4.0.0 
-## Date of fit:     Mon May 11 04:41:12 2020 
-## Date of summary: Mon May 11 04:41:12 2020 
+## Date of fit:     Tue May 26 17:01:07 2020 
+## Date of summary: Tue May 26 17:01:07 2020 
 ## 
 ## Equations:
 ## d_parent/dt = - k_parent * parent
diff --git a/vignettes/FOCUS_L.html b/vignettes/FOCUS_L.html
index 968ebf0c..7573ef58 100644
--- a/vignettes/FOCUS_L.html
+++ b/vignettes/FOCUS_L.html
@@ -1,17 +1,17 @@
 
 
-
+
 
 
 
 
-
 
+
 
 
 
 
-
+
 
 Example evaluation of FOCUS Laboratory Data L1 to L3
 
@@ -69,8 +69,6 @@ overflow: auto;
 margin-left: 2%;
 position: fixed;
 border: 1px solid #ccc;
-webkit-border-radius: 6px;
-moz-border-radius: 6px;
 border-radius: 6px;
 }
 
@@ -98,10 +96,15 @@ font-size: 12px;
 .tocify-subheader .tocify-subheader {
 text-indent: 30px;
 }
-
 .tocify-subheader .tocify-subheader .tocify-subheader {
 text-indent: 40px;
 }
+.tocify-subheader .tocify-subheader .tocify-subheader .tocify-subheader {
+text-indent: 50px;
+}
+.tocify-subheader .tocify-subheader .tocify-subheader .tocify-subheader .tocify-subheader {
+text-indent: 60px;
+}
 
 .tocify .tocify-item > a, .tocify .nav-list .nav-header {
 margin: 0px;
@@ -504,13 +507,13 @@ float: none;
 
           item.append($("", {
 
-            "text": self.text()
+            "html": self.html()
 
           }));
 
         } else {
 
-          item.text(self.text());
+          item.html(self.html());
 
         }
 
@@ -1341,7 +1344,6 @@ code {
 }
 img {
   max-width:100%;
-  height: auto;
 }
 .tabbed-pane {
   padding-top: 12px;
@@ -1403,6 +1405,7 @@ summary {
   border: none;
   display: inline-block;
   border-radius: 4px;
+  background-color: transparent;
 }
 
 .tabset-dropdown > .nav-tabs.nav-tabs-open > li {
@@ -1415,49 +1418,10 @@ summary {
 }
 
 
-
-
 
 
 
 
-
-
 
 
 
@@ -217,7 +215,7 @@ code > span.fu { color: #900; font-weight: bold; }  code > span.er { color: #a61
 
 
    Calculation of time weighted average concentrations with mkin
    
 
    Johannes Ranke
    
-
    2019-09-18
    
+
    2020-05-26
    
 
 
 
@@ -261,6 +259,9 @@ code > span.fu { color: #900; font-weight: bold; }  code > span.er { color: #a61
 
 
 
+
+
+