diff options
| -rw-r--r-- | NEWS.md | 2 | ||||
| -rw-r--r-- | R/mhmkin.R | 91 | ||||
| -rw-r--r-- | log/check.log | 8 | ||||
| -rw-r--r-- | log/test.log | 48 | ||||
| -rw-r--r-- | man/mhmkin.Rd | 50 | ||||
| -rw-r--r-- | tests/testthat/illparms_hfits_synth.txt | 10 | ||||
| -rw-r--r-- | tests/testthat/print_fits_synth_const.txt | 4 | ||||
| -rw-r--r-- | tests/testthat/summary_hfit_sfo_tc.txt | 34 | ||||
| -rw-r--r-- | tests/testthat/test_mhmkin.R | 34 |
9 files changed, 187 insertions, 94 deletions
@@ -1,6 +1,6 @@ # mkin 1.2.2 -- 'R/mhmkin.R': Allow an 'illparms.mhmkin' object as value of the argument 'no_random_effects', making it possible to exclude random effects that were ill-defined in simpler variants of the set of degradation models. Remove the possibility to exclude random effects based on separate fits, as it did not work well. +- 'R/mhmkin.R': Allow an 'illparms.mhmkin' object or a list with suitable dimensions as value of the argument 'no_random_effects', making it possible to exclude random effects that were ill-defined in simpler variants of the set of degradation models. Remove the possibility to exclude random effects based on separate fits, as it did not work well. - 'R/summary.saem.mmkin.R': List all initial parameter values in the summary, including random effects and error model parameters @@ -14,11 +14,12 @@ #' supported #' @param no_random_effect Default is NULL and will be passed to [saem]. If a #' character vector is supplied, it will be passed to all calls to [saem], -#' regardless if the corresponding parameter is in the model. Alternatively, -#' an object of class [illparms.mhmkin] can be specified. This has to have -#' the same dimensions as the return object of the current call. In this way, -#' ill-defined parameters found in corresponding simpler versions of the -#' degradation model can be specified. +#' which will exclude random effects for all matching parameters. Alternatively, +#' a list of character vectors or an object of class [illparms.mhmkin] can be +#' specified. They have to have the same dimensions that the return object of +#' the current call will have, i.e. the number of rows must match the number +#' of degradation models in the mmkin object(s), and the number of columns must +#' match the number of error models used in the mmkin object(s). #' @param algorithm The algorithm to be used for fitting (currently not used) #' @param \dots Further arguments that will be passed to the nonlinear mixed-effects #' model fitting function. @@ -50,6 +51,42 @@ mhmkin.mmkin <- function(objects, ...) { #' @export #' @rdname mhmkin +#' @examples +#' \dontrun{ +#' # We start with separate evaluations of all the first six datasets with two +#' # degradation models and two error models +#' f_sep_const <- mmkin(c("SFO", "FOMC"), ds_fomc[1:6], cores = 2, quiet = TRUE) +#' f_sep_tc <- update(f_sep_const, error_model = "tc") +#' # The mhmkin function sets up hierarchical degradation models aka +#' # nonlinear mixed-effects models for all four combinations, specifying +#' # uncorrelated random effects for all degradation parameters +#' f_saem_1 <- mhmkin(list(f_sep_const, f_sep_tc), cores = 2) +#' status(f_saem_1) +#' # The 'illparms' function shows that in all hierarchical fits, at least +#' # one random effect is ill-defined (the confidence interval for the +#' # random effect expressed as standard deviation includes zero) +#' illparms(f_saem_1) +#' # Therefore we repeat the fits, excluding the ill-defined random effects +#' f_saem_2 <- update(f_saem_1, no_random_effect = illparms(f_saem_1)) +#' status(f_saem_2) +#' illparms(f_saem_2) +#' # Model comparisons show that FOMC with two-component error is preferable, +#' # and confirms our reduction of the default parameter model +#' anova(f_saem_1) +#' anova(f_saem_2) +#' # The convergence plot for the selected model looks fine +#' saemix::plot(f_saem_2[["FOMC", "tc"]]$so, plot.type = "convergence") +#' # The plot of predictions versus data shows that we have a pretty data-rich +#' # situation with homogeneous distribution of residuals, because we used the +#' # same degradation model, error model and parameter distribution model that +#' # was used in the data generation. +#' plot(f_saem_2[["FOMC", "tc"]]) +#' # We can specify the same parameter model reductions manually +#' no_ranef <- list("parent_0", "log_beta", "parent_0", c("parent_0", "log_beta")) +#' dim(no_ranef) <- c(2, 2) +#' f_saem_2m <- update(f_saem_1, no_random_effect = no_ranef) +#' anova(f_saem_2m) +#' } mhmkin.list <- function(objects, backend = "saemix", algorithm = "saem", no_random_effect = NULL, ..., @@ -97,25 +134,38 @@ mhmkin.list <- function(objects, backend = "saemix", algorithm = "saem", dimnames(fit_indices) <- list(degradation = names(deg_models), error = error_models) - fit_function <- function(fit_index) { - w <- which(fit_indices == fit_index, arr.ind = TRUE) - deg_index <- w[1] - error_index <- w[2] - mmkin_row <- objects[[error_index]][deg_index, ] + if (is.null(no_random_effect) || length(dim(no_random_effect)) == 1) { + no_ranef <- rep(list(no_random_effect), n.fits) + dim(no_ranef) <- dim(fit_indices) + } else { + if (!identical(dim(no_random_effect), dim(fit_indices))) { + stop("Dimensions of argument 'no_random_effect' are not suitable") + } if (is(no_random_effect, "illparms.mhmkin")) { - if (identical(dim(no_random_effect), dim(fit_indices))) { - no_ranef_split <- strsplit(no_random_effect[[fit_index]], ", ") - no_ranef <- sapply(no_ranef_split, function(x) { - gsub("sd\\((.*)\\)", "\\1", x) + no_ranef_dim <- dim(no_random_effect) + no_ranef <- lapply(no_random_effect, function(x) { + no_ranef_split <- strsplit(x, ", ") + ret <- sapply(no_ranef_split, function(y) { + gsub("sd\\((.*)\\)", "\\1", y) }) - } else { - stop("Dimensions of illparms.mhmkin object given as 'no_random_effect' are not suitable") - } + return(ret) + }) + dim(no_ranef) <- no_ranef_dim } else { no_ranef <- no_random_effect } + } + + fit_function <- function(fit_index) { + w <- which(fit_indices == fit_index, arr.ind = TRUE) + deg_index <- w[1] + error_index <- w[2] + mmkin_row <- objects[[error_index]][deg_index, ] res <- try(do.call(backend_function, - args = c(list(mmkin_row), dot_args, list(no_random_effect = no_ranef)))) + args = c( + list(mmkin_row), + dot_args, + list(no_random_effect = no_ranef[[deg_index, error_index]])))) return(res) } @@ -145,15 +195,16 @@ mhmkin.list <- function(objects, backend = "saemix", algorithm = "saem", #' @param j Column index selecting the fits to specific datasets #' @param drop If FALSE, the method always returns an mhmkin object, otherwise #' either a list of fit objects or a single fit object. -#' @return An object of class \code{\link{mhmkin}}. +#' @return An object inheriting from \code{\link{mhmkin}}. #' @rdname mhmkin #' @export `[.mhmkin` <- function(x, i, j, ..., drop = FALSE) { + original_class <- class(x) class(x) <- NULL x_sub <- x[i, j, drop = drop] if (!drop) { - class(x_sub) <- "mhmkin" + class(x_sub) <- original_class } return(x_sub) } diff --git a/log/check.log b/log/check.log index 31fc31eb..42365918 100644 --- a/log/check.log +++ b/log/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 ‘1.2.1’ +* this is package ‘mkin’ version ‘1.2.2’ * package encoding: UTF-8 * checking CRAN incoming feasibility ... Note_to_CRAN_maintainers Maintainer: ‘Johannes Ranke <johannes.ranke@jrwb.de>’ @@ -18,7 +18,7 @@ Maintainer: ‘Johannes Ranke <johannes.ranke@jrwb.de>’ * checking for portable file names ... OK * checking for sufficient/correct file permissions ... OK * checking serialization versions ... OK -* checking whether package ‘mkin’ can be installed ... OK +* checking whether package ‘mkin’ can be installed ... [11s/11s] OK * checking installed package size ... OK * checking package directory ... OK * checking for future file timestamps ... OK @@ -41,7 +41,7 @@ Maintainer: ‘Johannes Ranke <johannes.ranke@jrwb.de>’ * checking S3 generic/method consistency ... OK * checking replacement functions ... OK * checking foreign function calls ... OK -* checking R code for possible problems ... [17s/17s] OK +* checking R code for possible problems ... [19s/19s] OK * checking Rd files ... OK * checking Rd metadata ... OK * checking Rd line widths ... OK @@ -57,7 +57,7 @@ Maintainer: ‘Johannes Ranke <johannes.ranke@jrwb.de>’ * checking data for ASCII and uncompressed saves ... OK * checking installed files from ‘inst/doc’ ... OK * checking files in ‘vignettes’ ... OK -* checking examples ... [20s/20s] OK +* checking examples ... [24s/24s] OK * checking for unstated dependencies in ‘tests’ ... OK * checking tests ... SKIPPED * checking for unstated dependencies in vignettes ... OK diff --git a/log/test.log b/log/test.log index e17ecc1f..84fa49b9 100644 --- a/log/test.log +++ b/log/test.log @@ -1,58 +1,58 @@ ℹ Testing mkin ✔ | F W S OK | Context ✔ | 5 | AIC calculation -✔ | 5 | Analytical solutions for coupled models [3.3s] +✔ | 5 | Analytical solutions for coupled models [4.2s] ✔ | 5 | Calculation of Akaike weights ✔ | 3 | Export dataset for reading into CAKE -✔ | 12 | Confidence intervals and p-values [1.1s] -✔ | 1 12 | Dimethenamid data from 2018 [32.2s] +✔ | 12 | Confidence intervals and p-values [1.2s] +✔ | 1 12 | Dimethenamid data from 2018 [42.0s] ──────────────────────────────────────────────────────────────────────────────── 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 [4.9s] +✔ | 14 | Error model fitting [6.5s] ✔ | 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.3s] -✔ | 1 | Fitting the logistic model [0.2s] -✔ | 8 | Batch fitting and diagnosing hierarchical kinetic models [14.5s] -✔ | 1 11 | Nonlinear mixed-effects models [13.1s] +✔ | 4 | Calculation of FOCUS chi2 error levels [0.7s] +✔ | 14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [1.1s] +✔ | 4 | Test fitting the decline of metabolites from their maximum [0.5s] +✔ | 1 | Fitting the logistic model [0.3s] +✔ | 10 | Batch fitting and diagnosing hierarchical kinetic models [54.1s] +✔ | 1 11 | Nonlinear mixed-effects models [14.3s] ──────────────────────────────────────────────────────────────────────────────── Skip ('test_mixed.R:78'): 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.4s] -✔ | 3 | mkinfit features [0.7s] -✔ | 8 | mkinmod model generation and printing [0.2s] +✔ | 10 | Special cases of mkinfit calls [0.8s] +✔ | 3 | mkinfit features [0.9s] +✔ | 8 | mkinmod model generation and printing [0.3s] ✔ | 3 | Model predictions with mkinpredict [0.3s] -✔ | 12 | Multistart method for saem.mmkin models [50.1s] -✔ | 16 | Evaluations according to 2015 NAFTA guidance [2.2s] -✔ | 9 | Nonlinear mixed-effects models with nlme [8.7s] -✔ | 15 | Plotting [10.2s] +✔ | 12 | Multistart method for saem.mmkin models [80.1s] +✔ | 16 | Evaluations according to 2015 NAFTA guidance [2.8s] +✔ | 9 | Nonlinear mixed-effects models with nlme [11.4s] +✔ | 15 | Plotting [12.1s] ✔ | 4 | Residuals extracted from mkinfit models -✔ | 1 36 | saemix parent models [103.8s] +✔ | 1 36 | saemix parent models [85.9s] ──────────────────────────────────────────────────────────────────────────────── 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] +✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.9s] ✔ | 11 | Processing of residue series -✔ | 10 | Fitting the SFORB model [3.8s] +✔ | 10 | Fitting the SFORB model [4.6s] ✔ | 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.1s] +✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.9s] +✔ | 9 | Hypothesis tests [11.0s] ✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s] ══ Results ═════════════════════════════════════════════════════════════════════ -Duration: 266.0 s +Duration: 342.6 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 270 ] Error while shutting down parallel: unable to terminate some child processes diff --git a/man/mhmkin.Rd b/man/mhmkin.Rd index 4230e44f..c77f4289 100644 --- a/man/mhmkin.Rd +++ b/man/mhmkin.Rd @@ -43,11 +43,12 @@ supported} \item{no_random_effect}{Default is NULL and will be passed to \link{saem}. If a character vector is supplied, it will be passed to all calls to \link{saem}, -regardless if the corresponding parameter is in the model. Alternatively, -an object of class \link{illparms.mhmkin} can be specified. This has to have -the same dimensions as the return object of the current call. In this way, -ill-defined parameters found in corresponding simpler versions of the -degradation model can be specified.} +which will exclude random effects for all matching parameters. Alternatively, +a list of character vectors or an object of class \link{illparms.mhmkin} can be +specified. They have to have the same dimensions that the return object of +the current call will have, i.e. the number of rows must match the number +of degradation models in the mmkin object(s), and the number of columns must +match the number of error models used in the mmkin object(s).} \item{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 @@ -74,7 +75,7 @@ be indexed using the degradation model names for the first index (row index) and the error model names for the second index (column index), with class attribute 'mhmkin'. -An object of class \code{\link{mhmkin}}. +An object inheriting from \code{\link{mhmkin}}. } \description{ The name of the methods expresses that (\strong{m}ultiple) \strong{h}ierarchichal @@ -82,6 +83,43 @@ The name of the methods expresses that (\strong{m}ultiple) \strong{h}ierarchicha fitted. Our kinetic models are nonlinear, so we can use various nonlinear mixed-effects model fitting functions. } +\examples{ +\dontrun{ +# We start with separate evaluations of all the first six datasets with two +# degradation models and two error models +f_sep_const <- mmkin(c("SFO", "FOMC"), ds_fomc[1:6], cores = 2, quiet = TRUE) +f_sep_tc <- update(f_sep_const, error_model = "tc") +# The mhmkin function sets up hierarchical degradation models aka +# nonlinear mixed-effects models for all four combinations, specifying +# uncorrelated random effects for all degradation parameters +f_saem_1 <- mhmkin(list(f_sep_const, f_sep_tc), cores = 2) +status(f_saem_1) +# The 'illparms' function shows that in all hierarchical fits, at least +# one random effect is ill-defined (the confidence interval for the +# random effect expressed as standard deviation includes zero) +illparms(f_saem_1) +# Therefore we repeat the fits, excluding the ill-defined random effects +f_saem_2 <- update(f_saem_1, no_random_effect = illparms(f_saem_1)) +status(f_saem_2) +illparms(f_saem_2) +# Model comparisons show that FOMC with two-component error is preferable, +# and confirms our reduction of the default parameter model +anova(f_saem_1) +anova(f_saem_2) +# The convergence plot for the selected model looks fine +saemix::plot(f_saem_2[["FOMC", "tc"]]$so, plot.type = "convergence") +# The plot of predictions versus data shows that we have a pretty data-rich +# situation with homogeneous distribution of residuals, because we used the +# same degradation model, error model and parameter distribution model that +# was used in the data generation. +plot(f_saem_2[["FOMC", "tc"]]) +# We can specify the same parameter model reductions manually +no_ranef <- list("parent_0", "log_beta", "parent_0", c("parent_0", "log_beta")) +dim(no_ranef) <- c(2, 2) +f_saem_2m <- update(f_saem_1, no_random_effect = no_ranef) +anova(f_saem_2m) +} +} \seealso{ \code{\link{[.mhmkin}} for subsetting \link{mhmkin} objects } diff --git a/tests/testthat/illparms_hfits_synth.txt b/tests/testthat/illparms_hfits_synth.txt index affd1318..7a69645b 100644 --- a/tests/testthat/illparms_hfits_synth.txt +++ b/tests/testthat/illparms_hfits_synth.txt @@ -1,8 +1,4 @@ error -degradation const - SFO - FOMC sd(log_alpha), sd(log_beta) - error -degradation tc - SFO sd(parent_0) - FOMC sd(parent_0), sd(log_alpha), sd(log_beta) +degradation const tc + SFO sd(parent_0) sd(parent_0) + FOMC sd(log_beta) sd(parent_0), sd(log_beta) diff --git a/tests/testthat/print_fits_synth_const.txt b/tests/testthat/print_fits_synth_const.txt index b4bbe6ca..5d076d3d 100644 --- a/tests/testthat/print_fits_synth_const.txt +++ b/tests/testthat/print_fits_synth_const.txt @@ -4,8 +4,6 @@ Status of individual fits: dataset model 1 2 3 4 5 6 SFO OK OK OK OK OK OK - FOMC C OK OK OK OK C + FOMC OK OK OK OK OK OK -C: Optimisation did not converge: -false convergence (8) OK: No warnings diff --git a/tests/testthat/summary_hfit_sfo_tc.txt b/tests/testthat/summary_hfit_sfo_tc.txt index 0a61f75f..0618c715 100644 --- a/tests/testthat/summary_hfit_sfo_tc.txt +++ b/tests/testthat/summary_hfit_sfo_tc.txt @@ -8,7 +8,7 @@ Equations: d_parent/dt = - k_parent * parent Data: -104 observations of 1 variable(s) grouped in 6 datasets +95 observations of 1 variable(s) grouped in 6 datasets Model predictions using solution type analytical @@ -19,7 +19,7 @@ Variance model: Two-component variance function Starting values for degradation parameters: parent_0 log_k_parent - 101 -3 + 94 -2 Fixed degradation parameter values: None @@ -27,7 +27,7 @@ None Starting values for random effects (square root of initial entries in omega): parent_0 log_k_parent parent_0 4 0.0 -log_k_parent 0 0.4 +log_k_parent 0 0.7 Starting values for error model parameters: a.1 b.1 @@ -37,15 +37,15 @@ Results: Likelihood computed by importance sampling AIC BIC logLik - 524 523 -257 + 542 541 -266 Optimised parameters: - est. lower upper -parent_0 100.68 99.27 102.08 -log_k_parent -3.38 -3.55 -3.21 -a.1 0.87 0.59 1.14 -b.1 0.05 0.04 0.06 -SD.log_k_parent 0.21 0.09 0.33 + est. lower upper +parent_0 92.52 89.11 95.9 +log_k_parent -1.66 -2.07 -1.3 +a.1 2.03 1.60 2.5 +b.1 0.09 0.07 0.1 +SD.log_k_parent 0.51 0.22 0.8 Correlation: pr_0 @@ -53,18 +53,18 @@ log_k_parent 0.1 Random effects: est. lower upper -SD.log_k_parent 0.2 0.09 0.3 +SD.log_k_parent 0.5 0.2 0.8 Variance model: est. lower upper -a.1 0.87 0.59 1.14 -b.1 0.05 0.04 0.06 +a.1 2.03 1.60 2.5 +b.1 0.09 0.07 0.1 Backtransformed parameters: - est. lower upper -parent_0 1e+02 99.27 1e+02 -k_parent 3e-02 0.03 4e-02 + est. lower upper +parent_0 92.5 89.1 95.9 +k_parent 0.2 0.1 0.3 Estimated disappearance times: DT50 DT90 -parent 20 68 +parent 4 12 diff --git a/tests/testthat/test_mhmkin.R b/tests/testthat/test_mhmkin.R index e2339f28..da063326 100644 --- a/tests/testthat/test_mhmkin.R +++ b/tests/testthat/test_mhmkin.R @@ -3,8 +3,11 @@ context("Batch fitting and diagnosing hierarchical kinetic models") test_that("Multiple hierarchical kinetic models can be fitted and diagnosed", { skip_on_cran() - fits_synth_const <- suppressWarnings( - mmkin(c("SFO", "FOMC"), ds_sfo[1:6], cores = n_cores, quiet = TRUE)) + fits_synth_const <- mmkin(c("SFO", "FOMC"), ds_fomc[1:6], cores = n_cores, quiet = TRUE) + + expect_known_output( + print(fits_synth_const), + "print_fits_synth_const.txt") fits_synth_tc <- suppressWarnings( update(fits_synth_const, error_model = "tc")) @@ -19,8 +22,8 @@ test_that("Multiple hierarchical kinetic models can be fitted and diagnosed", { print(illparms(hfits)), "illparms_hfits_synth.txt") - expect_equal(which.min(AIC(hfits)), 3) - expect_equal(which.min(BIC(hfits)), 3) + expect_equal(which.min(AIC(hfits)), 4) + expect_equal(which.min(BIC(hfits)), 4) hfit_sfo_tc <- update(hfits[["SFO", "tc"]], covariance.model = diag(c(0, 1))) @@ -38,12 +41,19 @@ test_that("Multiple hierarchical kinetic models can be fitted and diagnosed", { expect_known_output(print(test_summary, digits = 1), "summary_hfit_sfo_tc.txt") - # It depends on the platform exactly which of the datasets fail to converge - # with FOMC, because they were generated to be SFO - skip_on_travis() - - expect_known_output( - print(fits_synth_const), - "print_fits_synth_const.txt") - + hfits_sfo_reduced <- update(hfits, + no_random_effect = illparms(hfits)) + expect_equal( + as.character(illparms(hfits_sfo_reduced)), + rep("", 4)) + + # We can also manually set up an object specifying random effects to be + # excluded. Entries in the inital list have to be by column + no_ranef <- list("parent_0", "log_beta", "parent_0", c("parent_0", "log_beta")) + dim(no_ranef) <- c(2, 2) + + hfits_sfo_reduced_2 <- update(hfits, + no_random_effect = no_ranef) + expect_equivalent(round(anova(hfits_sfo_reduced), 0), + round(anova(hfits_sfo_reduced_2), 0)) }) |