From 06df20703a9390692ab1ece3ae1702a71fff05ae Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 26 Oct 2022 08:40:58 +0200 Subject: Remove kernel density estimate from llhist --- NAMESPACE | 1 - R/llhist.R | 20 +++---- log/test.log | 30 +++++----- man/llhist.Rd | 5 +- .../multistart/llhist-for-biphasic-saemix-fit.svg | 33 ++++++----- .../_snaps/multistart/llhist-for-sfo-fit.svg | 65 ++++++++++++++++++++++ tests/testthat/test_multistart.R | 5 ++ 7 files changed, 112 insertions(+), 47 deletions(-) create mode 100644 tests/testthat/_snaps/multistart/llhist-for-sfo-fit.svg diff --git a/NAMESPACE b/NAMESPACE index aaaa6484..8ffa0425 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -143,7 +143,6 @@ export(which.best) import(deSolve) import(graphics) import(nlme) -importFrom(KernSmooth,bkde) importFrom(R6,R6Class) importFrom(grDevices,dev.cur) importFrom(lmtest,lrtest) diff --git a/R/llhist.R b/R/llhist.R index 9ddf5b10..22e3aa08 100644 --- a/R/llhist.R +++ b/R/llhist.R @@ -1,8 +1,7 @@ #' Plot the distribution of log likelihoods from multistart objects #' -#' Produces a histogram of log-likelihoods, and an overlayed kernel density -#' estimate. In addition, the likelihood of the original fit is shown as -#' a red vertical line. +#' Produces a histogram of log-likelihoods. In addition, the likelihood of the +#' original fit is shown as a red vertical line. #' #' @param object The [multistart] object #' @param breaks Passed to [hist] @@ -10,9 +9,10 @@ #' @param main Title of the plot #' @param \dots Passed to [hist] #' @seealso [multistart] -#' @importFrom KernSmooth bkde #' @export -llhist <- function(object, breaks = "Sturges", lpos = "topleft", main = "", ...) { +llhist <- function(object, breaks = "Sturges", lpos = "topleft", main = "", + ...) +{ oldpar <- par(no.readonly = TRUE) on.exit(par(oldpar, no.readonly = TRUE)) @@ -26,19 +26,17 @@ llhist <- function(object, breaks = "Sturges", lpos = "topleft", main = "", ...) } ll <- stats::na.omit(sapply(object, llfunc)) - kde <- KernSmooth::bkde(ll) par(las = 1) h <- hist(ll, freq = TRUE, - xlim = range(kde$x), xlab = "", main = main, ylab = "Frequency of log likelihoods", breaks = breaks, ...) freq_factor <- h$counts[1] / h$density[1] - lines(kde$x, freq_factor * kde$y) + abline(v = logLik(attr(object, "orig")), col = 2) + legend(lpos, inset = c(0.05, 0.05), bty = "n", - lty = 1, col = c(2, 1), - legend = c("original log likelihood", - "kernel density estimate")) + lty = 1, col = c(2), + legend = "original fit") } diff --git a/log/test.log b/log/test.log index 4515f383..c5bc0e44 100644 --- a/log/test.log +++ b/log/test.log @@ -1,22 +1,22 @@ ℹ Testing mkin ✔ | F W S OK | Context ✔ | 5 | AIC calculation -✔ | 5 | Analytical solutions for coupled models [3.4s] +✔ | 5 | Analytical solutions for coupled models [3.2s] ✔ | 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 [32.3s] +✔ | 1 12 | Dimethenamid data from 2018 [32.0s] ──────────────────────────────────────────────────────────────────────────────── 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 [5.0s] +✔ | 14 | Error model fitting [4.9s] ✔ | 5 | Time step normalisation ✔ | 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] +✔ | 4 | Test fitting the decline of metabolites from their maximum [0.3s] ✔ | 1 | Fitting the logistic model [0.2s] -✔ | 7 | Batch fitting and diagnosing hierarchical kinetic models [14.3s] +✔ | 7 | Batch fitting and diagnosing hierarchical kinetic models [14.1s] ✔ | 1 12 | Nonlinear mixed-effects models [0.3s] ──────────────────────────────────────────────────────────────────────────────── Skip (test_mixed.R:74:3): saemix results are reproducible for biphasic fits @@ -27,31 +27,31 @@ Reason: Fitting with saemix takes around 10 minutes when using deSolve ✔ | 3 | mkinfit features [0.7s] ✔ | 8 | mkinmod model generation and printing [0.2s] ✔ | 3 | Model predictions with mkinpredict [0.3s] -✔ | 3 | Multistart method for saem.mmkin models [18.4s] -✔ | 16 | Evaluations according to 2015 NAFTA guidance [1.8s] +✔ | 4 | Multistart method for saem.mmkin models [33.5s] +✔ | 16 | Evaluations according to 2015 NAFTA guidance [2.3s] ✔ | 9 | Nonlinear mixed-effects models with nlme [8.5s] -✔ | 16 | Plotting [10.2s] +✔ | 16 | Plotting [9.9s] ✔ | 4 | Residuals extracted from mkinfit models -✔ | 1 35 | saemix parent models [72.5s] +✔ | 1 35 | saemix parent models [71.4s] ──────────────────────────────────────────────────────────────────────────────── Skip (test_saemix_parent.R:143: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] +✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.5s] ✔ | 11 | Processing of residue series ✔ | 7 | Fitting the SFORB model [3.8s] ✔ | 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] -✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s] +✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.3s] +✔ | 9 | Hypothesis tests [8.2s] +✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.1s] ══ Results ═════════════════════════════════════════════════════════════════════ -Duration: 189.5 s +Duration: 203.0 s ── Skipped tests ────────────────────────────────────────────────────────────── • Fitting this ODE model with saemix takes about 15 minutes on my system (1) • Fitting with saemix takes around 10 minutes when using deSolve (1) • This still takes almost 2.5 minutes although we do not solve ODEs (1) -[ FAIL 0 | WARN 0 | SKIP 3 | PASS 256 ] +[ FAIL 0 | WARN 0 | SKIP 3 | PASS 257 ] diff --git a/man/llhist.Rd b/man/llhist.Rd index 09dc85d4..bec30ecf 100644 --- a/man/llhist.Rd +++ b/man/llhist.Rd @@ -18,9 +18,8 @@ llhist(object, breaks = "Sturges", lpos = "topleft", main = "", ...) \item{\dots}{Passed to \link{hist}} } \description{ -Produces a histogram of log-likelihoods, and an overlayed kernel density -estimate. In addition, the likelihood of the original fit is shown as -a red vertical line. +Produces a histogram of log-likelihoods. In addition, the likelihood of the +original fit is shown as a red vertical line. } \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 4767e229..fa92c92a 100644 --- a/tests/testthat/_snaps/multistart/llhist-for-biphasic-saemix-fit.svg +++ b/tests/testthat/_snaps/multistart/llhist-for-biphasic-saemix-fit.svg @@ -19,15 +19,17 @@ Frequency of log likelihoods - - - - - --1190 --1180 --1170 --1160 + + + + + + +-1185 +-1180 +-1175 +-1170 +-1165 @@ -44,14 +46,11 @@ - - - - - + + + + - -original log likelihood -kernel density estimate +original fit diff --git a/tests/testthat/_snaps/multistart/llhist-for-sfo-fit.svg b/tests/testthat/_snaps/multistart/llhist-for-sfo-fit.svg new file mode 100644 index 00000000..21fa8af8 --- /dev/null +++ b/tests/testthat/_snaps/multistart/llhist-for-sfo-fit.svg @@ -0,0 +1,65 @@ + + + + + + + + + + + + +Frequency of log likelihoods + + + + + + + + +-644.0 +-644.5 +-645.0 +-645.5 +-646.0 +-646.5 +-647.0 + + + + + + +0 +1 +2 +3 +4 + + + + + + + + + + + + + + + +original fit + + diff --git a/tests/testthat/test_multistart.R b/tests/testthat/test_multistart.R index 2c7a2bc0..60fdefdc 100644 --- a/tests/testthat/test_multistart.R +++ b/tests/testthat/test_multistart.R @@ -1,6 +1,11 @@ context("Multistart method for saem.mmkin models") test_that("multistart works for saem.mmkin models", { + saem_sfo_s_multi <- multistart(sfo_saem_1, n = 8, cores = n_cores) + + llhist_sfo <- function() llhist(saem_sfo_s_multi, xlim = c(-644, -647)) + vdiffr::expect_doppelganger("llhist for sfo fit", llhist_sfo) + set.seed(123456) saem_biphasic_m_multi <- multistart(saem_biphasic_m, n = 8, cores = n_cores) -- cgit v1.2.1