From f820bf5b91be0f589de16c3e3250f5f79672df75 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Fri, 28 Oct 2022 13:39:15 +0200 Subject: Rename parhist to parplot and make it generic That parhist name was not the brightest idea, as it does not show histograms. --- NAMESPACE | 3 +- NEWS.md | 4 +- R/llhist.R | 3 +- R/multistart.R | 6 +- R/parhist.R | 99 ----------- R/parplot.R | 105 +++++++++++ R/saem.R | 2 +- log/test.log | 42 ++--- man/logLik.saem.mmkin.Rd | 2 +- man/multistart.Rd | 6 +- man/parhist.Rd | 43 ----- man/parplot.Rd | 46 +++++ .../multistart/llhist-for-biphasic-saemix-fit.svg | 1 + .../_snaps/multistart/llhist-for-sfo-fit.svg | 2 +- .../multistart/parhist-for-biphasic-saemix-fit.svg | 195 --------------------- .../_snaps/multistart/parhist-for-sfo-fit.svg | 106 ----------- .../multistart/parplot-for-biphasic-saemix-fit.svg | 195 +++++++++++++++++++++ .../_snaps/multistart/parplot-for-sfo-fit.svg | 106 +++++++++++ tests/testthat/test_multistart.R | 8 +- vignettes/web_only/multistart.rmd | 6 +- 20 files changed, 496 insertions(+), 484 deletions(-) delete mode 100644 R/parhist.R create mode 100644 R/parplot.R delete mode 100644 man/parhist.Rd create mode 100644 man/parplot.Rd delete mode 100644 tests/testthat/_snaps/multistart/parhist-for-biphasic-saemix-fit.svg delete mode 100644 tests/testthat/_snaps/multistart/parhist-for-sfo-fit.svg create mode 100644 tests/testthat/_snaps/multistart/parplot-for-biphasic-saemix-fit.svg create mode 100644 tests/testthat/_snaps/multistart/parplot-for-sfo-fit.svg diff --git a/NAMESPACE b/NAMESPACE index a5691eee..4755f631 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -36,6 +36,7 @@ S3method(parms,mkinfit) S3method(parms,mmkin) S3method(parms,multistart) S3method(parms,saem.mmkin) +S3method(parplot,multistart.saem.mmkin) S3method(plot,mixed.mmkin) S3method(plot,mkinfit) S3method(plot,mmkin) @@ -126,8 +127,8 @@ export(nafta) export(nlme) export(nlme_data) export(nlme_function) -export(parhist) export(parms) +export(parplot) export(plot_err) export(plot_res) export(plot_sep) diff --git a/NEWS.md b/NEWS.md index d288c32f..c7927559 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,6 +1,6 @@ # mkin 1.1.2 -- 'R/multistart.R': New method for testing multiple start parameters for hierarchical model fits, with diagnostic plotting functions 'llhist' and 'parhist'. +- 'R/multistart.R': New method for testing multiple start parameters for hierarchical model fits, with diagnostic plotting functions 'llhist' and 'parplot'. - 'R/mhmkin.R': New method for performing multiple hierarchical mkin fits in one function call, optionally in parallel. @@ -8,7 +8,7 @@ - 'R/saem.R': 'logLik' and 'update' methods for 'saem.mmkin' objects. -- 'R/convergence.R': New generic to show convergence information with methods for 'mmkin' and 'mhmkin' objects. +- 'R/status.R': New generic to show status information for fit array objects with methods for 'mmkin' and 'mhmkin' objects. - 'R/illparms.R': New generic to show ill-defined parameters with methods for 'mkinfit', 'mmkin', 'saem.mmkin' and 'mhmkin' objects. diff --git a/R/llhist.R b/R/llhist.R index 22e3aa08..e158495d 100644 --- a/R/llhist.R +++ b/R/llhist.R @@ -25,6 +25,7 @@ llhist <- function(object, breaks = "Sturges", lpos = "topleft", main = "", stop("llhist is only implemented for multistart.saem.mmkin objects") } + ll_orig <- logLik(attr(object, "orig")) ll <- stats::na.omit(sapply(object, llfunc)) par(las = 1) @@ -34,7 +35,7 @@ llhist <- function(object, breaks = "Sturges", lpos = "topleft", main = "", freq_factor <- h$counts[1] / h$density[1] - abline(v = logLik(attr(object, "orig")), col = 2) + abline(v = ll_orig, col = 2) legend(lpos, inset = c(0.05, 0.05), bty = "n", lty = 1, col = c(2), diff --git a/R/multistart.R b/R/multistart.R index 29ccdc44..14683c11 100644 --- a/R/multistart.R +++ b/R/multistart.R @@ -17,7 +17,7 @@ #' @param x The multistart object to print #' @return A list of [saem.mmkin] objects, with class attributes #' 'multistart.saem.mmkin' and 'multistart'. -#' @seealso [parhist], [llhist] +#' @seealso [parplot], [llhist] #' #' @references Duchesne R, Guillemin A, Gandrillon O, Crauste F. Practical #' identifiability in the frame of nonlinear mixed effects models: the example @@ -40,7 +40,7 @@ #' f_mmkin <- mmkin("DFOP", dmta_ds, error_model = "tc", cores = 7, quiet = TRUE) #' f_saem_full <- saem(f_mmkin) #' f_saem_full_multi <- multistart(f_saem_full, n = 16, cores = 16) -#' parhist(f_saem_full_multi, lpos = "bottomleft") +#' parplot(f_saem_full_multi, lpos = "bottomleft") #' illparms(f_saem_full) #' #' f_saem_reduced <- update(f_saem_full, no_random_effect = "log_k2") @@ -52,7 +52,7 @@ #' cl <- makePSOCKcluster(12) #' clusterExport(cl, "f_mmkin") #' f_saem_reduced_multi <- multistart(f_saem_reduced, n = 16, cluster = cl) -#' parhist(f_saem_reduced_multi, lpos = "topright") +#' parplot(f_saem_reduced_multi, lpos = "topright") #' } multistart <- function(object, n = 50, cores = if (Sys.info()["sysname"] == "Windows") 1 else parallel::detectCores(), diff --git a/R/parhist.R b/R/parhist.R deleted file mode 100644 index 0f9c9964..00000000 --- a/R/parhist.R +++ /dev/null @@ -1,99 +0,0 @@ -#' Plot parameter distributions from multistart objects -#' -#' Produces a boxplot with all parameters from the multiple runs, scaled -#' either by the parameters of the run with the highest likelihood, -#' or by their medians as proposed in the paper by Duchesne et al. (2021). -#' -#' @param object The [multistart] object -#' @param llmin The minimum likelihood of objects to be shown -#' @param scale By default, scale parameters using the best available fit. -#' If 'median', parameters are scaled using the median parameters from all fits. -#' @param main Title of the plot -#' @param lpos Positioning of the legend. -#' @param \dots Passed to [boxplot] -#' @references Duchesne R, Guillemin A, Gandrillon O, Crauste F. Practical -#' identifiability in the frame of nonlinear mixed effects models: the example -#' of the in vitro erythropoiesis. BMC Bioinformatics. 2021 Oct 4;22(1):478. -#' doi: 10.1186/s12859-021-04373-4. -#' @seealso [multistart] -#' @importFrom stats median -#' @export -parhist <- function(object, llmin = -Inf, scale = c("best", "median"), - lpos = "bottomleft", main = "", ...) -{ - oldpar <- par(no.readonly = TRUE) - on.exit(par(oldpar, no.readonly = TRUE)) - - orig <- attr(object, "orig") - orig_parms <- parms(orig) - start_parms <- orig$mean_dp_start - all_parms <- parms(object) - - if (inherits(object, "multistart.saem.mmkin")) { - llfunc <- function(object) { - if (inherits(object$so, "try-error")) return(NA) - else return(logLik(object$so)) - } - } else { - stop("parhist is only implemented for multistart.saem.mmkin objects") - } - ll <- sapply(object, llfunc) - selected <- which(ll > llmin) - selected_parms <- all_parms[selected, ] - - par(las = 1) - if (orig$transformations == "mkin") { - degparm_names_transformed <- names(start_parms) - degparm_index <- which(names(orig_parms) %in% degparm_names_transformed) - orig_parms[degparm_names_transformed] <- backtransform_odeparms( - orig_parms[degparm_names_transformed], - orig$mmkin[[1]]$mkinmod, - transform_rates = orig$mmkin[[1]]$transform_rates, - transform_fractions = orig$mmkin[[1]]$transform_fractions) - start_parms <- backtransform_odeparms(start_parms, - orig$mmkin[[1]]$mkinmod, - transform_rates = orig$mmkin[[1]]$transform_rates, - transform_fractions = orig$mmkin[[1]]$transform_fractions) - degparm_names <- names(start_parms) - - names(orig_parms) <- c(degparm_names, names(orig_parms[-degparm_index])) - - selected_parms[, degparm_names_transformed] <- - t(apply(selected_parms[, degparm_names_transformed], 1, backtransform_odeparms, - orig$mmkin[[1]]$mkinmod, - transform_rates = orig$mmkin[[1]]$transform_rates, - transform_fractions = orig$mmkin[[1]]$transform_fractions)) - colnames(selected_parms)[1:length(degparm_names)] <- degparm_names - } - - scale <- match.arg(scale) - parm_scale <- switch(scale, - best = selected_parms[which.best(object[selected]), ], - median = apply(selected_parms, 2, median) - ) - - # Boxplots of all scaled parameters - selected_scaled_parms <- t(apply(selected_parms, 1, function(x) x / parm_scale)) - boxplot(selected_scaled_parms, log = "y", main = main, , - ylab = "Normalised parameters", ...) - - # Show starting parameters - start_scaled_parms <- rep(NA_real_, length(orig_parms)) - names(start_scaled_parms) <- names(orig_parms) - start_scaled_parms[names(start_parms)] <- - start_parms / parm_scale[names(start_parms)] - points(start_scaled_parms, col = 3, cex = 3) - - # Show parameters of original run - orig_scaled_parms <- orig_parms / parm_scale - points(orig_scaled_parms, col = 2, cex = 2) - - abline(h = 1, lty = 2) - - legend(lpos, inset = c(0.05, 0.05), bty = "n", - pch = 1, col = 3:1, lty = c(NA, NA, 1), - legend = c( - "Starting parameters", - "Original run", - "Multistart runs")) -} diff --git a/R/parplot.R b/R/parplot.R new file mode 100644 index 00000000..627a4eb8 --- /dev/null +++ b/R/parplot.R @@ -0,0 +1,105 @@ +#' Plot parameter variability of multistart objects +#' +#' Produces a boxplot with all parameters from the multiple runs, scaled +#' either by the parameters of the run with the highest likelihood, +#' or by their medians as proposed in the paper by Duchesne et al. (2021). +#' +#' @param object The [multistart] object +#' @param llmin The minimum likelihood of objects to be shown +#' @param scale By default, scale parameters using the best available fit. +#' If 'median', parameters are scaled using the median parameters from all fits. +#' @param main Title of the plot +#' @param lpos Positioning of the legend. +#' @param \dots Passed to [boxplot] +#' @references Duchesne R, Guillemin A, Gandrillon O, Crauste F. Practical +#' identifiability in the frame of nonlinear mixed effects models: the example +#' of the in vitro erythropoiesis. BMC Bioinformatics. 2021 Oct 4;22(1):478. +#' doi: 10.1186/s12859-021-04373-4. +#' @seealso [multistart] +#' @importFrom stats median +#' @export +parplot <- function(object, ...) { + UseMethod("parplot") +} + +#' @rdname parplot +#' @export +parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, scale = c("best", "median"), + lpos = "bottomleft", main = "", ...) +{ + oldpar <- par(no.readonly = TRUE) + on.exit(par(oldpar, no.readonly = TRUE)) + + orig <- attr(object, "orig") + orig_parms <- parms(orig) + start_parms <- orig$mean_dp_start + all_parms <- parms(object) + + if (inherits(object, "multistart.saem.mmkin")) { + llfunc <- function(object) { + if (inherits(object$so, "try-error")) return(NA) + else return(logLik(object$so)) + } + } else { + stop("parplot is only implemented for multistart.saem.mmkin objects") + } + ll <- sapply(object, llfunc) + selected <- which(ll > llmin) + selected_parms <- all_parms[selected, ] + + par(las = 1) + if (orig$transformations == "mkin") { + degparm_names_transformed <- names(start_parms) + degparm_index <- which(names(orig_parms) %in% degparm_names_transformed) + orig_parms[degparm_names_transformed] <- backtransform_odeparms( + orig_parms[degparm_names_transformed], + orig$mmkin[[1]]$mkinmod, + transform_rates = orig$mmkin[[1]]$transform_rates, + transform_fractions = orig$mmkin[[1]]$transform_fractions) + start_parms <- backtransform_odeparms(start_parms, + orig$mmkin[[1]]$mkinmod, + transform_rates = orig$mmkin[[1]]$transform_rates, + transform_fractions = orig$mmkin[[1]]$transform_fractions) + degparm_names <- names(start_parms) + + names(orig_parms) <- c(degparm_names, names(orig_parms[-degparm_index])) + + selected_parms[, degparm_names_transformed] <- + t(apply(selected_parms[, degparm_names_transformed], 1, backtransform_odeparms, + orig$mmkin[[1]]$mkinmod, + transform_rates = orig$mmkin[[1]]$transform_rates, + transform_fractions = orig$mmkin[[1]]$transform_fractions)) + colnames(selected_parms)[1:length(degparm_names)] <- degparm_names + } + + scale <- match.arg(scale) + parm_scale <- switch(scale, + best = selected_parms[which.best(object[selected]), ], + median = apply(selected_parms, 2, median) + ) + + # Boxplots of all scaled parameters + selected_scaled_parms <- t(apply(selected_parms, 1, function(x) x / parm_scale)) + boxplot(selected_scaled_parms, log = "y", main = main, , + ylab = "Normalised parameters", ...) + + # Show starting parameters + start_scaled_parms <- rep(NA_real_, length(orig_parms)) + names(start_scaled_parms) <- names(orig_parms) + start_scaled_parms[names(start_parms)] <- + start_parms / parm_scale[names(start_parms)] + points(start_scaled_parms, col = 3, cex = 3) + + # Show parameters of original run + orig_scaled_parms <- orig_parms / parm_scale + points(orig_scaled_parms, col = 2, cex = 2) + + abline(h = 1, lty = 2) + + legend(lpos, inset = c(0.05, 0.05), bty = "n", + pch = 1, col = 3:1, lty = c(NA, NA, 1), + legend = c( + "Starting parameters", + "Original run", + "Multistart runs")) +} diff --git a/R/saem.R b/R/saem.R index cf67b8e1..090ed3bf 100644 --- a/R/saem.R +++ b/R/saem.R @@ -726,7 +726,7 @@ saemix_data <- function(object, covariates = NULL, verbose = FALSE, ...) { #' @param \dots Passed to [saemix::logLik.SaemixObject] #' @param method Passed to [saemix::logLik.SaemixObject] #' @export -logLik.saem.mmkin <- function(object, ..., method = c("lin", "is", "gq")) { +logLik.saem.mmkin <- function(object, ..., method = c("is", "lin", "gq")) { method <- match.arg(method) return(logLik(object$so, method = method)) } diff --git a/log/test.log b/log/test.log index 429a2e02..0e2ca6b2 100644 --- a/log/test.log +++ b/log/test.log @@ -1,53 +1,53 @@ ℹ Testing mkin ✔ | F W S OK | Context ✔ | 5 | AIC calculation -✔ | 5 | Analytical solutions for coupled models [3.3s] +✔ | 5 | Analytical solutions for coupled models [3.4s] ✔ | 5 | Calculation of Akaike weights ✔ | 3 | Export dataset for reading into CAKE ✔ | 12 | Confidence intervals and p-values [1.0s] -✔ | 1 12 | Dimethenamid data from 2018 [43.2s] +✔ | 1 12 | Dimethenamid data from 2018 [63.1s] ──────────────────────────────────────────────────────────────────────────────── Skip (test_dmta.R:98:3): Different backends get consistent results for SFO-SFO3+, dimethenamid data Reason: Fitting this ODE model with saemix takes about 15 minutes on my system ──────────────────────────────────────────────────────────────────────────────── -✔ | 14 | Error model fitting [11.5s] +✔ | 14 | Error model fitting [5.5s] ✔ | 5 | Time step normalisation -✔ | 4 | Calculation of FOCUS chi2 error levels [1.3s] -✔ | 14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [1.9s] -✔ | 4 | Test fitting the decline of metabolites from their maximum [0.9s] -✔ | 1 | Fitting the logistic model [0.6s] -✔ | 7 | Batch fitting and diagnosing hierarchical kinetic models [30.9s] -✔ | 1 12 | Nonlinear mixed-effects models [0.4s] +✔ | 4 | Calculation of FOCUS chi2 error levels [0.6s] +✔ | 14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [0.8s] +✔ | 4 | Test fitting the decline of metabolites from their maximum [0.4s] +✔ | 1 | Fitting the logistic model [0.2s] +✔ | 7 | Batch fitting and diagnosing hierarchical kinetic models [28.6s] +✔ | 1 12 | Nonlinear mixed-effects models [0.6s] ──────────────────────────────────────────────────────────────────────────────── Skip (test_mixed.R:74:3): saemix results are reproducible for biphasic fits Reason: Fitting with saemix takes around 10 minutes when using deSolve ──────────────────────────────────────────────────────────────────────────────── ✔ | 3 | Test dataset classes mkinds and mkindsg -✔ | 10 | Special cases of mkinfit calls [0.7s] -✔ | 3 | mkinfit features [1.0s] -✔ | 8 | mkinmod model generation and printing [0.3s] -✔ | 3 | Model predictions with mkinpredict [0.3s] -✔ | 7 | Multistart method for saem.mmkin models [39.8s] -✔ | 16 | Evaluations according to 2015 NAFTA guidance [2.6s] -✔ | 9 | Nonlinear mixed-effects models with nlme [8.8s] -✔ | 16 | Plotting [10.1s] +✔ | 10 | Special cases of mkinfit calls [1.1s] +✔ | 3 | mkinfit features [1.7s] +✔ | 8 | mkinmod model generation and printing [0.4s] +✔ | 3 | Model predictions with mkinpredict [0.6s] +✔ | 7 | Multistart method for saem.mmkin models [93.9s] +✔ | 16 | Evaluations according to 2015 NAFTA guidance [7.0s] +✔ | 9 | Nonlinear mixed-effects models with nlme [13.0s] +✔ | 16 | Plotting [11.3s] ✔ | 4 | Residuals extracted from mkinfit models -✔ | 1 36 | saemix parent models [72.3s] +✔ | 1 36 | saemix parent models [72.6s] ──────────────────────────────────────────────────────────────────────────────── Skip (test_saemix_parent.R:152:3): We can also use mkin solution methods for saem Reason: This still takes almost 2.5 minutes although we do not solve ODEs ──────────────────────────────────────────────────────────────────────────────── ✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.4s] ✔ | 11 | Processing of residue series -✔ | 7 | Fitting the SFORB model [3.8s] +✔ | 7 | Fitting the SFORB model [3.7s] ✔ | 1 | Summaries of old mkinfit objects ✔ | 5 | Summary [0.2s] ✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.2s] -✔ | 9 | Hypothesis tests [8.0s] +✔ | 9 | Hypothesis tests [8.1s] ✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s] ══ Results ═════════════════════════════════════════════════════════════════════ -Duration: 249.3 s +Duration: 324.3 s ── Skipped tests ────────────────────────────────────────────────────────────── • Fitting this ODE model with saemix takes about 15 minutes on my system (1) diff --git a/man/logLik.saem.mmkin.Rd b/man/logLik.saem.mmkin.Rd index 603f4607..bd0bb72e 100644 --- a/man/logLik.saem.mmkin.Rd +++ b/man/logLik.saem.mmkin.Rd @@ -4,7 +4,7 @@ \alias{logLik.saem.mmkin} \title{logLik method for saem.mmkin objects} \usage{ -\method{logLik}{saem.mmkin}(object, ..., method = c("lin", "is", "gq")) +\method{logLik}{saem.mmkin}(object, ..., method = c("is", "lin", "gq")) } \arguments{ \item{object}{The fitted \link{saem.mmkin} object} diff --git a/man/multistart.Rd b/man/multistart.Rd index ebcc7d80..1f6773fe 100644 --- a/man/multistart.Rd +++ b/man/multistart.Rd @@ -77,7 +77,7 @@ dmta_ds[["Elliot 1"]] <- dmta_ds[["Elliot 2"]] <- NULL f_mmkin <- mmkin("DFOP", dmta_ds, error_model = "tc", cores = 7, quiet = TRUE) f_saem_full <- saem(f_mmkin) f_saem_full_multi <- multistart(f_saem_full, n = 16, cores = 16) -parhist(f_saem_full_multi, lpos = "bottomleft") +parplot(f_saem_full_multi, lpos = "bottomleft") illparms(f_saem_full) f_saem_reduced <- update(f_saem_full, no_random_effect = "log_k2") @@ -89,7 +89,7 @@ library(parallel) cl <- makePSOCKcluster(12) clusterExport(cl, "f_mmkin") f_saem_reduced_multi <- multistart(f_saem_reduced, n = 16, cluster = cl) -parhist(f_saem_reduced_multi, lpos = "topright") +parplot(f_saem_reduced_multi, lpos = "topright") } } \references{ @@ -99,5 +99,5 @@ of the in vitro erythropoiesis. BMC Bioinformatics. 2021 Oct 4;22(1):478. doi: 10.1186/s12859-021-04373-4. } \seealso{ -\link{parhist}, \link{llhist} +\link{parplot}, \link{llhist} } diff --git a/man/parhist.Rd b/man/parhist.Rd deleted file mode 100644 index f86ff734..00000000 --- a/man/parhist.Rd +++ /dev/null @@ -1,43 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/parhist.R -\name{parhist} -\alias{parhist} -\title{Plot parameter distributions from multistart objects} -\usage{ -parhist( - object, - llmin = -Inf, - scale = c("best", "median"), - lpos = "bottomleft", - main = "", - ... -) -} -\arguments{ -\item{object}{The \link{multistart} object} - -\item{llmin}{The minimum likelihood of objects to be shown} - -\item{scale}{By default, scale parameters using the best available fit. -If 'median', parameters are scaled using the median parameters from all fits.} - -\item{lpos}{Positioning of the legend.} - -\item{main}{Title of the plot} - -\item{\dots}{Passed to \link{boxplot}} -} -\description{ -Produces a boxplot with all parameters from the multiple runs, scaled -either by the parameters of the run with the highest likelihood, -or by their medians as proposed in the paper by Duchesne et al. (2021). -} -\references{ -Duchesne R, Guillemin A, Gandrillon O, Crauste F. Practical -identifiability in the frame of nonlinear mixed effects models: the example -of the in vitro erythropoiesis. BMC Bioinformatics. 2021 Oct 4;22(1):478. -doi: 10.1186/s12859-021-04373-4. -} -\seealso{ -\link{multistart} -} diff --git a/man/parplot.Rd b/man/parplot.Rd new file mode 100644 index 00000000..37c5841d --- /dev/null +++ b/man/parplot.Rd @@ -0,0 +1,46 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/parplot.R +\name{parplot} +\alias{parplot} +\alias{parplot.multistart.saem.mmkin} +\title{Plot parameter variability of multistart objects} +\usage{ +parplot(object, ...) + +\method{parplot}{multistart.saem.mmkin}( + object, + llmin = -Inf, + scale = c("best", "median"), + lpos = "bottomleft", + main = "", + ... +) +} +\arguments{ +\item{object}{The \link{multistart} object} + +\item{\dots}{Passed to \link{boxplot}} + +\item{llmin}{The minimum likelihood of objects to be shown} + +\item{scale}{By default, scale parameters using the best available fit. +If 'median', parameters are scaled using the median parameters from all fits.} + +\item{lpos}{Positioning of the legend.} + +\item{main}{Title of the plot} +} +\description{ +Produces a boxplot with all parameters from the multiple runs, scaled +either by the parameters of the run with the highest likelihood, +or by their medians as proposed in the paper by Duchesne et al. (2021). +} +\references{ +Duchesne R, Guillemin A, Gandrillon O, Crauste F. Practical +identifiability in the frame of nonlinear mixed effects models: the example +of the in vitro erythropoiesis. BMC Bioinformatics. 2021 Oct 4;22(1):478. +doi: 10.1186/s12859-021-04373-4. +} +\seealso{ +\link{multistart} +} diff --git a/tests/testthat/_snaps/multistart/llhist-for-biphasic-saemix-fit.svg b/tests/testthat/_snaps/multistart/llhist-for-biphasic-saemix-fit.svg index fa92c92a..fe38865d 100644 --- a/tests/testthat/_snaps/multistart/llhist-for-biphasic-saemix-fit.svg +++ b/tests/testthat/_snaps/multistart/llhist-for-biphasic-saemix-fit.svg @@ -50,6 +50,7 @@ + original fit diff --git a/tests/testthat/_snaps/multistart/llhist-for-sfo-fit.svg b/tests/testthat/_snaps/multistart/llhist-for-sfo-fit.svg index 4342df9b..98513d06 100644 --- a/tests/testthat/_snaps/multistart/llhist-for-sfo-fit.svg +++ b/tests/testthat/_snaps/multistart/llhist-for-sfo-fit.svg @@ -50,7 +50,7 @@ - + original fit diff --git a/tests/testthat/_snaps/multistart/parhist-for-biphasic-saemix-fit.svg b/tests/testthat/_snaps/multistart/parhist-for-biphasic-saemix-fit.svg deleted file mode 100644 index 75519b07..00000000 --- a/tests/testthat/_snaps/multistart/parhist-for-biphasic-saemix-fit.svg +++ /dev/null @@ -1,195 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -parent_0 -f_parent_to_m1 -k2 -g -a.1 -b.1 -SD.log_k_m1 -SD.log_k1 -SD.g_qlogis - - - - - -0.5 -1.0 -1.5 -2.0 -Normalised parameters - - - - - - - - - - - - - - - - - - - - - - - - - - - - -Starting parameters -Original run -Multistart runs - - diff --git a/tests/testthat/_snaps/multistart/parhist-for-sfo-fit.svg b/tests/testthat/_snaps/multistart/parhist-for-sfo-fit.svg deleted file mode 100644 index f3373901..00000000 --- a/tests/testthat/_snaps/multistart/parhist-for-sfo-fit.svg +++ /dev/null @@ -1,106 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -parent_0 -k_parent -a.1 -b.1 -SD.k_parent - - - - - -0.5 -1.0 -1.5 -2.0 -Normalised parameters - - - - - - - - - - - - - - - -Starting parameters -Original run -Multistart runs - - diff --git a/tests/testthat/_snaps/multistart/parplot-for-biphasic-saemix-fit.svg b/tests/testthat/_snaps/multistart/parplot-for-biphasic-saemix-fit.svg new file mode 100644 index 00000000..75519b07 --- /dev/null +++ b/tests/testthat/_snaps/multistart/parplot-for-biphasic-saemix-fit.svg @@ -0,0 +1,195 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +parent_0 +f_parent_to_m1 +k2 +g +a.1 +b.1 +SD.log_k_m1 +SD.log_k1 +SD.g_qlogis + + + + + +0.5 +1.0 +1.5 +2.0 +Normalised parameters + + + + + + + + + + + + + + + + + + + + + + + + + + + + +Starting parameters +Original run +Multistart runs + + diff --git a/tests/testthat/_snaps/multistart/parplot-for-sfo-fit.svg b/tests/testthat/_snaps/multistart/parplot-for-sfo-fit.svg new file mode 100644 index 00000000..f3373901 --- /dev/null +++ b/tests/testthat/_snaps/multistart/parplot-for-sfo-fit.svg @@ -0,0 +1,106 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +parent_0 +k_parent +a.1 +b.1 +SD.k_parent + + + + + +0.5 +1.0 +1.5 +2.0 +Normalised parameters + + + + + + + + + + + + + + + +Starting parameters +Original run +Multistart runs + + diff --git a/tests/testthat/test_multistart.R b/tests/testthat/test_multistart.R index e4885191..56eb140c 100644 --- a/tests/testthat/test_multistart.R +++ b/tests/testthat/test_multistart.R @@ -25,14 +25,14 @@ test_that("multistart works for saem.mmkin models", { skip_on_travis() # Plots are platform dependent llhist_sfo <- function() llhist(saem_sfo_s_multi) - parhist_sfo <- function() parhist(saem_sfo_s_multi, ylim = c(0.5, 2)) + parplot_sfo <- function() parplot(saem_sfo_s_multi, ylim = c(0.5, 2)) vdiffr::expect_doppelganger("llhist for sfo fit", llhist_sfo) - vdiffr::expect_doppelganger("parhist for sfo fit", parhist_sfo) + vdiffr::expect_doppelganger("parplot for sfo fit", parplot_sfo) llhist_biphasic <- function() llhist(saem_biphasic_m_multi) - parhist_biphasic <- function() parhist(saem_biphasic_m_multi, + parplot_biphasic <- function() parplot(saem_biphasic_m_multi, ylim = c(0.5, 2)) vdiffr::expect_doppelganger("llhist for biphasic saemix fit", llhist_biphasic) - vdiffr::expect_doppelganger("parhist for biphasic saemix fit", parhist_biphasic) + vdiffr::expect_doppelganger("parplot for biphasic saemix fit", parplot_biphasic) }) diff --git a/vignettes/web_only/multistart.rmd b/vignettes/web_only/multistart.rmd index 94f4ef98..27a8a96a 100644 --- a/vignettes/web_only/multistart.rmd +++ b/vignettes/web_only/multistart.rmd @@ -40,7 +40,7 @@ with different starting values. ```{r} f_saem_full_multi <- multistart(f_saem_full, n = 16, cores = 16) -parhist(f_saem_full_multi) +parplot(f_saem_full_multi) ``` This confirms that the variance of k2 is the most problematic parameter, so we @@ -51,7 +51,7 @@ for k2. f_saem_reduced <- update(f_saem_full, no_random_effect = "log_k2") illparms(f_saem_reduced) f_saem_reduced_multi <- multistart(f_saem_reduced, n = 16, cores = 16) -parhist(f_saem_reduced_multi, lpos = "topright") +parplot(f_saem_reduced_multi, lpos = "topright") ``` The results confirm that all remaining parameters can be determined with sufficient @@ -67,7 +67,7 @@ The parameter histograms can be further improved by excluding the result with the low likelihood. ```{r} -parhist(f_saem_reduced_multi, lpos = "topright", llmin = -326, ylim = c(0.5, 2)) +parplot(f_saem_reduced_multi, lpos = "topright", llmin = -326, ylim = c(0.5, 2)) ``` We can use the `anova` method to compare the models, including a likelihood ratio -- cgit v1.2.1