From 21ad91256dc29423bd905de5c298fd23862b1f3b Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Mon, 14 Nov 2022 20:03:42 +0100 Subject: Automatic starting parameters for saem.mmkin For the case of mkin transformations. This gives faster convergence, and appears to avoid problems with numeric ODE solutions --- GNUmakefile | 2 +- NEWS.md | 2 + R/saem.R | 29 +- log/check.log | 4 +- log/test.log | 32 +- man/saem.Rd | 16 +- .../multistart/llhist-for-biphasic-saemix-fit.svg | 41 +- .../multistart/parplot-for-biphasic-saemix-fit.svg | 242 +- ...t-for-saem-object-with-mkin-transformations.svg | 2332 ++++++++++---------- tests/testthat/anova_sfo_saem.txt | 2 +- tests/testthat/print_dfop_saemix_1.txt | 16 +- tests/testthat/summary_hfit_sfo_tc.txt | 6 +- tests/testthat/test_multistart.R | 5 +- tests/testthat/test_saemix_parent.R | 40 +- 14 files changed, 1407 insertions(+), 1362 deletions(-) diff --git a/GNUmakefile b/GNUmakefile index 114f7eef..6e75d666 100644 --- a/GNUmakefile +++ b/GNUmakefile @@ -103,7 +103,7 @@ vignettes/%.html: vignettes/mkin_vignettes.css vignettes/references.bib vignette vignettes: vignettes/mkin.html vignettes/FOCUS_D.html vignettes/FOCUS_L.html vignettes/twa.html vignettes/web_only/%.html: vignettes/references.bib vignettes/web_only/%.rmd - "$(RBIN)/Rscript" -e "tools::buildVignette(file = 'vignettes/web_only/$*.rmd', dir = 'vignettes/web_only', keep='mkin_benchmarks.rda')" + "$(RBIN)/Rscript" -e "tools::buildVignette(file = 'vignettes/web_only/$*.rmd', dir = 'vignettes/web_only', keep=c('mkin_benchmarks.rda', 'saem_benchmarks.rda'))" articles: vignettes/web_only/FOCUS_Z.html vignettes/web_only/compiled_models.html vignettes/web_only/benchmarks.html vignettes/web_only/dimethenamid_2018.html vignettes/web_only/multistart.html diff --git a/NEWS.md b/NEWS.md index ebd9a9f5..846c7c50 100644 --- a/NEWS.md +++ b/NEWS.md @@ -10,6 +10,8 @@ - 'R/saem.R': 'logLik', 'update' and 'anova' methods for 'saem.mmkin' objects. +- 'R/saem.R': Automatic estimation of start parameters for random effects for the case of mkin transformations, nicely improving convergence and reducing problems with iterative ODE solutions. + - 'R/status.R': New generic to show status information for fit array objects with methods for 'mmkin', 'mhmkin' and 'multistart' objects. - 'R/summary.mmkin.R': Summary method for mmkin objects. diff --git a/R/saem.R b/R/saem.R index 87718603..e980ed1f 100644 --- a/R/saem.R +++ b/R/saem.R @@ -37,14 +37,22 @@ utils::globalVariables(c("predicted", "std")) #' @param no_random_effect Character vector of degradation parameters for #' which there should be no variability over the groups. Only used #' if the covariance model is not explicitly specified. -#' @param covariance.model Will be passed to [saemix::SaemixModel()]. Per +#' @param covariance.model Will be passed to [saemix::saemixModel()]. Per #' default, uncorrelated random effects are specified for all degradation #' parameters. +#' @param omega.init Will be passed to [saemix::saemixModel()]. If using +#' mkin transformations and the default covariance model with optionally +#' excluded random effects, the variances of the degradation parameters +#' are estimated using [mean_degparms], with testing of untransformed +#' log parameters for significant difference from zero. If not using +#' mkin transformations or a custom covariance model, the default +#' initialisation of [saemix::saemixModel] is used for omega.init. #' @param covariates A data frame with covariate data for use in #' 'covariate_models', with dataset names as row names. #' @param covariate_models A list containing linear model formulas with one explanatory #' variable, i.e. of the type 'parameter ~ covariate'. Covariates must be available #' in the 'covariates' data frame. +#' @param error.init Will be passed to [saemix::saemixModel()]. #' @param quiet Should we suppress the messages saemix prints at the beginning #' and the end of the optimisation process? #' @param nbiter.saemix Convenience option to increase the number of @@ -136,9 +144,11 @@ saem.mmkin <- function(object, conf.level = 0.6, solution_type = "auto", covariance.model = "auto", + omega.init = "auto", covariates = NULL, covariate_models = NULL, no_random_effect = NULL, + error.init = c(3, 0.1), nbiter.saemix = c(300, 100), control = list(displayProgress = FALSE, print = FALSE, nbiter.saemix = nbiter.saemix, @@ -154,8 +164,10 @@ saem.mmkin <- function(object, solution_type = solution_type, transformations = transformations, covariance.model = covariance.model, + omega.init = omega.init, covariates = covariates, covariate_models = covariate_models, + error.init = error.init, no_random_effect = no_random_effect, ...) d_saemix <- saemix_data(object, covariates = covariates, verbose = verbose) @@ -283,7 +295,9 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mkin", "saemix"), error_model = "auto", degparms_start = numeric(), covariance.model = "auto", no_random_effect = NULL, + omega.init = "auto", covariates = NULL, covariate_models = NULL, + error.init = numeric(), test_log_parms = FALSE, conf.level = 0.6, verbose = FALSE, ...) { if (nrow(object) > 1) stop("Only row objects allowed") @@ -637,6 +651,18 @@ saemix_model <- function(object, solution_type = "auto", covariance.model = diag(covariance_diagonal) } + if (omega.init[1] == "auto") { + if (transformations == "mkin") { + degparms_eta_ini <- as.numeric( # remove names + mean_degparms(object, + random = TRUE, test_log_parms = TRUE)$eta) + + omega.init <- 2 * diag(degparms_eta_ini^2) + } else { + omega.init <- matrix(nrow = 0, ncol = 0) + } + } + if (is.null(covariate_models)) { covariate.model <- matrix(nrow = 0, ncol = 0) # default in saemixModel() } else { @@ -680,6 +706,7 @@ saemix_model <- function(object, solution_type = "auto", verbose = verbose, covariance.model = covariance.model, covariate.model = covariate.model, + omega.init = omega.init, ... ) attr(res, "mean_dp_start") <- degparms_optim diff --git a/log/check.log b/log/check.log index d8225158..3fea2ec6 100644 --- a/log/check.log +++ b/log/check.log @@ -57,12 +57,12 @@ Maintainer: ‘Johannes Ranke ’ * checking data for ASCII and uncompressed saves ... OK * checking installed files from ‘inst/doc’ ... OK * checking files in ‘vignettes’ ... OK -* checking examples ... [36s/46s] OK +* checking examples ... [15s/15s] OK * checking for unstated dependencies in ‘tests’ ... OK * checking tests ... SKIPPED * checking for unstated dependencies in vignettes ... OK * checking package vignettes in ‘inst/doc’ ... OK -* checking re-building of vignette outputs ... [26s/41s] OK +* checking re-building of vignette outputs ... OK * checking PDF version of manual ... OK * checking HTML version of manual ... OK * checking for non-standard things in the check directory ... OK diff --git a/log/test.log b/log/test.log index 874d6e24..f20a7f18 100644 --- a/log/test.log +++ b/log/test.log @@ -1,45 +1,45 @@ ℹ 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.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 [31.5s] +✔ | 1 12 | Dimethenamid data from 2018 [32.1s] ──────────────────────────────────────────────────────────────────────────────── -Skip (test_dmta.R:98:3): Different backends get consistent results for SFO-SFO3+, dimethenamid data +Skip ('test_dmta.R:98'): 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] ✔ | 10 | Batch fitting and diagnosing hierarchical kinetic models [24.0s] ✔ | 1 12 | Nonlinear mixed-effects models [0.3s] ──────────────────────────────────────────────────────────────────────────────── -Skip (test_mixed.R:74:3): saemix results are reproducible for biphasic fits +Skip ('test_mixed.R:74'): 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.5s] ✔ | 3 | mkinfit features [0.7s] ✔ | 8 | mkinmod model generation and printing [0.2s] -✔ | 3 | Model predictions with mkinpredict [0.3s] -✔ | 7 | Multistart method for saem.mmkin models [36.3s] -✔ | 16 | Evaluations according to 2015 NAFTA guidance [2.4s] -✔ | 9 | Nonlinear mixed-effects models with nlme [8.8s] -✔ | 16 | Plotting [10.1s] +✔ | 3 | Model predictions with mkinpredict [0.4s] +✔ | 7 | Multistart method for saem.mmkin models [36.7s] +✔ | 16 | Evaluations according to 2015 NAFTA guidance [2.3s] +✔ | 9 | Nonlinear mixed-effects models with nlme [8.5s] +✔ | 16 | Plotting [9.8s] ✔ | 4 | Residuals extracted from mkinfit models -✔ | 1 37 | saemix parent models [71.5s] +✔ | 1 36 | saemix parent models [66.4s] ──────────────────────────────────────────────────────────────────────────────── -Skip (test_saemix_parent.R:153:3): We can also use mkin solution methods for saem +Skip ('test_saemix_parent.R:143'): 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 -✔ | 10 | Fitting the SFORB model [3.8s] +✔ | 10 | 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] @@ -47,11 +47,11 @@ Reason: This still takes almost 2.5 minutes although we do not solve ODEs ✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s] ══ Results ═════════════════════════════════════════════════════════════════════ -Duration: 216.0 s +Duration: 211.1 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 268 ] +[ FAIL 0 | WARN 0 | SKIP 3 | PASS 267 ] diff --git a/man/saem.Rd b/man/saem.Rd index 3dc8001a..11463351 100644 --- a/man/saem.Rd +++ b/man/saem.Rd @@ -20,9 +20,11 @@ saem(object, ...) conf.level = 0.6, solution_type = "auto", covariance.model = "auto", + omega.init = "auto", covariates = NULL, covariate_models = NULL, no_random_effect = NULL, + error.init = c(3, 0.1), nbiter.saemix = c(300, 100), control = list(displayProgress = FALSE, print = FALSE, nbiter.saemix = nbiter.saemix, save = FALSE, save.graphs = FALSE), @@ -41,8 +43,10 @@ saemix_model( degparms_start = numeric(), covariance.model = "auto", no_random_effect = NULL, + omega.init = "auto", covariates = NULL, covariate_models = NULL, + error.init = numeric(), test_log_parms = FALSE, conf.level = 0.6, verbose = FALSE, @@ -81,10 +85,18 @@ for parameter that are tested if requested by 'test_log_parms'.} \item{solution_type}{Possibility to specify the solution type in case the automatic choice is not desired} -\item{covariance.model}{Will be passed to \code{\link[saemix:SaemixModel-class]{saemix::SaemixModel()}}. Per +\item{covariance.model}{Will be passed to \code{\link[saemix:saemixModel]{saemix::saemixModel()}}. Per default, uncorrelated random effects are specified for all degradation parameters.} +\item{omega.init}{Will be passed to \code{\link[saemix:saemixModel]{saemix::saemixModel()}}. If using +mkin transformations and the default covariance model with optionally +excluded random effects, the variances of the degradation parameters +are estimated using \link{mean_degparms}, with testing of untransformed +log parameters for significant difference from zero. If not using +mkin transformations or a custom covariance model, the default +initialisation of \link[saemix:saemixModel]{saemix::saemixModel} is used for omega.init.} + \item{covariates}{A data frame with covariate data for use in 'covariate_models', with dataset names as row names.} @@ -96,6 +108,8 @@ in the 'covariates' data frame.} which there should be no variability over the groups. Only used if the covariance model is not explicitly specified.} +\item{error.init}{Will be passed to \code{\link[saemix:saemixModel]{saemix::saemixModel()}}.} + \item{nbiter.saemix}{Convenience option to increase the number of iterations} 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 fe38865d..6015aed8 100644 --- a/tests/testthat/_snaps/multistart/llhist-for-biphasic-saemix-fit.svg +++ b/tests/testthat/_snaps/multistart/llhist-for-biphasic-saemix-fit.svg @@ -21,24 +21,28 @@ Frequency of log likelihoods - - - + + + + --1185 --1180 --1175 --1170 --1165 +-1149.5 +-1149.4 +-1149.3 +-1149.2 +-1149.1 +-1149.0 - - + + + 0 -1 -2 -3 +1 +2 +3 +4 @@ -46,11 +50,12 @@ - - - - - + + + + + + original fit diff --git a/tests/testthat/_snaps/multistart/parplot-for-biphasic-saemix-fit.svg b/tests/testthat/_snaps/multistart/parplot-for-biphasic-saemix-fit.svg index 75519b07..c0332fd5 100644 --- a/tests/testthat/_snaps/multistart/parplot-for-biphasic-saemix-fit.svg +++ b/tests/testthat/_snaps/multistart/parplot-for-biphasic-saemix-fit.svg @@ -25,105 +25,109 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -163,26 +167,26 @@ - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + diff --git a/tests/testthat/_snaps/plot/mixed-model-fit-for-saem-object-with-mkin-transformations.svg b/tests/testthat/_snaps/plot/mixed-model-fit-for-saem-object-with-mkin-transformations.svg index 375ab089..69fa6a4d 100644 --- a/tests/testthat/_snaps/plot/mixed-model-fit-for-saem-object-with-mkin-transformations.svg +++ b/tests/testthat/_snaps/plot/mixed-model-fit-for-saem-object-with-mkin-transformations.svg @@ -156,7 +156,7 @@ - + @@ -175,7 +175,7 @@ - + @@ -212,7 +212,7 @@ - + @@ -249,7 +249,7 @@ - + @@ -268,7 +268,7 @@ - + @@ -287,7 +287,7 @@ - + @@ -342,7 +342,7 @@ - + @@ -415,7 +415,7 @@ - + @@ -470,7 +470,7 @@ - + @@ -525,7 +525,7 @@ - + @@ -562,7 +562,7 @@ - + @@ -617,7 +617,7 @@ - + @@ -672,7 +672,7 @@ - + @@ -709,7 +709,7 @@ - + @@ -728,8 +728,8 @@ - - + + @@ -739,34 +739,34 @@ - + - - - - - + + + + + 0 -20 -40 -60 -80 -100 - - - - +20 +40 +60 +80 +100 + + + + - - - --3 --2 --1 + + + +-3 +-2 +-1 0 -1 -2 -3 +1 +2 +3 @@ -780,582 +780,582 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -1419,7 +1419,7 @@ - + @@ -1436,7 +1436,7 @@ - + @@ -1469,7 +1469,7 @@ - + @@ -1502,7 +1502,7 @@ - + @@ -1519,7 +1519,7 @@ - + @@ -1536,7 +1536,7 @@ - + @@ -1585,7 +1585,7 @@ - + @@ -1650,7 +1650,7 @@ - + @@ -1699,7 +1699,7 @@ - + @@ -1748,7 +1748,7 @@ - + @@ -1781,7 +1781,7 @@ - + @@ -1830,7 +1830,7 @@ - + @@ -1879,7 +1879,7 @@ - + @@ -1912,7 +1912,7 @@ - + @@ -1929,8 +1929,8 @@ - - + + @@ -1940,32 +1940,32 @@ - + - - - - + + + + 0 -10 -20 -30 -40 - - - - +10 +20 +30 +40 + + + + - - - --3 --2 --1 + + + +-3 +-2 +-1 0 -1 -2 -3 +1 +2 +3 @@ -1979,518 +1979,518 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/testthat/anova_sfo_saem.txt b/tests/testthat/anova_sfo_saem.txt index 4a21e81f..9e4bf71f 100644 --- a/tests/testthat/anova_sfo_saem.txt +++ b/tests/testthat/anova_sfo_saem.txt @@ -4,4 +4,4 @@ Data: 262 observations of 1 variable(s) grouped in 15 datasets sfo_saem_1_reduced 5 1310 1313 -650 sfo_saem_1_reduced_mkin 5 1310 1313 -650 0 0 sfo_saem_1 6 1312 1316 -650 0 1 1 -sfo_saem_1_mkin 6 1310 1315 -649 1 0 +sfo_saem_1_mkin 6 1312 1316 -650 0 0 diff --git a/tests/testthat/print_dfop_saemix_1.txt b/tests/testthat/print_dfop_saemix_1.txt index f427b3e6..f6fda37c 100644 --- a/tests/testthat/print_dfop_saemix_1.txt +++ b/tests/testthat/print_dfop_saemix_1.txt @@ -13,12 +13,12 @@ Likelihood computed by importance sampling Fitted parameters: estimate lower upper -parent_0 100.04 98.89 101.20 -log_k1 -4.12 -4.24 -4.00 -log_k2 -2.67 -2.90 -2.44 -g_qlogis 0.43 0.22 0.64 -a.1 0.92 0.67 1.16 +parent_0 100.09 98.94 101.25 +log_k1 -2.68 -2.91 -2.45 +log_k2 -4.12 -4.24 -4.00 +g_qlogis -0.41 -0.63 -0.20 +a.1 0.91 0.67 1.15 b.1 0.05 0.04 0.06 -SD.log_k1 0.22 0.14 0.30 -SD.log_k2 0.36 0.21 0.51 -SD.g_qlogis 0.14 -0.11 0.39 +SD.log_k1 0.36 0.21 0.50 +SD.log_k2 0.22 0.13 0.30 +SD.g_qlogis 0.15 -0.09 0.40 diff --git a/tests/testthat/summary_hfit_sfo_tc.txt b/tests/testthat/summary_hfit_sfo_tc.txt index 49765187..41743091 100644 --- a/tests/testthat/summary_hfit_sfo_tc.txt +++ b/tests/testthat/summary_hfit_sfo_tc.txt @@ -32,9 +32,9 @@ Likelihood computed by importance sampling Optimised parameters: est. lower upper -parent_0 101.01 99.57 102.45 +parent_0 101.02 99.58 102.46 log_k_parent -3.32 -3.53 -3.11 -a.1 0.90 0.64 1.17 +a.1 0.91 0.64 1.17 b.1 0.05 0.04 0.06 SD.log_k_parent 0.27 0.11 0.42 @@ -48,7 +48,7 @@ SD.log_k_parent 0.3 0.1 0.4 Variance model: est. lower upper -a.1 0.90 0.64 1.17 +a.1 0.91 0.64 1.17 b.1 0.05 0.04 0.06 Backtransformed parameters: diff --git a/tests/testthat/test_multistart.R b/tests/testthat/test_multistart.R index 56eb140c..c1a10d10 100644 --- a/tests/testthat/test_multistart.R +++ b/tests/testthat/test_multistart.R @@ -21,7 +21,10 @@ test_that("multistart works for saem.mmkin models", { anova_biphasic <- anova(saem_biphasic_m, best(saem_biphasic_m_multi)) - expect_true(anova_biphasic[2, "AIC"] < anova_biphasic[1, "AIC"]) + # With the new starting parameters we do not improve + # with multistart any more + expect_equal(anova_biphasic[2, "AIC"], anova_biphasic[1, "AIC"], + tolerance = 1e-4) skip_on_travis() # Plots are platform dependent llhist_sfo <- function() llhist(saem_sfo_s_multi) diff --git a/tests/testthat/test_saemix_parent.R b/tests/testthat/test_saemix_parent.R index 79c5fa69..20889c6c 100644 --- a/tests/testthat/test_saemix_parent.R +++ b/tests/testthat/test_saemix_parent.R @@ -73,40 +73,36 @@ test_that("Parent fits using saemix are correctly implemented", { # FOMC mmkin_fomc_1 <- mmkin("FOMC", ds_fomc, quiet = TRUE, error_model = "tc", cores = n_cores) fomc_saem_1 <- saem(mmkin_fomc_1, quiet = TRUE, transformations = "saemix", no_random_effect = "parent_0") - fomc_saem_2 <- update(fomc_saem_1, transformations = "mkin") - ci_fomc_s1 <- summary(fomc_saem_1)$confint_back fomc_pop <- as.numeric(fomc_pop) + ci_fomc_s1 <- summary(fomc_saem_1)$confint_back expect_true(all(ci_fomc_s1[, "lower"] < fomc_pop)) expect_true(all(ci_fomc_s1[, "upper"] > fomc_pop)) - expect_equal(endpoints(fomc_saem_1), endpoints(fomc_saem_2), tol = 0.01) - mmkin_fomc_2 <- update(mmkin_fomc_1, state.ini = 100, fixed_initials = "parent") - fomc_saem_2 <- saem(mmkin_fomc_2, quiet = TRUE, transformations = "mkin") + fomc_saem_2 <- update(fomc_saem_1, transformations = "mkin") ci_fomc_s2 <- summary(fomc_saem_2)$confint_back + expect_true(all(ci_fomc_s2[, "lower"] < fomc_pop)) + expect_true(all(ci_fomc_s2[, "upper"] > fomc_pop)) - expect_true(all(ci_fomc_s2[, "lower"] < fomc_pop[2:3])) - expect_true(all(ci_fomc_s2[, "upper"] > fomc_pop[2:3])) + expect_equal(endpoints(fomc_saem_1), endpoints(fomc_saem_2), tol = 0.01) # DFOP dfop_saemix_2 <- saem(mmkin_dfop_1, quiet = TRUE, transformations = "saemix", no_random_effect = "parent_0") - s_dfop_s1 <- summary(dfop_saemix_1) - s_dfop_s2 <- summary(dfop_saemix_2) + s_dfop_s1 <- summary(dfop_saemix_1) # mkin transformations + s_dfop_s2 <- summary(dfop_saemix_2) # saemix transformations s_dfop_n <- summary(dfop_nlme_1) dfop_pop <- as.numeric(dfop_pop) - # When using DFOP with mkin transformations, k1 and k2 are sometimes swapped - swap_k1_k2 <- function(p) c(p[1], p[3], p[2], 1 - p[4]) - expect_true(all(s_dfop_s1$confint_back[, "lower"] < swap_k1_k2(dfop_pop))) - expect_true(all(s_dfop_s1$confint_back[, "upper"] > swap_k1_k2(dfop_pop))) + expect_true(all(s_dfop_s1$confint_back[, "lower"] < dfop_pop)) + expect_true(all(s_dfop_s1$confint_back[, "upper"] > dfop_pop)) expect_true(all(s_dfop_s2$confint_back[, "lower"] < dfop_pop)) expect_true(all(s_dfop_s2$confint_back[, "upper"] > dfop_pop)) - # We get < 20% deviations with transformations made in mkin (need to swap k1 and k2) - rel_diff_1 <- (swap_k1_k2(s_dfop_s1$confint_back[, "est."]) - dfop_pop) / dfop_pop + # We get < 20% deviations with transformations made in mkin + rel_diff_1 <- (s_dfop_s1$confint_back[, "est."] - dfop_pop) / dfop_pop expect_true(all(rel_diff_1 < 0.20)) # We get < 20% deviations with transformations made in saemix @@ -123,27 +119,21 @@ test_that("Parent fits using saemix are correctly implemented", { transformations = "saemix") expect_equal( log(endpoints(dfop_saemix_1)$distimes[1:2]), - log(endpoints(sforb_saemix_1)$distimes[1:2]), tolerance = 0.03) + log(endpoints(sforb_saemix_1)$distimes[1:2]), tolerance = 0.01) expect_equal( log(endpoints(sforb_saemix_1)$distimes[1:2]), log(endpoints(sforb_saemix_2)$distimes[1:2]), tolerance = 0.01) mmkin_hs_1 <- mmkin("HS", ds_hs, quiet = TRUE, error_model = "const", cores = n_cores) - hs_saem_1 <- saem(mmkin_hs_1, quiet = TRUE) - hs_saem_2 <- saem(mmkin_hs_1, quiet = TRUE, transformations = "mkin") + hs_saem_1 <- saem(mmkin_hs_1, quiet = TRUE, no_random_effect = "parent_0") + hs_saem_2 <- saem(mmkin_hs_1, quiet = TRUE, transformations = "mkin", + no_random_effect = "parent_0") expect_equal(endpoints(hs_saem_1), endpoints(hs_saem_2), tol = 0.01) ci_hs_s1 <- summary(hs_saem_1)$confint_back hs_pop <- as.numeric(hs_pop) #expect_true(all(ci_hs_s1[, "lower"] < hs_pop)) # k1 is overestimated expect_true(all(ci_hs_s1[, "upper"] > hs_pop)) - - mmkin_hs_2 <- update(mmkin_hs_1, state.ini = 100, fixed_initials = "parent") - hs_saem_3 <- saem(mmkin_hs_2, quiet = TRUE) - ci_hs_s3 <- summary(hs_saem_3)$confint_back - - #expect_true(all(ci_hs_s3[, "lower"] < hs_pop[2:4])) # k1 again overestimated - expect_true(all(ci_hs_s3[, "upper"] > hs_pop[2:4])) }) test_that("We can also use mkin solution methods for saem", { -- cgit v1.2.1